# -*- coding: utf-8 -*- # Preamble import ee import matplotlib.pyplot as plt import matplotlib.cm as cm import numpy as np import pandas as pd from datetime import datetime import pickle import pdb import pytz from pylab import ginput import scipy.io as sio import scipy.interpolate import os # image processing modules import skimage.filters as filters import skimage.exposure as exposure import skimage.transform as transform import sklearn.decomposition as decomposition import skimage.morphology as morphology import skimage.measure as measure # my functions import functions.utils as utils import functions.sds as sds np.seterr(all='ignore') # raise/ignore divisions by 0 and nans ee.Initialize() au_tz = pytz.timezone('Australia/Sydney') #%% # load SDS shorelines with open('data\data_gt_l8.pkl', 'rb') as f: data = pickle.load(f) # load quadbike dates and convert from datenum to datetime suffix = '.mat' dir_name = os.getcwd() file_name = 'data\quadbike_dates' file_path = os.path.join(dir_name, file_name + suffix) quad_dates = sio.loadmat(file_path)['dates'] dt_quad = [] for i in range(quad_dates.shape[0]): dt_quad.append(datetime(quad_dates[i,0], quad_dates[i,1], quad_dates[i,2], tzinfo=au_tz)) # remove overlapping images, keep the one with lowest cloud_cover n = len(data['cloud_cover']) idx_worst = [] for i in range(n): date_im = data['date_acquired'][i] idx_double = np.isin(data['date_acquired'], date_im) if sum(idx_double.astype(int)) > 1: idx_worst.append(np.where(idx_double)[0][np.argmax(np.array(data['cloud_cover'])[idx_double])]) dt_sat = [] new_meta = {'contours':[], 'cloud_cover':[], 'geom_rmse_model':[], 'gcp_model':[], 'quality':[], 'sun_azimuth':[], 'sun_elevation':[]} for i in range(n): if not np.isin(i,idx_worst): dt_sat.append(data['dt'][i].astimezone(au_tz)) new_meta['contours'].append(data['contours'][i]) new_meta['cloud_cover'].append(data['cloud_cover'][i]) new_meta['geom_rmse_model'].append(data['geom_rmse_model'][i]) new_meta['gcp_model'].append(data['gcp_model'][i]) new_meta['quality'].append(data['quality'][i]) new_meta['sun_azimuth'].append(data['sun_azimuth'][i]) new_meta['sun_elevation'].append(data['sun_elevation'][i]) # calculate difference between days diff_days = [ [(x - _).days for _ in dt_quad] for x in dt_sat] day_thresh = 15 idx_close = [utils.find_indices(_, lambda e: abs(e) < day_thresh) for _ in diff_days] # put everything in a dictionnary and save it wl_comp = [] for i in range(len(dt_sat)): wl_comp.append({'sat dt': dt_sat[i], 'quad dt': [dt_quad[_] for _ in idx_close[i]], 'days diff': [diff_days[i][_] for _ in idx_close[i]], 'contours': new_meta['contours'][i], 'cloud_cover': new_meta['cloud_cover'][i], 'geom_rmse_model': new_meta['geom_rmse_model'][i], 'gcp_model': new_meta['gcp_model'][i], 'quality': new_meta['quality'][i], 'sun_azimuth': new_meta['sun_azimuth'][i], 'sun_elevation': new_meta['sun_elevation'][i]}) with open('wl_l8_comparison.pkl', 'wb') as f: pickle.dump(wl_comp, f) #%% with open('data\wl_l8_comparison.pkl', 'rb') as f: wl = pickle.load(f) # load quadbike dates and convert from datenum to datetime suffix = '.mat' dir_name = os.getcwd() subfolder_name = 'data\quadbike_surveys' file_path = os.path.join(dir_name, subfolder_name) file_names = os.listdir(file_path) for i in range(len(file_names)): fn_mat = os.path.join(file_path, file_names[i]) years = int(file_names[i][6:10]) months = int(file_names[i][11:13]) days = int(file_names[i][14:16]) for j in range(len(wl)): if wl[j]['quad dt'][0] == datetime(years, months, days, tzinfo=au_tz): quad_mat = sio.loadmat(fn_mat) wl[j].update({'quad_data':{'x':quad_mat['x'], 'y':quad_mat['y'], 'z':quad_mat['z'], 'dt': datetime(years, months, days, tzinfo=au_tz)}}) with open('data\wl_final.pkl', 'wb') as f: pickle.dump(wl, f) #%% with open('data\wl_final.pkl', 'rb') as f: wl = pickle.load(f) i = 0 x = wl[i]['quad_data']['x'] y = wl[i]['quad_data']['y'] z = wl[i]['quad_data']['z'] x = x.reshape(x.shape[0] * x.shape[1]) y = y.reshape(y.shape[0] * y.shape[1]) z = z.reshape(z.shape[0] * z.shape[1]) idx_nan = np.isnan(z) x_nan = x[idx_nan] y_nan = y[idx_nan] z_nan = z[idx_nan] x_nonan = x[~idx_nan] y_nonan = y[~idx_nan] z_nonan = z[~idx_nan] xs = x_nonan[::10] ys = y_nonan[::10] zs = z_nonan[::10] xq = wl[i]['contours'][:,0] yq = wl[i]['contours'][:,1] # cut xq around xs np.min(xs) np.max(xs) np.min(ys) np.max(ys) idx_x = np.logical_and(xq < np.max(xs), xq > np.min(xs)) idx_y = np.logical_and(yq < np.max(ys), yq > np.min(ys)) idx_in = np.logical_and(idx_x, idx_y) xq = xq[idx_in] yq = yq[idx_in] for i in range(len(xq)): idx_x = np.logical_and(xs < xq[i] + 10, xs > xq[i] - 10) idx_y = np.logical_and(ys < yq[i] + 10, ys > yq[i] - 10) xint = xs[idx_x] yint = ys[idx_y] f = interpolate.interp2d(xs, ys, zs, kind='linear') zq = f(xq,yq) plt.figure() plt.grid() plt.scatter(xs, ys, s=10, c=zs, marker='o', cmap=cm.get_cmap('jet'), label='quad data') plt.plot(xq,yq,'r-o', markersize=5, label='SDS') plt.axis('equal') plt.legend() plt.colorbar(label='mAHD') plt.xlabel('Eastings [m]') plt.ylabel('Northings [m]') plt.show() plt.figure() plt.plot(zq[:,0]) plt.show() plt.figure() plt.grid() plt.scatter(x_nonan, y_nonan, s=10, c=z_nonan, marker='o', cmap=cm.get_cmap('jet'), label='quad data') #plt.plot(x_nan, y_nan, 'k.', label='nans') plt.plot(xq,yq,'r-o', markersize=5, label='SDS') plt.axis('equal') plt.legend() plt.colorbar(label='mAHD') plt.xlabel('Eastings [m]') plt.ylabel('Northings [m]') plt.show() z2 = scipy.interpolate.griddata([x, y], z, [xq, yq], method='linear') f_interp = scipy.interpolate.interp2d(x1,y1,z1, kind='linear') sio.savemat('shoreline1.mat', {'x':xq, 'y':yq}) from scipy import interpolate x = np.arange(-5.01, 5.01, 0.01) y = np.arange(-5.01, 5.01, 0.01) xx, yy = np.meshgrid(x, y) z = np.sin(xx**2+yy**2) f = interpolate.interp2d(x, y, z, kind='cubic') xnew = np.arange(-5.01, 5.01, 1e-2) ynew = np.arange(-5.01, 5.01, 1e-2) znew = f(xnew, ynew) plt.plot(x, z[:, 0], 'ro-', xnew, znew[:, 0], 'b-') plt.show()