nsw-beach-profile-database/extract_contours.py

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 for interpolation
                    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)

                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')