From 1882c98d9b615bb9c1712353375776224d45553a Mon Sep 17 00:00:00 2001 From: kvos Date: Mon, 9 Apr 2018 13:38:50 +1000 Subject: [PATCH] update sds.py module --- .gitignore | 6 + .gitingnore | 1 - functions/sds.py | 25 ++- functions/utils.py | 5 +- sds_extract.py | 237 ++++++--------------------- sds_extract_allcode.py | 359 ----------------------------------------- 6 files changed, 84 insertions(+), 549 deletions(-) create mode 100644 .gitignore delete mode 100644 .gitingnore delete mode 100644 sds_extract_allcode.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d088f68 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.pyc +*.mat +*.tif +*.png +*.mp4 +*.gif \ No newline at end of file diff --git a/.gitingnore b/.gitingnore deleted file mode 100644 index 7e99e36..0000000 --- a/.gitingnore +++ /dev/null @@ -1 +0,0 @@ -*.pyc \ No newline at end of file diff --git a/functions/sds.py b/functions/sds.py index a4cea6c..ec62f0e 100644 --- a/functions/sds.py +++ b/functions/sds.py @@ -25,6 +25,7 @@ import skimage.exposure as exposure import skimage.transform as transform import sklearn.decomposition as decomposition import skimage.measure as measure +import skimage.morphology as morphology # import own modules @@ -79,7 +80,7 @@ def load_image(image, polygon, bandsId): bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)] return np.stack(bands, 2), georef -def create_cloud_mask(im_qa): +def create_cloud_mask(im_qa, satname, plot_bool): """ Creates a cloud mask from the image containing the QA band information @@ -89,6 +90,10 @@ def create_cloud_mask(im_qa): ----------- im_qa: np.ndarray Image containing the QA band + satname: string + short name for the satellite (L8, L7, S2) + plot_bool: boolean + True if plot is wanted Returns: ----------- @@ -97,15 +102,27 @@ def create_cloud_mask(im_qa): """ # convert QA bits - cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908] + if satname == 'L8': + cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908] + elif satname == 'L7': + cloud_values = [752, 756, 760, 764] + cloud_mask = np.isin(im_qa, cloud_values) + # remove isolated cloud pixels (there are some in the swash and they cause problems) + if sum(sum(cloud_mask)) > 0: + morphology.remove_small_objects(cloud_mask, min_size=10, connectivity=1, in_place=True) + + if plot_bool: + plt.figure() + plt.imshow(cloud_mask, cmap='gray') + plt.draw() #cloud_shadow_values = [2976, 2980, 2984, 2988, 3008, 3012, 3016, 3020] #cloud_shadow_mask = np.isin(im_qa, cloud_shadow_values) return cloud_mask -def read_eeimage(im, polygon, plot_bool): +def read_eeimage(im, polygon, sat_name, plot_bool): """ Read an ee.Image() object and returns the panchromatic band, multispectral bands (B, G, R, NIR, SWIR) and a cloud mask. All outputs are at 15m resolution (bilinear interpolation for the multispectral bands) @@ -160,7 +177,7 @@ def read_eeimage(im, polygon, plot_bool): qa_band = [im_bands[11]] im_qa, crs_qa = load_image(im, polygon, qa_band) im_qa = im_qa[:,:,0] - im_cloud = create_cloud_mask(im_qa) + im_cloud = create_cloud_mask(im_qa, sat_name, plot_bool) im_cloud = transform.resize(im_cloud, (im_pan.shape[0], im_pan.shape[1]), order=0, preserve_range=True, mode='constant').astype('bool_') diff --git a/functions/utils.py b/functions/utils.py index 35aecca..019c0d4 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 @@ -56,3 +55,7 @@ def compare_images(im1, im2): def find_indices(lst, condition): "imitation of MATLAB find function" return [i for i, elem in enumerate(lst) if condition(elem)] + +def reject_outliers(data, m=2): + "rejects outliers in a numpy array" + return data[abs(data - np.mean(data)) < m * np.std(data)] diff --git a/sds_extract.py b/sds_extract.py index e04675d..8de4f40 100644 --- a/sds_extract.py +++ b/sds_extract.py @@ -1,13 +1,10 @@ # -*- coding: utf-8 -*- -""" -Created on Thu Mar 1 14:32:08 2018 -@author: z5030440 - -Main code to extract shorelines from Landsat imagery -""" -# Preamble +#==========================================================# +# Extract shorelines from Landsat images +#==========================================================# +# Initial settings import ee import matplotlib.pyplot as plt import matplotlib.cm as cm @@ -27,206 +24,78 @@ import sklearn.decomposition as decomposition import skimage.morphology as morphology import skimage.measure as measure -# my functions +# my modules 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'] = False +plt.rcParams['figure.max_open_warning'] = 100 ee.Initialize() -#%% Select images - # parameters -plot_bool = False # if you want the plots -prob_high = 99.9 # upper probability to clip and rescale pixel intensity -min_contour_points = 100 # minimum number of points contained in each water line -output_epsg = 28356 # GDA94 / MGA Zone 56 -cloud_threshold = 0.8 +cloud_thresh = 0.5 # threshold for cloud cover +plot_bool = True # if you want the plots +prob_high = 99.9 # upper probability to clip and rescale pixel intensity +min_contour_points = 100# minimum number of points contained in each water line +output_epsg = 28356 # GDA94 / MGA Zone 56 # select collection -input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA') +satname = 'L8' +input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA') # Landsat 8 Tier 1 TOA # location (Narrabeen-Collaroy beach) -rect_narra = [[[151.3473129272461,-33.69035274454718], - [151.2820816040039,-33.68206818063878], - [151.27281188964844,-33.74775138989556], - [151.3425064086914,-33.75231878701767], - [151.3473129272461,-33.69035274454718]]]; - -with open('data/narra_beach.pkl', 'rb') as f: - pts_beach = pickle.load(f) - -#rect_narra = [[[151.301454, -33.700754], -# [151.311453, -33.702075], -# [151.307237, -33.739761], -# [151.294220, -33.736329], -# [151.301454, -33.700754]]]; - -# Dates -start_date = '2016-01-01' -end_date = '2016-12-31' -# filter by location -flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra)).filterDate(start_date, end_date) +polygon = [[[151.3473129272461,-33.69035274454718], + [151.2820816040039,-33.68206818063878], + [151.27281188964844,-33.74775138989556], + [151.3425064086914,-33.75231878701767], + [151.3473129272461,-33.69035274454718]]]; + +# dates +start_date = '2013-01-01' +end_date = '2018-12-31' + +# filter by location and date +flt_col = input_col.filterBounds(ee.Geometry.Polygon(polygon)).filterDate(start_date, end_date) n_img = flt_col.size().getInfo() -print('Number of images covering Narrabeen:', n_img) +print('Number of images covering the polygon:', n_img) im_all = flt_col.getInfo().get('features') -#%% Extract shorelines -metadata = {'timestamp':[], - 'date_acquired':[], - 'cloud_cover':[], - 'geom_rmse_model':[], - 'gcp_model':[], - 'quality':[], - 'sun_azimuth':[], - 'sun_elevation':[]} -skipped_images = np.zeros((n_img,1)).astype(bool) -output_wl = [] -# loop through all images -for i in range(n_img): - - # find each image in ee database - im = ee.Image(im_all[i].get('id')) - # load image as np.array - im_pan, im_ms, im_cloud, crs, meta = sds.read_eeimage(im, rect_narra, plot_bool) - - # if clouds -> skip the image - if sum(sum(im_cloud.astype(int)))/(im_cloud.shape[0]*im_cloud.shape[1]) > cloud_threshold: - skipped_images[i] = True - continue - - # rescale intensities - im_ms = sds.rescale_image_intensity(im_ms, im_cloud, prob_high, plot_bool) - im_pan = sds.rescale_image_intensity(im_pan, im_cloud, prob_high, plot_bool) - - # pansharpen rgb image - im_ms_ps = sds.pansharpen(im_ms[:,:,[0,1,2]], im_pan, im_cloud, plot_bool) - - # add down-sized bands for NIR and SWIR (since pansharpening is not possible) - im_ms_ps = np.append(im_ms_ps, im_ms[:,:,[3,4]], axis=2) - - # calculate NDWI - im_ndwi = sds.nd_index(im_ms_ps[:,:,3], im_ms_ps[:,:,1], im_cloud, plot_bool) - - # edge detection - wl_pix = sds.find_wl_contours(im_ndwi, im_cloud, min_contour_points, plot_bool) - - plt.figure() - plt.imshow(im_ms_ps[:,:,[2,1,0]]) - for i,contour in enumerate(wl_pix): plt.plot(contour[:, 1], contour[:, 0], linewidth=2) - plt.axis('image') - plt.title('Detected water lines') - plt.show() - - # convert from pixels to world coordinates - wl_coords = sds.convert_pix2world(wl_pix, crs['crs_15m']) - # convert to output epsg spatial reference - wl = sds.convert_epsg(wl_coords, crs['epsg_code'], output_epsg) - - # find contour closest to narrabeen beach - sum_dist = np.zeros(len(wl)) - for k,contour in enumerate(wl): - min_dist = np.zeros(len(pts_beach)) - for j,pt in enumerate(pts_beach): - min_dist[j] = np.min(np.linalg.norm(contour - pt, axis=1)) - sum_dist[k] = np.sum(min_dist)/len(min_dist) - try: - wl_beach = wl[np.argmin(sum_dist)] -# plt.figure() -# plt.axis('equal') -# plt.plot(pts_beach[:,0], pts_beach[:,1], 'ko') -# plt.plot(wl_beach[:,0], wl_beach[:,1], 'r') -# plt.show() - except: - wl_beach = [] - - # plot for QA - plt.figure() - plt.imshow(im_ms_ps[:,:,[2,1,0]]) - for k,contour in enumerate(wl_pix): plt.plot(contour[:, 1], contour[:, 0], linewidth=2) - if len(wl_beach) > 0: - plt.plot(wl_pix[np.argmin(sum_dist)][:,1], wl_pix[np.argmin(sum_dist)][:,0], linewidth=3, color='w') - plt.axis('image') - plt.title('im ' + str(i) + ' : ' + datetime.strftime(datetime - .fromtimestamp(meta['timestamp']/1000, tz=pytz.utc) - .astimezone(pytz.timezone('Australia/Sydney')), '%Y-%m-%d %H:%M:%S %Z%z')) - plt.show() +i = 0 # first image - # store metadata of each image in dict - metadata['timestamp'].append(meta['timestamp']) - metadata['date_acquired'].append(meta['date_acquired']) - metadata['cloud_cover'].append(sum(sum(im_cloud.astype(int)))/(im_cloud.shape[0]*im_cloud.shape[1])) - metadata['geom_rmse_model'].append(meta['geom_rmse_model']) - metadata['gcp_model'].append(meta['gcp_model']) - metadata['quality'].append(meta['quality']) - metadata['sun_azimuth'].append(meta['sun_azimuth']) - metadata['sun_elevation'].append(meta['sun_elevation']) - # store water lines - output_wl.append(wl_beach) - - print(i) - -# generate datetimes -#fmt = '%Y-%m-%d %H:%M:%S %Z%z' -#au_tz = pytz.timezone('Australia/Sydney') -dt = []; -t = metadata['timestamp'] -for k in range(len(t)): dt.append(datetime.fromtimestamp(t[k]/1000, tz=pytz.utc)) - -# save outputs -data = metadata.copy() -data.update({'dt':dt}) -data.update({'contours':output_wl}) - -#with open('data_2016.pkl', 'wb') as f: -# pickle.dump(data, f) -#%% Load data - -#with open('data_2016.pkl', 'rb') as f: -# data = pickle.load(f) - - -# load backgroud image -i = 0 +# find image in ee database im = ee.Image(im_all[i].get('id')) -im_pan, im_ms, im_cloud, crs, meta = sds.read_eeimage(im, rect_narra, plot_bool) + +# load image as np.array +im_pan, im_ms, im_cloud, crs, meta = sds.read_eeimage(im, polygon, satname, plot_bool) + +# rescale intensities im_ms = sds.rescale_image_intensity(im_ms, im_cloud, prob_high, plot_bool) im_pan = sds.rescale_image_intensity(im_pan, im_cloud, prob_high, plot_bool) + +# pansharpen rgb image im_ms_ps = sds.pansharpen(im_ms[:,:,[0,1,2]], im_pan, im_cloud, plot_bool) +# add down-sized bands for NIR and SWIR (since pansharpening is not possible) +im_ms_ps = np.append(im_ms_ps, im_ms[:,:,[3,4]], axis=2) + +# calculate NDWI +im_ndwi = sds.nd_index(im_ms_ps[:,:,3], im_ms_ps[:,:,1], im_cloud, plot_bool) + +# edge detection +wl_pix = sds.find_wl_contours(im_ndwi, im_cloud, min_contour_points, plot_bool) + plt.figure() plt.imshow(im_ms_ps[:,:,[2,1,0]]) +for i,contour in enumerate(wl_pix): plt.plot(contour[:, 1], contour[:, 0], linewidth=2) plt.axis('image') -plt.title('2016 shorelines') - -n = len(data['cloud_cover']) -idx_best = [] -# remove overlapping images, based on cloud cover -for i in range(n): - date_im = data['date_acquired'][i] - idx = np.isin(data['date_acquired'], date_im) - best = np.where(idx)[0][np.argmin(np.array(data['cloud_cover'])[idx])] - if ~np.isin(best, idx_best): - idx_best.append(best) - -point_narra = np.array([342500, 6266990]) -plt.figure() -plt.axis('equal') -plt.grid() -cmap = cm.get_cmap('jet') -colours = cmap(np.linspace(0, 1, num=len(idx_best))) -for i, idx in enumerate(idx_best): - for j in range(len(data['contours'][i])): - if np.any(np.linalg.norm(data['contours'][i][j][:,[0,1]] - point_narra, axis=1) < 200): - plt.plot(data['contours'][i][j][:,0], data['contours'][i][j][:,1], - label=str(data['date_acquired'][i]), - linewidth=2, color=colours[i,:]) - -plt.legend() -plt.show() +plt.title('Detected water lines') +plt.show() - - - +# convert from pixels to world coordinates +wl_coords = sds.convert_pix2world(wl_pix, crs['crs_15m']) +# convert to output epsg spatial reference +wl = sds.convert_epsg(wl_coords, crs['epsg_code'], output_epsg) diff --git a/sds_extract_allcode.py b/sds_extract_allcode.py deleted file mode 100644 index 9ef6659..0000000 --- a/sds_extract_allcode.py +++ /dev/null @@ -1,359 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Feb 21 18:05:01 2018 - -@author: z5030440 -""" - -#%% Initial settings - -# import packages -import ee -from IPython import display -import math -import matplotlib.pyplot as plt -import numpy as np -from osgeo import gdal -import tempfile -import tensorflow as tf -import urllib -from urllib.request import urlretrieve -import zipfile -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 scripts -from GEEImageFunctions import * - -np.seterr(all='ignore') # raise divisions by 0 and nans -ee.Initialize() - -# Load image collection and filter it based on location (Narrabeen) - -input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA') -#n_img = input_col.size().getInfo() -#print('Number of images in collection:', n_img) - -# filter based on location (Narrabeen-Collaroy) -rect_narra = [[[151.3473129272461,-33.69035274454718], - [151.2820816040039,-33.68206818063878], - [151.27281188964844,-33.74775138989556], - [151.3425064086914,-33.75231878701767], - [151.3473129272461,-33.69035274454718]]]; - -flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra)) -n_img = flt_col.size().getInfo() -print('Number of images covering Narrabeen:', n_img) - -# Select the most recent image and download it - -im = ee.Image(flt_col.sort('SENSING_TIME',False).first()) -im_dic = im.getInfo() -image_prop = im_dic.get('properties') -im_bands = im_dic.get('bands') -for i in range(len(im_bands)): del im_bands[i]['dimensions'] # delete the dimensions key - -# download the panchromatic band (B8) -pan_band = [im_bands[7]] -im_pan = load_image(im, rect_narra, pan_band) -im_pan = im_pan[:,:,0] -size_pan = im_pan.shape -vec_pan = im_pan.reshape(size_pan[0] * size_pan[1]) -# download the QA band (BQA) -qa_band = [im_bands[11]] -im_qa = load_image(im, rect_narra, qa_band) -im_qa = im_qa[:,:,0] - -# convert QA bits -cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908] -cloud_shadow_values = [2976, 2980, 2984, 2988, 3008, 3012, 3016, 3020] - -# Create cloud mask (resized to be applied to the Pan band) -im_cloud = np.isin(im_qa, cloud_values) -im_cloud_shadow = np.isin(im_qa, cloud_shadow_values) -im_cloud_res = transform.resize(im_cloud,(im_pan.shape[0], im_pan.shape[1]), order=0, preserve_range=True).astype('bool_') -vec_cloud = im_cloud.reshape(im_cloud.shape[0] * im_cloud.shape[1]) -vec_cloud_res = im_cloud_res.reshape(size_pan[0] * size_pan[1]) - - -# download the other bands (B2,B3,B4,B5,B6) = (blue,green,red,nir,swir1) -ms_bands = [im_bands[1], im_bands[2], im_bands[3], im_bands[4], im_bands[5]] -im_ms = load_image(im, rect_narra, ms_bands) -size_ms = im_ms.shape -vec_ms = im_ms.reshape(size_ms[0] * size_ms[1], size_ms[2]) - -# Plot the RGB image and cloud masks -plt.figure() -ax1 = plt.subplot(121) -plt.imshow(im_ms[:,:,[2,1,0]]) -plt.title('RGB') -ax2 = plt.subplot(122) -plt.imshow(im_cloud, cmap='gray') -plt.title('Cloud mask') -#ax3 = plt.subplot(133, sharex=ax1, sharey=ax1) -#plt.imshow(im_cloud_shadow) -#plt.title('Cloud mask shadow') -plt.show() - -# Resize multispectral bands (30m) to the size of the pan band (15m) using bilinear interpolation -im_ms_res = transform.resize(im_ms,(size_pan[0], size_pan[1]), order=1, preserve_range=True, mode='constant') -vec_ms_res = im_ms_res.reshape(size_pan[0] * size_pan[1], size_ms[2]) - -# Adjust intensities (set cloud pixels to 0 intensity) -cloud_value = np.nan - -prc_low = 0 # lower percentile -prob_high = 99.9 # upper percentile probability to clip - -# Rescale intensities between 0 and 1 -vec_ms_adj = np.ones((len(vec_cloud_res),size_ms[2])) * np.nan -for i in range(im_ms.shape[2]): - prc_high = np.percentile(vec_ms_res[~vec_cloud_res,i], prob_high) - vec_rescaled = exposure.rescale_intensity(vec_ms_res[~vec_cloud_res,i], in_range=(prc_low,prc_high)) - plt.figure() - plt.hist(vec_rescaled, bins = 300) - plt.show() - vec_ms_adj[~vec_cloud_res,i] = vec_rescaled - - -im_ms_adj = vec_ms_adj.reshape(size_pan[0], size_pan[1], size_ms[2]) - -# same for the pan band -vec_pan_adj = np.ones(len(vec_cloud_res)) * np.nan -prc_high = np.percentile(vec_pan[~vec_cloud_res],prob_high) -vec_rescaled = exposure.rescale_intensity(vec_pan[~vec_cloud_res], in_range=(prc_low,prc_high)) -plt.figure() -plt.hist(vec_rescaled, bins = 300) -plt.show() -vec_pan_adj[~vec_cloud_res] = vec_rescaled -im_pan_adj = vec_pan_adj.reshape(size_pan[0], size_pan[1]) - -# Plot adjusted images -plt.figure() -plt.subplot(131) -plt.imshow(im_pan_adj, cmap='gray') -plt.title('PANCHROMATIC (15 m pixel)') - -plt.subplot(132) -plt.imshow(im_ms_adj[:,:,[2,1,0]]) -plt.title('RGB (30 m pixel)') -plt.show() - -plt.subplot(133) -plt.imshow(im_ms_adj[:,:,[3,1,0]]) -plt.title('NIR-GB (30 m pixel)') -plt.show() - - -#%% Pansharpening (PCA) -# Run PCA on selected bands - -sel_bands = [0,1,2] -temp = vec_ms_adj[:,sel_bands] -vec_ms_adj_nocloud = temp[~vec_cloud_res,:] -pca = decomposition.PCA() -vec_pcs = pca.fit_transform(vec_ms_adj_nocloud) -vec_pcs_all = np.ones((len(vec_cloud_res),len(sel_bands))) * np.nan -vec_pcs_all[~vec_cloud_res,:] = vec_pcs -im_pcs = vec_pcs_all.reshape(size_pan[0], size_pan[1], vec_pcs.shape[1]) - -plt.figure() -plt.subplot(221) -plt.imshow(im_pcs[:,:,0], cmap='gray') -plt.title('Component 1') -plt.subplot(222) -plt.imshow(im_pcs[:,:,1], cmap='gray') -plt.title('Component 2') -plt.subplot(223) -plt.imshow(im_pcs[:,:,2], cmap='gray') -plt.title('Component 3') -plt.show() - -# Compare the Pan image with the 1st Principal component -compare_images(im_pan_adj,im_pcs[:,:,0]) -intensity_histogram(im_pan_adj) -intensity_histogram(im_pcs[:,:,0]) - -# Match histogram of the pan image with the 1st principal component and replace the 1st component -vec_pcs[:,0] = hist_match(vec_pan_adj[~vec_cloud_res], vec_pcs[:,0]) -vec_ms_ps = pca.inverse_transform(vec_pcs) - -# normalise between 0 and 1 -for i in range(vec_pcs.shape[1]): - vec_ms_ps[:,i] = np.divide(vec_ms_ps[:,i] - np.min(vec_ms_ps[:,i]), - np.max(vec_ms_ps[:,i]) - np.min(vec_ms_ps[:,i])) - -vec_ms_ps_all = np.ones((len(vec_cloud_res),len(sel_bands))) * np.nan -vec_ms_ps_all[~vec_cloud_res,:] = vec_ms_ps -im_ms_ps = vec_ms_ps_all.reshape(size_pan[0], size_pan[1], len(sel_bands)) -vec_ms_ps_all = np.append(vec_ms_ps_all, vec_ms_adj[:,[3,4]], axis=1) -im_ms_ps = np.append(im_ms_ps, im_ms_adj[:,:,[3,4]], axis=2) - -# Plot adjusted images -plt.figure() -plt.subplot(121) -plt.imshow(im_ms_adj[:,:,[2,1,0]]) -plt.title('Original RGB') -plt.show() - -plt.subplot(122) -plt.imshow(im_ms_ps[:,:,[2,1,0]]) -plt.title('Pansharpened RGB') -plt.show() - -plt.figure() -plt.subplot(121) -plt.imshow(im_ms_adj[:,:,[3,1,0]]) -plt.title('Original NIR-GB') -plt.show() - -plt.subplot(122) -plt.imshow(im_ms_ps[:,:,[3,1,0]]) -plt.title('Pansharpened NIR-GB') -plt.show() -#%% Compute Normalized Difference Water Index (NDWI) - -# With NIR -vec_ndwi_nir = np.ones(len(vec_cloud_res)) * np.nan -temp = np.divide(vec_ms_ps_all[~vec_cloud_res,3] - vec_ms_ps_all[~vec_cloud_res,1], - vec_ms_ps_all[~vec_cloud_res,3] + vec_ms_ps_all[~vec_cloud_res,1]) -vec_ndwi_nir[~vec_cloud_res] = temp -im_ndwi_nir = vec_ndwi_nir.reshape(size_pan[0], size_pan[1]) - -# With SWIR_1 -vec_ndwi_swir = np.ones(len(vec_cloud_res)) * np.nan -temp = np.divide(vec_ms_ps_all[~vec_cloud_res,4] - vec_ms_ps_all[~vec_cloud_res,1], - vec_ms_ps_all[~vec_cloud_res,4] + vec_ms_ps_all[~vec_cloud_res,1]) -vec_ndwi_swir[~vec_cloud_res] = temp -im_ndwi_swir = vec_ndwi_swir.reshape(size_pan[0], size_pan[1]) - -plt.figure() -ax1 = plt.subplot(211) -plt.hist(vec_ndwi_nir[~vec_cloud_res], bins=300, label='NIR') -plt.hist(vec_ndwi_swir[~vec_cloud_res], bins=300, label='SWIR', alpha=0.5) -plt.legend() -ax2 = plt.subplot(212, sharex=ax1) -plt.hist(vec_ndwi_nir[~vec_cloud_res], bins=300, cumulative=True, histtype='step', label='NIR') -plt.hist(vec_ndwi_swir[~vec_cloud_res], bins=300, cumulative=True, histtype='step', label='SWIR') -plt.legend() -plt.show() - -compare_images(im_ndwi_nir,im_ndwi_swir) - -plt.figure() -plt.imshow(im_ndwi_nir, cmap='seismic') -plt.title('Water Index') -plt.colorbar() -plt.show() - -#%% Extract shorelines (NIR) - -ndwi_nir = vec_ndwi_nir[~vec_cloud_res] - -t_otsu = filters.threshold_otsu(ndwi_nir) -t_min = filters.threshold_minimum(ndwi_nir) -t_mean = filters.threshold_mean(ndwi_nir) -t_li = filters.threshold_li(ndwi_nir) -# try all thresholding algorithms - -plt.figure() -plt.hist(ndwi_nir, bins=300) -plt.plot([t_otsu, t_otsu],[0, 15000], 'r-', label='Otsu threshold') -#plt.plot([t_min, t_min],[0, 15000], 'g-', label='min') -#plt.plot([t_mean, t_mean],[0, 15000], 'y-', label='mean') -#plt.plot([t_li, t_li],[0, 15000], 'm-', label='li') -plt.legend() -plt.show() - -plt.figure() -plt.imshow(im_ndwi_nir > t_otsu, cmap='gray') -plt.title('Binary image') -plt.show() - -im_bin = im_ndwi_nir > t_otsu -im_open = morphology.binary_opening(im_bin,morphology.disk(1)) -im_close = morphology.binary_closing(im_open,morphology.disk(1)) -im_bin_coast_in = im_close ^ morphology.erosion(im_close,morphology.disk(1)) -im_bin_sl_in = morphology.remove_small_objects(im_bin_coast_in,100,8) - -compare_images(im_close,im_bin_sl_in) - -plt.figure() -plt.subplot(121) -plt.imshow(im_close, cmap='gray') -plt.title('morphological closing') -plt.subplot(122) -plt.imshow(im_bin_sl_in, cmap='gray') -plt.title('Water mark') -plt.show() - - -im_bin_coast_out = morphology.dilation(im_close,morphology.disk(1)) ^ im_close -im_bin_sl_out = morphology.remove_small_objects(im_bin_coast_out,100,8) - - -# Plot shorelines on top of RGB image -im_rgb_sl = np.copy(im_ms_ps[:,:,[2,1,0]]) - -im_rgb_sl[im_bin_sl_in,0] = 0 -im_rgb_sl[im_bin_sl_in,1] = 1 -im_rgb_sl[im_bin_sl_in,2] = 1 - -im_rgb_sl[im_bin_sl_out,0] = 1 -im_rgb_sl[im_bin_sl_out,1] = 0 -im_rgb_sl[im_bin_sl_out,2] = 1 - -plt.figure() -plt.imshow(im_rgb_sl) -plt.title('Pansharpened RGB') -plt.show() - -#%% Extract shorelines SWIR - -ndwi_swir = vec_ndwi_swir[~vec_cloud_res] -t_otsu = filters.threshold_otsu(ndwi_swir) - -plt.figure() -plt.hist(ndwi_swir, bins=300) -plt.plot([t_otsu, t_otsu],[0, 15000], 'r-', label='Otsu threshold') -#plt.plot([t_min, t_min],[0, 15000], 'g-', label='min') -#plt.plot([t_mean, t_mean],[0, 15000], 'y-', label='mean') -#plt.plot([t_li, t_li],[0, 15000], 'm-', label='li') -plt.legend() -plt.show() - -plt.figure() -plt.imshow(im_ndwi_swir > t_otsu, cmap='gray') -plt.title('Binary image') -plt.show() - -im_bin = im_ndwi_swir > t_otsu -im_open = morphology.binary_opening(im_bin,morphology.disk(1)) -im_close = morphology.binary_closing(im_open,morphology.disk(1)) -im_bin_coast_in = im_close ^ morphology.erosion(im_close,morphology.disk(1)) -im_bin_sl_in = morphology.remove_small_objects(im_bin_coast_in,100,8) - -compare_images(im_close,im_bin_sl_in) - - -im_bin_coast_out = morphology.dilation(im_close,morphology.disk(1)) ^ im_close -im_bin_sl_out = morphology.remove_small_objects(im_bin_coast_out,100,8) - - -# Plot shorelines on top of RGB image -im_rgb_sl = np.copy(im_ms_ps[:,:,[2,1,0]]) - -im_rgb_sl[im_bin_sl_in,0] = 0 -im_rgb_sl[im_bin_sl_in,1] = 1 -im_rgb_sl[im_bin_sl_in,2] = 1 - -im_rgb_sl[im_bin_sl_out,0] = 1 -im_rgb_sl[im_bin_sl_out,1] = 0 -im_rgb_sl[im_bin_sl_out,2] = 1 - -plt.figure() -plt.imshow(im_rgb_sl) -plt.title('Pansharpened RGB') -plt.show() \ No newline at end of file