Update functions to output geojson for QGIS

src/cli.py

@@ -29,7 +29,7 @@ if __name__ == "__main__":
     cli.add_command(parse_mat.create_tides_csv)
     cli.add_command(parse_mat.create_profile_features)
     # cli.add_command(profile_features.create_profile_features)
-    cli.add_command(csv_to_shp.sites_csv_to_shp)
+    cli.add_command(csv_to_shp.sites_csv_to_geojson)
     cli.add_command(forecast_twl.create_twl_forecast)
     cli.add_command(forecasted_storm_impacts.create_forecasted_impacts)
     cli.add_command(observed_storm_impacts.create_observed_impacts)

src/data/csv_to_shp.py

@@ -2,36 +2,219 @@
 Converts .csv files to .shape files
 """
 import os
-
+import numpy.ma as ma
 import click
 import fiona
 import numpy as np
 import pandas as pd
 from fiona.crs import from_epsg
 from shapely.geometry import Point, mapping, LineString
-
+from collections import OrderedDict
 from data.parse_shp import convert_coord_systems
 from utils import setup_logging
 
 logger = setup_logging()
 
 
+
+def R_high_to_geojson(sites_csv, profiles_csv, impacts_csv, output_geojson):
+
+    sites_csv = './data/interim/sites.csv'
+    profiles_csv = './data/interim/profiles.csv'
+    impacts_csv = './data/interim/impacts_forecasted_mean_slope_sto06.csv'
+    output_geojson = './data/interim/R_high_forecasted_mean_slope_sto06.geojson'
+
+    df_sites = pd.read_csv(sites_csv, index_col=[0])
+    df_profiles = pd.read_csv(profiles_csv, index_col=[0,1,2])
+    df_impacts = pd.read_csv(impacts_csv, index_col=[0])
+
+    # Create geojson file
+    schema = {
+        'geometry': 'Point',
+        'properties': OrderedDict([
+            ('beach', 'str'),
+            ('site_id', 'str'),
+            ('elevation', 'float'),
+        ])
+    }
+
+    with fiona.open(output_geojson, 'w', driver='GeoJSON', crs=from_epsg(4326), schema=schema) as output:
+        for index, row in df_impacts.iterrows():
+
+            site_id = index
+            beach = index[:-4]
+
+            # Find lat/lon of R_high position
+            R_high_z = row['R_high']
+
+
+            # Get poststorm profile (or should this be prestorm?)
+            df_profile = df_profiles.query('site_id=="{}" & profile_type=="prestorm"'.format(index))
+            int_x = crossings(df_profile.index.get_level_values('x').tolist(),
+                              df_profile.z.tolist(),
+                              R_high_z)
+
+            # Take most landward interesection. Continue to next site if there is no intersection
+            try:
+                int_x = max(int_x)
+            except:
+                continue
+
+            # Get lat/lon on intercept position
+            site = df_sites.loc[site_id]
+            center_profile_x = site["profile_x_lat_lon"]
+            orientation = site["orientation"]
+            center_lat_lon = Point(site['lon'], site['lat'])  # Get lat/lon of center of profile
+            center_xy = convert_coord_systems(center_lat_lon)
+            center_x, center_y = center_xy.xy
+
+            # Calculate xy position of point and convert to lat/lon
+            point_x = center_x + (center_profile_x - int_x) * np.cos(np.deg2rad(orientation))
+            point_y = center_y + (center_profile_x - int_x) * np.sin(np.deg2rad(orientation))
+            point_xy = Point(point_x, point_y)
+            point_lat_lon = convert_coord_systems(point_xy,
+                                                  in_coord_system="EPSG:28356",
+                                                  out_coord_system="EPSG:4326")
+
+            prop = OrderedDict([
+                ('beach', beach),
+                ('site_id', site_id),
+                ('elevation', R_high_z),
+            ])
+            output.write({"geometry": mapping(point_lat_lon), "properties": prop})
+
+
+
+def crossings(profile_x, profile_z, constant_z):
+    """
+    Finds the x coordinate of a z elevation for a beach profile. Much faster than using shapely to calculate
+    intersections since we are only interested in intersections of a constant value. Will return multiple
+    intersections if found. Used in calculating beach slope.
+    Adapted from https://stackoverflow.com/a/34745789
+    :param profile_x: List of x coordinates for the beach profile section
+    :param profile_z: List of z coordinates for the beach profile section
+    :param constant_z: Float of the elevation to find corresponding x coordinates
+    :return: List of x coordinates which correspond to the constant_z
+    """
+
+    # Remove nans to suppress warning messages
+    valid = ~ma.masked_invalid(profile_z).mask
+    profile_z = np.array(profile_z)[valid]
+    profile_x = np.array(profile_x)[valid]
+
+    # Normalize the 'signal' to zero.
+    # Use np.subtract rather than a list comprehension for performance reasons
+    z = np.subtract(profile_z, constant_z)
+
+    # Find all indices right before any crossing.
+    # TODO Sometimes this can give a runtime warning https://stackoverflow.com/a/36489085
+    indicies = np.where(z[:-1] * z[1:] < 0)[0]
+
+    # Use linear interpolation to find intersample crossings.
+    return [profile_x[i] - (profile_x[i] - profile_x[i + 1]) / (z[i] - z[i + 1]) * (z[i]) for i in indicies]
+
+
+
+@click.command()
+@click.option("--sites-csv", required=True, help=".csv file to convert")
+@click.option("--profile-features-csv", required=True, help=".csv file to convert")
+@click.option("--output-geojson", required=True, help="where to store .geojson file")
+def profile_features_to_geojson(sites_csv, profile_features_csv, output_geojson):
+    """
+    Converts profile_features containing dune toes and crest locations to a geojson we can load into QGIS
+    :param sites_csv:
+    :param profiles_csv:
+    :param profile_features_csv:
+    :param output_geojson:
+    :return:
+    """
+    logger.info("Creating profile features geojson")
+
+    # Read files from interim folder
+    df_sites = pd.read_csv(sites_csv, index_col=[0])
+    df_profile_features = pd.read_csv(profile_features_csv, index_col=[0])
+
+    # Create geojson file
+    schema = {
+        'geometry': 'Point',
+        'properties': OrderedDict([
+            ('beach', 'str'),
+            ('site_id', 'str'),
+            ('point_type', 'str'),  # prestorm_dune_toe, prestorm_dune_crest, poststorm_dune_toe, poststorm_dune_crest
+            ('profile_type', 'str'),
+            ('elevation', 'float'),
+        ])
+    }
+
+    with fiona.open(output_geojson, 'w', driver='GeoJSON', crs=from_epsg(4326), schema=schema) as output:
+        for index, row in df_profile_features.iterrows():
+            beach = index[:-4]
+            site_id = index
+            profile_type = row['profile_type']
+
+            for point_type in ['crest', 'toe']:
+                # point_type='crest'
+                elevation = row['dune_{}_z'.format(point_type)]
+                x = row['dune_{}_x'.format(point_type)]
+
+                if np.isnan(x):
+                    continue
+
+                # Geojsons need to use 'null' instead of 'nan'
+                if np.isnan(elevation):
+                    elevation = None
+
+                # Convert x position to lat/lon
+                site = df_sites.loc[site_id]
+                center_profile_x = site["profile_x_lat_lon"]
+                orientation = site["orientation"]
+                center_lat_lon = Point(site['lon'], site['lat'])  # Get lat/lon of center of profile
+                center_xy = convert_coord_systems(center_lat_lon)
+                center_x, center_y = center_xy.xy
+
+                # Calculate xy position of point and convert to lat/lon
+                point_x = center_x + (center_profile_x - x) * np.cos(np.deg2rad(orientation))
+                point_y = center_y + (center_profile_x - x) * np.sin(np.deg2rad(orientation))
+                point_xy = Point(point_x, point_y)
+                point_lat_lon = convert_coord_systems(point_xy,
+                                                     in_coord_system="EPSG:28356",
+                                                     out_coord_system="EPSG:4326")
+
+                prop = OrderedDict([
+                    ('beach', beach),
+                    ('site_id',site_id),
+                    ('point_type', point_type),
+                    ('profile_type', profile_type),
+                    ('elevation', elevation),
+                ])
+                output.write({"geometry": mapping(point_lat_lon), "properties": prop})
+
+
+
 @click.command()
 @click.option("--input-csv", required=True, help=".csv file to convert")
-@click.option("--output-shp", required=True, help="where to store .shp file")
+@click.option("--output-geojson", required=True, help="where to store .geojson file")
-def sites_csv_to_shp(input_csv, output_shp):
+def sites_csv_to_geojson(input_csv, output_geojson):
     """
-    Converts our dataframe of sites to .shp to load in QGis
+    Converts our dataframe of sites to .geojson to load in QGis. Sites are loaded as linestrings of the profile
+    cross-sections
     :param input_csv:
-    :param output_shp:
+    :param output_geojson:
     :return:
     """
-    logger.info("Converting %s to %s", input_csv, output_shp)
+    logger.info("Converting %s to %s", input_csv, output_geojson)
     df_sites = pd.read_csv(input_csv, index_col=[0])
     logger.info(os.environ.get("GDAL_DATA", None))
 
-    schema = {"geometry": "LineString", "properties": {"beach": "str", "site_id": "str"}}
+    schema = {
-    with fiona.open(output_shp, "w", crs=from_epsg(4326), driver="ESRI Shapefile", schema=schema) as output:
+        'geometry': 'LineString',
+        'properties': OrderedDict([
+            ('beach','str'),
+            ('site_id', 'str'),
+        ])
+    }
+
+    with fiona.open(output_geojson, 'w', driver='GeoJSON', crs=from_epsg(4326), schema=schema) as output:
         for index, row in df_sites.iterrows():
             # Work out where center of profile is
             orientation = row["orientation"]
@@ -46,7 +229,8 @@ def sites_csv_to_shp(input_csv, output_shp):
             land_x = center_x + center_profile_x * np.cos(np.deg2rad(orientation))
             land_y = center_y + center_profile_x * np.sin(np.deg2rad(orientation))
             land_xy = Point(land_x, land_y)
-            land_lat_lon = convert_coord_systems(land_xy, in_coord_system="EPSG:28356", out_coord_system="EPSG:4326")
+            land_lat_lon = convert_coord_systems(land_xy, in_coord_system="EPSG:28356",
+                                                 out_coord_system="EPSG:4326")
 
             # Work out where seaward profile limit is
             sea_x = center_x - center_profile_x * np.cos(np.deg2rad(orientation))
@@ -55,7 +239,95 @@ def sites_csv_to_shp(input_csv, output_shp):
             sea_lat_lon = convert_coord_systems(sea_xy, in_coord_system="EPSG:28356", out_coord_system="EPSG:4326")
 
             line_string = LineString([land_lat_lon, center_lat_lon, sea_lat_lon])
-            prop = {"beach": row["beach"], "site_id": index}
+            prop = OrderedDict([("beach", row["beach"]),
+                                ("site_id", index)])
             output.write({"geometry": mapping(line_string), "properties": prop})
 
     logger.info("Done!")
+
+
+@click.command()
+@click.option("--sites-csv", required=True, help="sites.csv file to convert")
+@click.option("--observed-impacts-csv", required=True, help="impacts-observed.csv file to convert")
+@click.option("--forecast-impacts-csv", required=True, help="impacts-forecast.csv file to convert")
+@click.option("--output-geojson", required=True, help="where to store .geojson file")
+def impacts_to_geojson(sites_csv, observed_impacts_csv, forecast_impacts_csv, output_geojson):
+    """
+    Converts impacts observed and forecasted to a geojson for visualization in QGIS
+    :param sites_csv:
+    :param observed_impacts_csv:
+    :param forecast_impacts_csv:
+    :param output_geojson:
+    :return:
+    """
+
+    # Get information from .csv and read into pandas dataframe
+    df_sites = pd.read_csv(sites_csv, index_col=[0])
+    df_observed = pd.read_csv(observed_impacts_csv, index_col=[0])
+    df_forecast = pd.read_csv(forecast_impacts_csv, index_col=[0]).rename({'storm_regime': 'forecast_storm_regime'})
+
+    # Rename columns, so we can distinguish between forecast and observed
+    df_observed = df_observed.rename(columns={'storm_regime': 'observed_storm_regime'})
+    df_forecast = df_forecast.rename(columns={'storm_regime': 'forecast_storm_regime'})
+
+    # Concat into one big dataframe
+    df = pd.concat([df_sites, df_observed, df_forecast], sort=True,axis=1)
+
+    # Make new column for accuracy of forecast. Use underpredict/correct/overpredict classes
+    df.loc[df.observed_storm_regime == df.forecast_storm_regime, 'forecast_accuray'] = 'correct'
+
+    # Observed/Forecasted/Class for each combination
+    classes = [('swash', 'collision', 'overpredict'),
+               ('swash', 'swash', 'correct'),
+               ('swash', 'overwash', 'overpredict'),
+               ('collision', 'swash', 'underpredict'),
+               ('collision', 'collision', 'correct'),
+               ('collision', 'overwash', 'overpredict'),
+               ('overwash', 'swash', 'underpredict'),
+               ('overwash', 'collision', 'underpredict'),
+               ('overwash', 'overwash', 'correct')]
+    for c in classes:
+        df.loc[(df.observed_storm_regime == c[0])&(df.forecast_storm_regime == c[1]), 'forecast_accuracy'] = c[2]
+
+    schema = {
+        'geometry': 'Point',
+        'properties': OrderedDict([
+            ('beach','str'),
+            ('site_id', 'str'),
+            ('forecast_storm_regime', 'str'),
+            ('observed_storm_regime', 'str',),
+            ('forecast_accuracy', 'str')
+        ])
+    }
+
+    # TODO Impact marker location should be at the seaward end of the profile
+    with fiona.open(output_geojson, 'w', driver='GeoJSON', crs=from_epsg(4326), schema=schema) as output:
+        for index, row in df.iterrows():
+
+            # Locate the marker at the seaward end of the profile to avoid cluttering the coastline.
+            # Work out where seaward profile limit is
+            orientation = row["orientation"]
+            center_profile_x = row["profile_x_lat_lon"]
+            center_lon = row["lon"]
+            center_lat = row["lat"]
+            center_lat_lon = Point(center_lon, center_lat)
+            center_xy = convert_coord_systems(center_lat_lon)
+            center_x, center_y = center_xy.xy
+
+            sea_x = center_x - center_profile_x * np.cos(np.deg2rad(orientation))
+            sea_y = center_y - center_profile_x * np.sin(np.deg2rad(orientation))
+            sea_xy = Point(sea_x, sea_y)
+            sea_lat_lon = convert_coord_systems(sea_xy, in_coord_system="EPSG:28356", out_coord_system="EPSG:4326")
+
+            prop = OrderedDict([
+                ('beach',row["beach"]),
+                ('site_id', index),
+                ('forecast_storm_regime', row["forecast_storm_regime"]),
+                ('observed_storm_regime', row["observed_storm_regime"],),
+                ('forecast_accuracy', row["forecast_accuracy"])
+            ])
+
+            output.write({"geometry": mapping(sea_lat_lon), "properties": prop})
+
+
+    logger.info("Done!")