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