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