From 39c5bb05e1d606a1e23f6efc7b071ae72d4a097b Mon Sep 17 00:00:00 2001 From: Kilian Vos Date: Wed, 4 Apr 2018 08:03:02 +1000 Subject: [PATCH] sds interpolation on 3d quadbike surface --- functions/utils.py | 1 - time_coverage.py | 131 ++++++++++++++++++++++++++++++--------------- 2 files changed, 87 insertions(+), 45 deletions(-) diff --git a/functions/utils.py b/functions/utils.py index 35aecca..2bd49dd 100644 --- a/functions/utils.py +++ b/functions/utils.py @@ -9,7 +9,6 @@ Contains all the utilities, convenience functions and small functions that do si import matplotlib.pyplot as plt import numpy as np -import datetime import pdb diff --git a/time_coverage.py b/time_coverage.py index e84ff07..e3620ad 100644 --- a/time_coverage.py +++ b/time_coverage.py @@ -17,15 +17,9 @@ from datetime import datetime, timedelta import pickle import pytz import scipy.io as sio +import scipy.interpolate as interpolate +import statsmodels.api as sm -# image processing modules -import skimage.filters as filters -import skimage.exposure as exposure -import skimage.transform as transform -import skimage.morphology as morphology -import skimage.measure as measure -import sklearn.decomposition as decomposition -from scipy import spatial # my functions import functions.utils as utils import functions.sds as sds @@ -65,6 +59,8 @@ dates_wave = [datetime(wave_data['dates'][i,0], wave_data['dates'][i,1], wave_data['dates'][i,4], wave_data['dates'][i,5], tzinfo=au_tz) for i in idx] #%% make a plot of all the dates +orange = [255/255,140/255,0] +blue = [0,191/255,255/255] f = plt.figure() months = mdates.MonthLocator() month_fmt = mdates.DateFormatter('%b %Y') @@ -80,18 +76,18 @@ for k in range(len(years)): for j in range(len(dates_quad)): if dates_quad[j].year == sel_year: if cbool: - plt.plot([dates_quad[j], dates_quad[j]], [0, hsigmax], color=[255/255,140/255,0], label='survey') + plt.plot([dates_quad[j], dates_quad[j]], [0, hsigmax], color=orange, label='survey') cbool = False else: - plt.plot([dates_quad[j], dates_quad[j]], [0, hsigmax], color=[255/255,140/255,0]) + plt.plot([dates_quad[j], dates_quad[j]], [0, hsigmax], color=orange) cbool = True for j in range(len(dates_l8)): if dates_l8[j].year == sel_year: if cbool: - plt.plot([dates_l8[j], dates_l8[j]], [0, hsigmax], color=[0,191/255,255/255], label='landsat8') + plt.plot([dates_l8[j], dates_l8[j]], [0, hsigmax], color=blue, label='landsat8') cbool = False else: - plt.plot([dates_l8[j], dates_l8[j]], [0, hsigmax], color=[0,191/255,255/255]) + plt.plot([dates_l8[j], dates_l8[j]], [0, hsigmax], color=blue) if k == 3: plt.legend() plt.xlim((datetime(sel_year,1,1), datetime(sel_year,12,31, tzinfo=au_tz))) @@ -101,9 +97,7 @@ for k in range(len(years)): ax.xaxis.set_major_formatter(month_fmt) f.subplots_adjust(hspace=0.2) plt.draw() -#%% - -# calculate days difference +#%% calculate days difference diff_days = [ [(x - _).days for _ in dates_quad] for x in dates_l8] max_diff = 5 idx_closest = [utils.find_indices(_, lambda e: abs(e) <= max_diff) for _ in diff_days] @@ -126,7 +120,11 @@ for i in range(len(idx_closest)): np.mean([ np.abs(_['days diff']) for _ in dates_diff]) -#%% +#%% compare shorelines + +dist_thresh = 200 # maximum distance between an sds point and a narrabeen point +frac_smooth = 1./12 # fraction of the data used for smoothing (the bigger the smoother) +dist_buffer = 50 # buffer of points selected for interpolation # load quadbike .mat files foldername = 'data\quadbike\surveys3D' @@ -136,43 +134,88 @@ filenames = os.listdir(folderpath) # load the satellite shorelines sl = output['shorelines'] +# load narrabeen beach points (manually digitized) +with open(os.path.join(os.getcwd(), 'olddata', 'narra_beach' + '.pkl'), 'rb') as f: + narrabeach = pickle.load(f) + dates_quad = [datetime(int(_[6:10]), int(_[11:13]), int(_[14:16]), tzinfo= au_tz) for _ in filenames] -# for each satellite shoreline, load the corresponding 3D survey +zav = [] for i in range(len(dates_diff)): + # select closest 3D survey idx_closest = np.argmin(np.abs(np.array([(dates_diff[i]['date sat'] - _).days for _ in dates_quad]))) survey3d = sio.loadmat(os.path.join(folderpath, filenames[idx_closest])) - + xs = survey3d['x'].reshape(survey3d['x'].shape[0] * survey3d['x'].shape[1]) + ys = survey3d['y'].reshape(survey3d['y'].shape[0] * survey3d['y'].shape[1]) + zs = survey3d['z'].reshape(survey3d['z'].shape[0] * survey3d['z'].shape[1]) + idx_nan = np.isnan(zs) + xs = xs[~idx_nan] + ys = ys[~idx_nan] + zs = zs[~idx_nan] + # smooth (LOWESS) satellite shoreline + idx_beach = [np.min(np.linalg.norm(sl[i][k,:] - narrabeach, axis=1)) < dist_thresh for k in range(sl[i].shape[0])] + sl_smooth = sm.nonparametric.lowess(sl[i][idx_beach,0],sl[i][idx_beach,1], frac=frac_smooth, it = 6) + sl_smooth = sl_smooth[:,[1,0]] + # make plot plt.figure() plt.axis('equal') - plt.scatter(survey3d['x'], survey3d['y'], s=10, c=survey3d['z'], marker='o', cmap=cm.get_cmap('jet'), + plt.scatter(xs, ys, s=10, c=zs, marker='o', cmap=cm.get_cmap('jet'), label='quad data') - plt.plot(sl[i][:,0], sl[i][:,1], 'ko') + plt.plot(sl[i][idx_beach,0], sl[i][idx_beach,1], 'ko-', markersize=3) + plt.plot(sl_smooth[:,0], sl_smooth[:,1], 'ro-', markersize=3) + plt.xlabel('Eastings [m]') + plt.ylabel('Northings [m]') + plt.title('Local weighted scatterplot smoothing (LOWESS)') plt.draw() -import statsmodels.api as sm -lowess = sm.nonparametric.lowess + zq = np.zeros((sl_smooth.shape[0], 1)) + for j in range(sl_smooth.shape[0]): + xq = sl_smooth[j,0] + yq = sl_smooth[j,1] + dist_q = np.linalg.norm(np.transpose(np.array([[xq - _ for _ in xs],[yq - _ for _ in ys]])), axis=1) + idx_buffer = dist_q <= dist_buffer + + +# plt.figure() +# plt.axis('equal') +# plt.scatter(xs, ys, s=10, c=zs, marker='o', cmap=cm.get_cmap('jet'), +# label='quad data') +# plt.plot(xs[idx_buffer], ys[idx_buffer], 'ko') +# plt.plot(xq,yq,'ro') +# plt.draw() + + tck = interpolate.bisplrep(xs[idx_buffer], ys[idx_buffer], zs[idx_buffer]) + zq[j] = interpolate.bisplev(xq, yq, tck) + + zav.append(np.median(zq)) + plt.figure() + plt.plot(sl_smooth[:,1], zq, 'ko-', markersize=5) + plt.plot([sl_smooth[0,1], sl_smooth[-1,1]], [zav[i], zav[i]], 'r--') + plt.xlabel('Northings [m]') + plt.ylabel('Elevation [mAHD]') + plt.title('Interpolated SDS elevation') + plt.draw() -# For the 1D case: -x = sl[i][:,0] -y = sl[i][:,1] -x0 = x -f_hat = lo.lowess(x, y, x) -fig,ax = plt.subplots(1) -ax.scatter(x,y) -ax.plot(x0,f_hat,'ro') -plt.show() - -# 2D case (and more...) -x = np.random.randn(2, 100) -f = -1 * np.sin(x[0]) + 0.5 * np.cos(x[1]) + 0.2*np.random.randn(100) -x0 = np.mgrid[-1:1:.1, -1:1:.1] -x0 = np.vstack([x0[0].ravel(), x0[1].ravel()]) -f_hat = lo.lowess(x, f, x0, kernel=lo.tri_cube) -from mpl_toolkits.mplot3d import Axes3D -fig = plt.figure() -ax = fig.add_subplot(111, projection='3d') -ax.scatter(x[0], x[1], f) -ax.scatter(x0[0], x0[1], f_hat, color='r') - + + +#%% +i = 0 +lowess = sm.nonparametric.lowess +x = sl[i][idx_beach,0] +y = sl[i][idx_beach,1] +sl_smooth = lowess(x,y, frac=1./15, it = 6) + +plt.figure() +plt.axis('equal') +plt.scatter +plt.plot(x,y,'bo-', linewidth=2, marker='o', + color='b', label='original') +plt.plot(sl_smooth[:,1], sl_smooth[:,0], linewidth=2, marker='o', + color='r', label='smooth') +plt.legend() +plt.xlabel('Eastings [m]') +plt.ylabel('Northings [m]') +plt.title('Local weighted scatterplot smoothing (LOWESS)') +plt.draw() +