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.
222 lines
6.5 KiB
Python
222 lines
6.5 KiB
Python
# -*- 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() |