nsw-2016-storm-impact/src/data/parse_shp.py

from functools import partial

import click
import fiona
import numpy as np
import pandas as pd
import pyproj
from shapely.geometry import LineString, Point
from shapely.geometry import shape
from shapely.ops import transform

from utils import setup_logging

logger = setup_logging()


def convert_coord_systems(g1, in_coord_system="EPSG:4326", out_coord_system="EPSG:28356"):
    """
    Converts coordinates from one coordinates system to another. Needed because shapefiles are usually defined in
    lat/lon but should be converted to GDA to calculated distances.
    https://gis.stackexchange.com/a/127432
    :param in_coord_system: Default is lat/lon WGS84
    :param out_coord_system: Default is GDA56 for NSW coastline
    :return:
    """
    project = partial(
        pyproj.transform,
        pyproj.Proj(init=in_coord_system),  # source coordinate system
        pyproj.Proj(init=out_coord_system),
    )  # destination coordinate system

    g2 = transform(project, g1)  # apply projection
    return g2


def shapes_from_shp(shp_file):
    """
    Parses a shape file and returns a list of shapely shapes, ids and properties
    :param shp_file:
    :return:
    """
    shapes = []
    ids = []
    properties = []
    for feat in fiona.open(shp_file, "r"):
        shapes.append(shape(feat["geometry"]))
        ids.append(feat["id"])
        properties.append(feat["properties"])
    return shapes, ids, properties


def distance_to_intersection(lat, lon, landward_orientation, beach, line_strings, line_properties):
    """
    Returns the distance at whjch a line drawn from a lat/lon at an orientation intersects a line stinrg
    :param lat:
    :param lon:
    :param landward_orientation: Angle, anticlockwise positive from east in degrees, towards the land
    towards the
    land.
    :param line_string:
    :return:
    """
    start_point = Point(lon, lat)
    start_point = convert_coord_systems(start_point)

    distance = 1000  # m look up to 1000m for an intersection
    landward_point = Point(
        start_point.coords.xy[0] + distance * np.cos(np.deg2rad(landward_orientation)),
        start_point.coords.xy[1] + distance * np.sin(np.deg2rad(landward_orientation)),
    )
    landward_line = LineString([start_point, landward_point])
    seaward_point = Point(
        start_point.coords.xy[0] - distance * np.cos(np.deg2rad(landward_orientation)),
        start_point.coords.xy[1] - distance * np.sin(np.deg2rad(landward_orientation)),
    )
    seaward_line = LineString([start_point, seaward_point])

    # Look at relevant line_strings which have the same beach property in order to reduce computation time
    line_strings = [s for s, p in zip(line_strings, line_properties) if p["beach"] == beach]

    # Check whether profile_line intersects with any lines in line_string. If intersection point is landwards,
    # consider this negative, otherwise seawards is positive.
    for line_string in line_strings:
        land_intersect_points = landward_line.intersection(line_string)
        if not land_intersect_points.is_empty:
            return -land_intersect_points.distance(start_point)

        sea_intersect_points = seaward_line.intersection(line_string)
        if not sea_intersect_points.is_empty:
            return sea_intersect_points.distance(start_point)

    # If no intersections are found, return nothing.
    return None


def beach_profile_elevation(x_coord, df_profiles, profile_type, site_id):
    """
    Returns the beach profile elevation at a particular x-coordinate
    :param x_coord:
    :param df_profiles:
    :param profile_type: "prestorm" or "poststorm"
    :param site_id:
    :return:
    """

    if np.isnan(x_coord):
        return None

    # Get profile
    df_profile = df_profiles.query('profile_type == "{}" and site_id =="{}"'.format(profile_type, site_id))
    return np.interp(x_coord, df_profile.index.get_level_values("x"), df_profile["z"])


def parse_profile_features(df_sites, df_profiles, dune_crest_shp, dune_toe_shp):
    """
    Reads dune crest and toe files and creates a pandas dataframe with crest/toe locations at each site
    :return:
    """

    # Get site information. Base our profile features on each site
    df_profile_features = df_sites

    features = {"dune_crest": {"file": dune_crest_shp}, "dune_toe": {"file": dune_toe_shp}}

    # Import our dune crest and toes
    for feat in features.keys():
        shapes, _, properties = shapes_from_shp(features[feat]["file"])
        shapes = [convert_coord_systems(x) for x in shapes]

        # Figure out the x coordinates of our crest and toes, by looking at where our beach sections intersect our
        # shape files.
        col_name = "{}_x".format(feat)
        df_profile_features[col_name] = df_profile_features["profile_x_lat_lon"] + df_profile_features.apply(
            lambda row: distance_to_intersection(
                row["lat"], row["lon"], row["orientation"], row["beach"], shapes, properties
            ),
            axis=1,
        )
        # Get the elevations of the crest and toe
        col_name = "{}_z".format(feat)
        df_profile_features[col_name] = df_profile_features.apply(
            lambda row: beach_profile_elevation(row["{}_x".format(feat)], df_profiles, "prestorm", row.name), axis=1
        )

    df_profile_features = df_profile_features.drop(columns=["beach", "lat", "lon", "orientation", "profile_x_lat_lon"])
    return df_profile_features


@click.command(short_help="create .csv of dune toe and crest positions")
@click.option("--dune-crest-shp", required=True, help=".csv file to convert")
@click.option("--dune-toe-shp", required=True, help="where to store .shp file")
@click.option("--sites-csv", required=True, help="where to store .shp file")
@click.option("--profiles-csv", required=True, help="where to store .shp file")
@click.option("--output-csv", required=True, help="where to store .shp file")
def create_profile_features(dune_crest_shp, dune_toe_shp, sites_csv, profiles_csv, output_csv):
    logger.info("Creating .csv of dune crests and toes")
    df_sites = pd.read_csv(sites_csv, index_col=[0])
    df_profiles = pd.read_csv(profiles_csv, index_col=[0, 1, 2])
    df_profile_features = parse_profile_features(df_sites, df_profiles, dune_crest_shp, dune_toe_shp)
    df_profile_features.to_csv(output_csv)
    logger.info("Done!")