|
|
@ -1,14 +1,15 @@
|
|
|
|
import pandas as pd
|
|
|
|
|
|
|
|
import os
|
|
|
|
import os
|
|
|
|
|
|
|
|
from functools import partial
|
|
|
|
|
|
|
|
|
|
|
|
import fiona
|
|
|
|
import fiona
|
|
|
|
|
|
|
|
import numpy as np
|
|
|
|
|
|
|
|
import pandas as pd
|
|
|
|
|
|
|
|
import pyproj
|
|
|
|
|
|
|
|
from scipy import stats
|
|
|
|
from shapely.geometry import LineString, Point
|
|
|
|
from shapely.geometry import LineString, Point
|
|
|
|
from shapely.geometry import shape
|
|
|
|
from shapely.geometry import shape
|
|
|
|
from shapely.ops import transform
|
|
|
|
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):
|
|
|
|
def shapes_from_shp(shp_file):
|
|
|
|
"""
|
|
|
|
"""
|
|
|
@ -33,8 +34,8 @@ def convert_coord_systems(g1, in_coord_system='EPSG:4326', out_coord_system='EPS
|
|
|
|
"""
|
|
|
|
"""
|
|
|
|
project = partial(
|
|
|
|
project = partial(
|
|
|
|
pyproj.transform,
|
|
|
|
pyproj.transform,
|
|
|
|
pyproj.Proj(init=in_coord_system), # source coordinate system
|
|
|
|
pyproj.Proj(init=in_coord_system), # source coordinate system
|
|
|
|
pyproj.Proj(init=out_coord_system)) # destination coordinate system
|
|
|
|
pyproj.Proj(init=out_coord_system)) # destination coordinate system
|
|
|
|
|
|
|
|
|
|
|
|
g2 = transform(project, g1) # apply projection
|
|
|
|
g2 = transform(project, g1) # apply projection
|
|
|
|
return g2
|
|
|
|
return g2
|
|
|
@ -51,7 +52,7 @@ def get_slope(x, z, top_elevation, btm_elevation, method='end_points'):
|
|
|
|
:return:
|
|
|
|
:return:
|
|
|
|
"""
|
|
|
|
"""
|
|
|
|
|
|
|
|
|
|
|
|
profile_line = LineString([[x,z] for x,z in zip(x,z)])
|
|
|
|
profile_line = LineString([[x, z] for x, z in zip(x, z)])
|
|
|
|
|
|
|
|
|
|
|
|
end_points = {
|
|
|
|
end_points = {
|
|
|
|
'top': {
|
|
|
|
'top': {
|
|
|
@ -88,14 +89,14 @@ def get_slope(x, z, top_elevation, btm_elevation, method='end_points'):
|
|
|
|
return -(z_top - z_btm) / (x_top - x_btm)
|
|
|
|
return -(z_top - z_btm) / (x_top - x_btm)
|
|
|
|
|
|
|
|
|
|
|
|
elif method == 'least_squares':
|
|
|
|
elif method == 'least_squares':
|
|
|
|
profile_mask = [True if end_points['top']['x'] < pts < end_points['btm']['x'] else False for pts in x ]
|
|
|
|
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_x = np.array(x)[profile_mask].tolist()
|
|
|
|
profile_z = np.array(z)[profile_mask].tolist()
|
|
|
|
profile_z = np.array(z)[profile_mask].tolist()
|
|
|
|
slope, _, _, _, _ = stats.linregress(profile_x, profile_z)
|
|
|
|
slope, _, _, _, _ = stats.linregress(profile_x, profile_z)
|
|
|
|
return -slope
|
|
|
|
return -slope
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def distance_to_intersection(lat,lon,orientation,line_strings):
|
|
|
|
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
|
|
|
|
Returns the distance at whjch a line drawn from a lat/lon at an orientation intersects a line stinrg
|
|
|
|
:param lat:
|
|
|
|
:param lat:
|
|
|
@ -106,12 +107,12 @@ def distance_to_intersection(lat,lon,orientation,line_strings):
|
|
|
|
:param line_string:
|
|
|
|
:param line_string:
|
|
|
|
:return:
|
|
|
|
:return:
|
|
|
|
"""
|
|
|
|
"""
|
|
|
|
start_point = Point(lon,lat)
|
|
|
|
start_point = Point(lon, lat)
|
|
|
|
start_point = convert_coord_systems(start_point)
|
|
|
|
start_point = convert_coord_systems(start_point)
|
|
|
|
|
|
|
|
|
|
|
|
distance = 1000 # m look up to 1000m for an intersection
|
|
|
|
distance = 1000 # m look up to 1000m for an intersection
|
|
|
|
new_point = Point(start_point.coords.xy[0]+distance*np.cos(np.deg2rad(orientation)),
|
|
|
|
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)))
|
|
|
|
start_point.coords.xy[1] + distance * np.sin(np.deg2rad(orientation)))
|
|
|
|
profile_line = LineString([start_point, new_point])
|
|
|
|
profile_line = LineString([start_point, new_point])
|
|
|
|
|
|
|
|
|
|
|
|
# Check whether profile_line intersects with any lines in line_string
|
|
|
|
# Check whether profile_line intersects with any lines in line_string
|
|
|
@ -123,9 +124,10 @@ def distance_to_intersection(lat,lon,orientation,line_strings):
|
|
|
|
# If no intersections are found, return nothing.
|
|
|
|
# If no intersections are found, return nothing.
|
|
|
|
return None
|
|
|
|
return None
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def get_sites_dune_crest_toe():
|
|
|
|
def get_sites_dune_crest_toe():
|
|
|
|
data_folder = './data/interim'
|
|
|
|
data_folder = './data/interim'
|
|
|
|
df_sites = pd.read_csv(os.path.join(data_folder, 'sites.csv'),index_col=[0])
|
|
|
|
df_sites = pd.read_csv(os.path.join(data_folder, 'sites.csv'), index_col=[0])
|
|
|
|
|
|
|
|
|
|
|
|
# Import
|
|
|
|
# Import
|
|
|
|
for f in ['./data/raw/profile_features/dune_crests.shp']:
|
|
|
|
for f in ['./data/raw/profile_features/dune_crests.shp']:
|
|
|
|