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.
136 lines
4.7 KiB
Python
136 lines
4.7 KiB
Python
import pandas as pd
|
|
import os
|
|
import fiona
|
|
from shapely.geometry import LineString, Point
|
|
from shapely.geometry import shape
|
|
from shapely.ops import transform
|
|
import pyproj
|
|
from functools import partial
|
|
import numpy as np
|
|
from scipy.interpolate import interp1d
|
|
from scipy import stats
|
|
|
|
def shapes_from_shp(shp_file):
|
|
"""
|
|
Parses a shape file and returns a list of shapely shapes
|
|
:param shp_file:
|
|
:return:
|
|
"""
|
|
shapes = []
|
|
for feat in fiona.open(shp_file):
|
|
shapes.append(shape(feat['geometry']))
|
|
return shapes
|
|
|
|
|
|
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 get_slope(x, z, top_elevation, btm_elevation, method='end_points'):
|
|
"""
|
|
Returns a slope from a bed profile
|
|
:param x: List of x bed profile coordinates
|
|
:param z: List of z bed profile coordinates
|
|
:param top_elevation: Top elevation of where to take the slope
|
|
:param btm_elevation: Bottom elevation of where to take the slope
|
|
:param method: Method used to calculate slope (end_points or least_squares)
|
|
:return:
|
|
"""
|
|
|
|
profile_line = LineString([[x,z] for x,z in zip(x,z)])
|
|
|
|
end_points = {
|
|
'top': {
|
|
'elevation': top_elevation,
|
|
},
|
|
'btm': {
|
|
'elevation': btm_elevation,
|
|
}}
|
|
|
|
# If there are multiple intersections, take most seaward of landward point?
|
|
multiple_points = 'landward'
|
|
|
|
for end_type in end_points.keys():
|
|
elevation = end_points[end_type]['elevation']
|
|
elevation_line = LineString([[min(x), elevation], [max(x), elevation]])
|
|
intersection_points = profile_line.intersection(elevation_line)
|
|
|
|
if intersection_points.is_empty:
|
|
return None
|
|
|
|
if multiple_points == 'landward':
|
|
intersection_point = list(intersection_points)[0]
|
|
elif multiple_points == 'seaward':
|
|
intersection_point = list(intersection_points)[-1]
|
|
|
|
end_points[end_type]['x'] = intersection_point.xy[0][0]
|
|
end_points[end_type]['z'] = intersection_point.xy[1][0]
|
|
|
|
if method == 'end_points':
|
|
x_top = end_points['top']['x']
|
|
x_btm = end_points['btm']['x']
|
|
z_top = end_points['top']['z']
|
|
z_btm = end_points['btm']['z']
|
|
return -(z_top - z_btm) / (x_top - x_btm)
|
|
|
|
elif method == 'least_squares':
|
|
profile_mask = [True if end_points['top']['x'] < pts < end_points['btm']['x'] else False for pts in x ]
|
|
profile_x = np.array(x)[profile_mask].tolist()
|
|
profile_z = np.array(z)[profile_mask].tolist()
|
|
slope, _, _, _, _ = stats.linregress(profile_x, profile_z)
|
|
return -slope
|
|
|
|
|
|
def distance_to_intersection(lat,lon,orientation,line_strings):
|
|
"""
|
|
Returns the distance at whjch a line drawn from a lat/lon at an orientation intersects a line stinrg
|
|
:param lat:
|
|
:param lon:
|
|
:param orientation: Angle, clockwise positive from true north in degrees, of the tangent to the shoreline facing
|
|
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
|
|
new_point = Point(start_point.coords.xy[0]+distance*np.cos(np.deg2rad(orientation)),
|
|
start_point.coords.xy[1]+distance*np.sin(np.deg2rad(orientation)))
|
|
profile_line = LineString([start_point, new_point])
|
|
|
|
# Check whether profile_line intersects with any lines in line_string
|
|
for line_string in line_strings:
|
|
intersection_points = profile_line.intersection(line_string)
|
|
if not intersection_points.is_empty:
|
|
return intersection_points.distance(start_point)
|
|
|
|
# If no intersections are found, return nothing.
|
|
return None
|
|
|
|
def get_sites_dune_crest_toe():
|
|
data_folder = './data/interim'
|
|
df_sites = pd.read_csv(os.path.join(data_folder, 'sites.csv'),index_col=[0])
|
|
|
|
# Import
|
|
for f in ['./data/raw/profile_features/dune_crests.shp']:
|
|
shapes = shapes_from_shp(f)
|
|
shapes = [convert_coord_systems(x) for x in shapes]
|
|
|
|
# Iterate through each site
|