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.
89 lines
2.9 KiB
Python
89 lines
2.9 KiB
Python
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 Point, 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:
|
|
points = []
|
|
for profile in profiles:
|
|
try:
|
|
# Get data for current block/profile (if it exists)
|
|
survey = df.loc[(beach, block, profile, date), :]
|
|
except KeyError:
|
|
continue
|
|
|
|
# Reverse profile direction so elevation is increasing
|
|
survey = survey[::-1]
|
|
survey = survey.set_index('Chainage')
|
|
|
|
# Ignore profile information past high point
|
|
high_idx = survey['Elevation'].idxmax()
|
|
survey = survey[:high_idx]
|
|
|
|
# Extract coordinates for current profile
|
|
elevation = survey['Elevation'].values
|
|
eastings = survey['Easting'].values
|
|
northings = survey['Northing'].values
|
|
|
|
# Check if entire profile is above reference elevation
|
|
if np.all(elevation > contour_z):
|
|
continue
|
|
|
|
# Find first points on either side of reference elevation
|
|
crossing_idx = np.where(elevation < contour_z)[0][-1]
|
|
idx = slice(crossing_idx, crossing_idx + 2)
|
|
|
|
# Calculate shoreline intersection coordinates
|
|
x = np.interp(contour_z, elevation[idx], eastings[idx])
|
|
y = np.interp(contour_z, elevation[idx], northings[idx])
|
|
if x and y:
|
|
points.append(Point(x, y))
|
|
|
|
# Join points from same date into a line
|
|
if len(points) > 1:
|
|
line = LineString(points)
|
|
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')
|