You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

144 lines
5.6 KiB
Python

"""
This file can probably be removed since we are not reading from shp files any more
"""
import click
6 years ago
import fiona
6 years ago
import numpy as np
import pandas as pd
6 years ago
from shapely.geometry import LineString, Point
from shapely.geometry import shape
from logs import setup_logging
from utils import convert_coord_systems
logger = setup_logging()
6 years ago
6 years ago
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):
6 years ago
"""
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
6 years ago
towards the
land.
:param line_string:
:return:
"""
6 years ago
start_point = Point(lon, lat)
6 years ago
start_point = convert_coord_systems(start_point)
6 years ago
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.
6 years ago
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)
6 years ago
6 years ago
# If no intersections are found, return nothing.
6 years ago
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}}
6 years ago
# Import our dune crest and toes
for feat in features.keys():
shapes, _, properties = shapes_from_shp(features[feat]["file"])
6 years ago
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
6 years ago
@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!")