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.
96 lines
2.8 KiB
Python
96 lines
2.8 KiB
Python
# Converts shoreline chainages to MGA coordinates
|
|
|
|
import os
|
|
import re
|
|
import sys
|
|
import ast
|
|
import json
|
|
import yaml
|
|
import argparse
|
|
import numpy as np
|
|
from scipy import interpolate
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from shapely.geometry import LineString
|
|
import geopandas as gpd
|
|
from tqdm import tqdm
|
|
|
|
sys.path.insert(0, '../lidar/')
|
|
from nielsen import Gridder # noqa
|
|
|
|
MGA55 = 28355 # GDA94, MGA Zone 55
|
|
BEACH = 'Roches'
|
|
|
|
# Load profile data
|
|
xlsx_path = '../lidar/Profiles 1 to 12 2019 DEM.xlsx'
|
|
workbook = pd.ExcelFile(xlsx_path)
|
|
|
|
with open('../lidar/settings.json', 'r') as f:
|
|
sheets = json.loads(f.read())
|
|
|
|
# Create grids to convert chainages to eastings and northings
|
|
grids = {}
|
|
for s in sheets:
|
|
names = ['Chainage', 'Easting', 'Northing', 'Elevation']
|
|
p = workbook.parse(s, names=names)
|
|
k = (sheets[s]['block'], sheets[s]['profile'])
|
|
grids[k] = g = Gridder(chainage=p['Chainage'],
|
|
elevation=p['Elevation'],
|
|
easting=p['Easting'],
|
|
northing=p['Northing'])
|
|
|
|
input_dir = 'output_csv'
|
|
output_dir = 'output_shp'
|
|
|
|
master = gpd.GeoDataFrame(columns=['name', 'year', 'ep', 'type'])
|
|
|
|
# Load probabilistic recession chainages
|
|
file_list = [f for f in os.listdir(input_dir) if f.startswith(BEACH)]
|
|
|
|
df = pd.DataFrame()
|
|
for f in file_list:
|
|
info = re.search(r'(\d{4}) (Z\w+)', f).groups()
|
|
data = pd.read_csv(os.path.join(input_dir, f))
|
|
|
|
ep_cols = [c for c in data.columns if c.startswith('ep_')]
|
|
for i, d in data.iterrows():
|
|
row = {}
|
|
# row['beach'] = BEACH
|
|
row['block'] = d['block']
|
|
row['profile'] = d['profile']
|
|
row['year'] = int(info[0])
|
|
row['type'] = info[1]
|
|
for e in ep_cols:
|
|
row['ep'] = e[3:]
|
|
ch = d[e]
|
|
g = grids[d['block'], d['profile']] # Get correct gridder
|
|
east, north = g.from_chainage(ch)
|
|
row['easting'] = east.round(3)
|
|
row['northing'] = north.round(3)
|
|
df = pd.concat([df, pd.DataFrame([row])])
|
|
|
|
# Convert floats to int
|
|
for col in ['block', 'profile', 'year']:
|
|
df[col] = df[col].astype(int)
|
|
|
|
df = df.set_index(['year', 'type', 'ep', 'block', 'profile']).sort_index()
|
|
years = df.index.unique(level='year')
|
|
types = df.index.unique(level='type')
|
|
eps = df.index.unique(level='ep')
|
|
|
|
lines = []
|
|
for y in years:
|
|
for t in types:
|
|
for e in eps:
|
|
line = {}
|
|
line['year'] = y
|
|
line['type'] = t
|
|
line['ep'] = str(e)
|
|
line['geometry'] = LineString(df.loc[y, t, e].to_numpy())
|
|
lines.append(line)
|
|
|
|
gdf = gpd.GeoDataFrame(lines).set_geometry('geometry')
|
|
gdf = gdf.set_crs(MGA55)
|
|
gdf.to_file(os.path.join(output_dir, 'hazard-lines.shp'),
|
|
driver='ESRI Shapefile')
|