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.

77 lines
2.4 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')
# Reverse survey chainage so elevation is increasing
survey = survey[::-1]
# Find largest chainage in profile above contour elevation
idx = survey.where(
survey['Elevation'] > contour_z).first_valid_index()
elevation = survey.loc[:idx, 'Elevation']
eastings = survey.loc[:idx, 'Easting']
northings = survey.loc[:idx, 'Northing']
x = np.interp(contour_z, elevation, eastings)
y = np.interp(contour_z, elevation, northings)
5 years ago
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')