forked from kilianv/CoastSat_WRL
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.7 KiB
Python
89 lines
2.7 KiB
Python
7 years ago
|
# -*- coding: utf-8 -*-
|
||
|
|
|
||
|
|
#==========================================================#
|
||
|
|
# Process shorelines (clipping and smoothing)
|
||
|
|
#==========================================================#
|
||
|
|
|
||
|
|
# Initial settings
|
||
|
|
|
||
|
|
import os
|
||
|
|
import numpy as np
|
||
|
|
import matplotlib.pyplot as plt
|
||
|
|
import pdb
|
||
|
|
import ee
|
||
|
|
|
||
|
|
import matplotlib.dates as mdates
|
||
|
|
import matplotlib.cm as cm
|
||
|
|
import matplotlib.colors as mcolor
|
||
|
|
from datetime import datetime, timedelta
|
||
|
|
import pickle
|
||
|
|
import pytz
|
||
|
|
import scipy.io as sio
|
||
|
|
import scipy.interpolate as interpolate
|
||
|
|
import statsmodels.api as sm
|
||
|
|
import skimage.measure as measure
|
||
|
|
import simplekml
|
||
|
|
|
||
|
|
# my functions
|
||
|
|
import functions.utils as utils
|
||
|
|
import functions.sds as sds
|
||
|
|
|
||
|
|
# some settings
|
||
|
|
np.seterr(all='ignore') # raise/ignore divisions by 0 and nans
|
||
|
|
plt.rcParams['axes.grid'] = True
|
||
|
|
plt.rcParams['figure.max_open_warning'] = 100
|
||
|
|
|
||
|
|
au_tz = pytz.timezone('Australia/Sydney')
|
||
|
|
au_epsg = 28356
|
||
|
|
# load the satellite-derived shorelines
|
||
|
|
satname = 'L8'
|
||
|
|
sitename = 'OLDBAR'
|
||
|
|
filepath = os.path.join(os.getcwd(), 'data', satname, sitename)
|
||
|
|
with open(os.path.join(filepath, sitename + '_output' + '.pkl'), 'rb') as f:
|
||
|
|
output = pickle.load(f)
|
||
|
|
|
||
|
|
sl = output['shorelines']
|
||
|
|
dates_sl = output['t']
|
||
|
|
# convert to AEST
|
||
|
|
dates_sl = [_.astimezone(au_tz) for _ in dates_sl]
|
||
|
|
|
||
|
|
# load the reference shoreline points
|
||
|
|
with open(os.path.join(os.getcwd(), 'data', satname, sitename, sitename + '_refpoints.pkl'), 'rb') as f:
|
||
|
|
refpoints = pickle.load(f)
|
||
|
|
|
||
|
|
dist_thresh = 200
|
||
|
|
frac_smooth = 1./15
|
||
|
|
plt.figure()
|
||
|
|
plt.axis('equal')
|
||
|
|
cmap = cm.get_cmap('brg')
|
||
|
|
colours = cmap(np.linspace(0, 1, num=len(sl)))
|
||
|
|
kml = simplekml.Kml()
|
||
|
|
for i in range(len(sl)):
|
||
|
|
# select points of SDS that are close to the manually digitized points
|
||
|
|
idx_ref = [np.min(np.linalg.norm(sl[i][k,:] - refpoints, axis=1)) < dist_thresh for k in range(sl[i].shape[0])]
|
||
|
|
# smooth (LOWESS) satellite shoreline
|
||
|
|
sl_smooth = sm.nonparametric.lowess(sl[i][idx_ref,0],sl[i][idx_ref,1], frac=frac_smooth, it = 10)
|
||
|
|
sl_smooth = sl_smooth[:,[1,0]]
|
||
|
|
# sl_smooth = sl[i][idx_ref,:]
|
||
|
|
|
||
|
|
# plt.plot(sl[i][idx_ref,0],sl[i][idx_ref,1], 'k-')
|
||
|
|
plt.plot(sl_smooth[:,0], sl_smooth[:,1], color=colours[i,:], label=dates_sl[i].strftime('%d-%b-%Y'))
|
||
|
|
|
||
|
|
# convert to wgs84 (epsg = 4326)
|
||
|
|
sl_wgs84 = sds.convert_epsg(sl_smooth, 28356, 4326)
|
||
|
|
# save in kml file
|
||
|
|
ln = kml.newlinestring(name=dates_sl[i].strftime('%d-%b-%Y'))
|
||
|
|
ln.coords = sl_wgs84
|
||
|
|
ln.style.labelstyle.color = mcolor.rgb2hex(colours[i,:3])
|
||
|
|
ln.style.linestyle.color = mcolor.rgb2hex(colours[i,:3])
|
||
|
|
|
||
|
|
plt.legend(ncol=3)
|
||
|
|
plt.xlabel('Eastings [m]')
|
||
|
|
plt.ylabel('Northings [m]')
|
||
|
|
plt.title('Oldbar inlet (South)')
|
||
|
|
plt.draw()
|
||
|
|
|
||
|
|
kml.save(satname + sitename + '_shorelines.kml')
|
||
|
|
|
||
|
|
|