valentinh
/
geetools_VH
forked from kilianv/CoastSat_WRL
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

update sds.py module

Browse Source
master
kvos 7 years ago
parent b964426cd4
commit 1882c98d9b
6 changed files with 84 additions and 549 deletions
Show Stats Download Patch File Download Diff File

6
.gitignore vendored
Unescape Escape View File

@ -0,0 +1,6 @@
*.pyc
*.mat
*.tif
*.png
*.mp4
*.gif

1
.gitingnore
Unescape Escape View File

@ -1 +0,0 @@
*.pyc

25
functions/sds.py
Unescape Escape View File

@ -25,6 +25,7 @@ import skimage.exposure as exposure
import skimage.transform as transform import skimage.transform as transform
import sklearn.decomposition as decomposition import sklearn.decomposition as decomposition
import skimage.measure as measure import skimage.measure as measure
import skimage.morphology as morphology
# import own modules # 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)] bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)]
return np.stack(bands, 2), georef 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 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 im_qa: np.ndarray
Image containing the QA band Image containing the QA band
satname: string
short name for the satellite (L8, L7, S2)
plot_bool: boolean
True if plot is wanted
Returns: Returns:
----------- -----------
@ -97,15 +102,27 @@ def create_cloud_mask(im_qa):
""" """
# convert QA bits # 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) 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_values = [2976, 2980, 2984, 2988, 3008, 3012, 3016, 3020]
#cloud_shadow_mask = np.isin(im_qa, cloud_shadow_values) #cloud_shadow_mask = np.isin(im_qa, cloud_shadow_values)
return cloud_mask 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) 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) 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]] qa_band = [im_bands[11]]
im_qa, crs_qa = load_image(im, polygon, qa_band) im_qa, crs_qa = load_image(im, polygon, qa_band)
im_qa = im_qa[:,:,0] 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]), im_cloud = transform.resize(im_cloud, (im_pan.shape[0], im_pan.shape[1]),
order=0, preserve_range=True, mode='constant').astype('bool_') order=0, preserve_range=True, mode='constant').astype('bool_')

5
functions/utils.py
Unescape Escape View File

@ -9,7 +9,6 @@ Contains all the utilities, convenience functions and small functions that do si
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import datetime
import pdb import pdb
@ -56,3 +55,7 @@ def compare_images(im1, im2):
def find_indices(lst, condition): def find_indices(lst, condition):
"imitation of MATLAB find function" "imitation of MATLAB find function"
return [i for i, elem in enumerate(lst) if condition(elem)] 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)]

237
sds_extract.py
Unescape Escape View File

@ -1,13 +1,10 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
"""
Created on Thu Mar 1 14:32:08 2018
@author: z5030440 #==========================================================#
# Extract shorelines from Landsat images
Main code to extract shorelines from Landsat imagery #==========================================================#
"""
# Preamble
# Initial settings
import ee import ee
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.cm as cm import matplotlib.cm as cm
@ -27,206 +24,78 @@ import sklearn.decomposition as decomposition
import skimage.morphology as morphology import skimage.morphology as morphology
import skimage.measure as measure import skimage.measure as measure
# my functions # my modules
import functions.utils as utils import functions.utils as utils
import functions.sds as sds import functions.sds as sds
# some settings
np.seterr(all='ignore') # raise/ignore divisions by 0 and nans 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() ee.Initialize()
#%% Select images
# parameters # parameters
plot_bool = False # if you want the plots cloud_thresh = 0.5 # threshold for cloud cover
prob_high = 99.9 # upper probability to clip and rescale pixel intensity plot_bool = True # if you want the plots
min_contour_points = 100 # minimum number of points contained in each water line prob_high = 99.9 # upper probability to clip and rescale pixel intensity
output_epsg = 28356 # GDA94 / MGA Zone 56 min_contour_points = 100# minimum number of points contained in each water line
cloud_threshold = 0.8 output_epsg = 28356 # GDA94 / MGA Zone 56
# select collection # 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) # location (Narrabeen-Collaroy beach)
rect_narra = [[[151.3473129272461,-33.69035274454718], polygon = [[[151.3473129272461,-33.69035274454718],
[151.2820816040039,-33.68206818063878], [151.2820816040039,-33.68206818063878],
[151.27281188964844,-33.74775138989556], [151.27281188964844,-33.74775138989556],
[151.3425064086914,-33.75231878701767], [151.3425064086914,-33.75231878701767],
[151.3473129272461,-33.69035274454718]]]; [151.3473129272461,-33.69035274454718]]];
with open('data/narra_beach.pkl', 'rb') as f: # dates
pts_beach = pickle.load(f) start_date = '2013-01-01'
end_date = '2018-12-31'
#rect_narra = [[[151.301454, -33.700754],
# [151.311453, -33.702075], # filter by location and date
# [151.307237, -33.739761], flt_col = input_col.filterBounds(ee.Geometry.Polygon(polygon)).filterDate(start_date, end_date)
# [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)
n_img = flt_col.size().getInfo() 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') im_all = flt_col.getInfo().get('features')
#%% Extract shorelines i = 0 # first image
metadata = {'timestamp':[],
'date_acquired':[], # find image in ee database
'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()
# 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
im = ee.Image(im_all[i].get('id')) 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_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) 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) 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.figure()
plt.imshow(im_ms_ps[:,:,[2,1,0]]) 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.axis('image')
plt.title('2016 shorelines') plt.title('Detected water lines')
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.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)

359
sds_extract_allcode.py
Unescape Escape View File

@ -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()
Write Preview
Loading…
Cancel
Save