From e50fadb35c8635ad5c0341c468ad26220a9722c5 Mon Sep 17 00:00:00 2001 From: Chris Leaman Date: Tue, 18 Dec 2018 14:42:57 +1100 Subject: [PATCH] Rename csv_to_shp to csv_to_geojson --- src/data/csv_to_geojson.py | 298 +++++++++++++++++++++++++++++++++ src/data/csv_to_shp.py | 333 ------------------------------------- 2 files changed, 298 insertions(+), 333 deletions(-) create mode 100644 src/data/csv_to_geojson.py delete mode 100644 src/data/csv_to_shp.py diff --git a/src/data/csv_to_geojson.py b/src/data/csv_to_geojson.py new file mode 100644 index 0000000..40516ef --- /dev/null +++ b/src/data/csv_to_geojson.py @@ -0,0 +1,298 @@ +""" +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 utils import crossings, convert_coord_systems +from logs import setup_logging + +logger = setup_logging() + + +def lat_lon_from_profile_x_coord(center_lat_lon, orientation, center_profile_x, x_coord): + """ + Returns the lat/lon of a point on a profile with the given x_coord + :param center_lat_lon: Shapely point of lat/lon of profile center + :param orientation: Orientation of the profile (positive east, counterclockwise) + :param center_profile_x: x value of the center of the profile + :param x_coord: X coordinate of the point to get a lat lon from + :return: + """ + center_xy = convert_coord_systems(center_lat_lon) + center_x, center_y = center_xy.xy + + point_x = center_x + (center_profile_x - x_coord) * np.cos(np.deg2rad(orientation)) + point_y = center_y + (center_profile_x - x_coord) * 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") + return point_lat_lon + + +@click.command() +@click.option("--sites-csv", required=True, help=".csv file to convert") +@click.option("--profile-csv", required=True, help=".csv file to convert") +@click.option("--impacts-csv", required=True, help=".csv file to convert") +@click.option("--output-geojson", required=True, help="where to store .geojson file") +def R_high_to_geojson(sites_csv, profiles_csv, impacts_csv, output_geojson): + """ + Converts impact R_high into a lat/lon geojson that we can plot in QGIS + :param sites_csv: + :param profiles_csv: + :param impacts_csv: + :param output_geojson: + :return: + """ + 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] + point_lat_lon = lat_lon_from_profile_x_coord( + center_lat_lon=Point(site["lon"], site["lat"]), + orientation=site["orientation"], + center_profile_x=site["profile_x_lat_lon"], + x_coord=int_x, + ) + + prop = OrderedDict([("beach", beach), ("site_id", site_id), ("elevation", R_high_z)]) + output.write({"geometry": mapping(point_lat_lon), "properties": prop}) + + +@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] + point_lat_lon = lat_lon_from_profile_x_coord( + center_lat_lon=Point(site["lon"], site["lat"]), + orientation=site["orientation"], + center_profile_x=site["profile_x_lat_lon"], + x_coord=x, + ) + + 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-geojson", required=True, help="where to store .geojson file") +def sites_csv_to_geojson(input_csv, output_geojson): + """ + 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_geojson: + :return: + """ + 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": 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(): + + center_lat_lon = Point(row["lon"], row["lat"]) + # Work out where landward profile limit is + land_lat_lon = lat_lon_from_profile_x_coord( + center_lat_lon=center_lat_lon, + orientation=row["orientation"], + center_profile_x=row["profile_x_lat_lon"], + x_coord=0, + ) + # Work out where seaward profile limit is + sea_lat_lon = lat_lon_from_profile_x_coord( + center_lat_lon=center_lat_lon, + orientation=row["orientation"], + center_profile_x=row["profile_x_lat_lon"], + x_coord=2 * row["profile_x_lat_lon"], + ) + + line_string = LineString([land_lat_lon, center_lat_lon, sea_lat_lon]) + 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"), + ] + ), + } + + 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 + sea_lat_lon = lat_lon_from_profile_x_coord( + center_lat_lon=Point(row["lon"], row["lat"]), + orientation=row["orientation"], + center_profile_x=row["profile_x_lat_lon"], + x_coord=2 * row["profile_x_lat_lon"], + ) + 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!") diff --git a/src/data/csv_to_shp.py b/src/data/csv_to_shp.py deleted file mode 100644 index 96624b5..0000000 --- a/src/data/csv_to_shp.py +++ /dev/null @@ -1,333 +0,0 @@ -""" -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-geojson", required=True, help="where to store .geojson file") -def sites_csv_to_geojson(input_csv, output_geojson): - """ - 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_geojson: - :return: - """ - 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': 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"] - 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 - - # Work out where landward profile limit is - 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") - - # Work out where seaward profile limit is - 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") - - line_string = LineString([land_lat_lon, center_lat_lon, sea_lat_lon]) - 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!")