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.

162 lines
6.2 KiB
Python

6 years ago
from functools import partial
import click
6 years ago
import fiona
6 years ago
import numpy as np
import pandas as pd
import pyproj
6 years ago
from shapely.geometry import LineString, Point
from shapely.geometry import shape
from shapely.ops import transform
from logs import setup_logging
logger = setup_logging()
6 years ago
6 years ago
def convert_coord_systems(g1, in_coord_system="EPSG:4326", out_coord_system="EPSG:28356"):
6 years ago
"""
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,
6 years ago
pyproj.Proj(init=in_coord_system), # source coordinate system
pyproj.Proj(init=out_coord_system),
) # destination coordinate system
6 years ago
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):
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!")