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.
74 lines
2.3 KiB
Python
74 lines
2.3 KiB
Python
5 years ago
|
import os
|
||
|
|
import re
|
||
|
|
import fiona
|
||
|
|
import numpy as np
|
||
|
|
import pandas as pd
|
||
|
|
import geopandas as gpd
|
||
|
|
import matplotlib.pyplot as plt
|
||
|
|
from shapely.geometry import LineString
|
||
|
|
|
||
|
|
# Set contour elevation
|
||
|
|
contour_z = 0.7
|
||
|
|
|
||
|
|
# Set coordinate reference system
|
||
|
|
epsg = 28356
|
||
|
|
|
||
|
|
csv_name = 'photogrammetry_Xsections_Wamberal and Terrigal.csv'
|
||
|
|
|
||
|
|
# Load csv
|
||
|
|
df = pd.read_csv(csv_name, header=5)
|
||
|
|
df = df.set_index(['Beach', 'Block', 'Profile', 'Year/Date'])
|
||
|
|
|
||
|
|
# Get unique survey parameters
|
||
|
|
beaches = df.index.get_level_values('Beach').unique()
|
||
|
|
blocks = df.index.get_level_values('Block').unique()
|
||
|
|
profiles = df.index.get_level_values('Profile').unique()
|
||
|
|
dates = df.index.get_level_values('Year/Date').unique()
|
||
|
|
|
||
|
|
# Sort dataframe
|
||
|
|
df = df.sort_index()
|
||
|
|
|
||
|
|
contours = []
|
||
|
|
for beach in beaches:
|
||
|
|
for block in blocks:
|
||
|
|
for date in dates:
|
||
|
|
x = []
|
||
|
|
y = []
|
||
|
|
for profile in profiles:
|
||
|
|
try:
|
||
|
|
# Get data for current profile
|
||
|
|
survey = df.loc[(beach, block, profile, date), :]
|
||
|
|
survey = survey.set_index('Chainage')
|
||
|
|
|
||
|
|
# Find last point in profile above ref contour elevation
|
||
|
|
last_idx = survey.where(
|
||
|
|
survey['Elevation'] > contour_z).last_valid_index()
|
||
|
|
elevation = survey.loc[last_idx:, 'Elevation']
|
||
|
|
eastings = survey.loc[last_idx:, 'Easting']
|
||
|
|
northings = survey.loc[last_idx:, 'Northing']
|
||
|
|
x.append(np.interp(contour_z, elevation, eastings))
|
||
|
|
y.append(np.interp(contour_z, elevation, northings))
|
||
|
|
|
||
|
|
except KeyError:
|
||
|
|
pass
|
||
|
|
|
||
|
|
if x and y:
|
||
|
|
# Join points from same date into a line
|
||
|
|
line = LineString([(x1, y1) for x1, y1 in zip(x, y)])
|
||
|
|
contours.append({
|
||
|
|
'date': date,
|
||
|
|
'beach': beach,
|
||
|
|
'block': block,
|
||
|
|
'geometry': line
|
||
|
|
})
|
||
|
|
|
||
|
|
# Make geodataframe
|
||
|
|
gdf = gpd.GeoDataFrame(contours)
|
||
|
|
|
||
|
|
# Set coordinate reference system
|
||
|
|
gdf.crs = fiona.crs.from_epsg(epsg)
|
||
|
|
|
||
|
|
beach_name = csv_name.replace('photogrammetry_Xsections_', '')
|
||
|
|
shp_name = beach_name.replace('.csv', f' {contour_z}m contour.shp')
|
||
|
|
gdf.to_file(shp_name, driver='ESRI Shapefile')
|