From 3f2cfa50aadfc773b0e654f257528fbeadc4e27d Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Tue, 18 Dec 2018 13:48:00 +1100 Subject: [PATCH] Update functions to output geojson for QGIS --- src/cli.py | 2 +- src/data/csv_to_shp.py | 294 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 284 insertions(+), 12 deletions(-) diff --git a/src/cli.py b/src/cli.py index 0323345..e227868 100644 --- a/src/cli.py +++ b/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) diff --git a/src/data/csv_to_shp.py b/src/data/csv_to_shp.py index 5f8ecc3..96624b5 100644 --- a/src/data/csv_to_shp.py +++ b/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") -def sites_csv_to_shp(input_csv, output_shp): +@click.option("--output-geojson", required=True, help="where to store .geojson file") +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"}} - with fiona.open(output_shp, "w", crs=from_epsg(4326), driver="ESRI Shapefile", schema=schema) as output: + schema = { + '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!")