CoastSat_WRL/sds_extract.py

# -*- coding: utf-8 -*-
"""
Created on Thu Mar  1 14:32:08 2018

@author: z5030440

Main code to extract shorelines from Landsat imagery
"""
# Preamble

import ee
from IPython import display
import math
import matplotlib.pyplot as plt
import numpy as np
import pdb
# 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
# my modules
from utils import *
from sds1 import *

np.seterr(all='ignore') # raise/ignore divisions by 0 and nans
ee.Initialize()
plot_bool = True

input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA')
# filter collection on location (Narrabeen-Collaroy rectangle)
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 of the filtered collection
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) and QA band (Q11)
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])

qa_band = [im_bands[11]]
im_qa = load_image(im, rect_narra, qa_band)
im_qa = im_qa[:,:,0]

# 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])

# create cloud mask
im_cloud = create_cloud_mask(im_qa)
im_cloud_res = transform.resize(im_cloud, (size_pan[0], size_pan[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])

# 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')
# Rescale image intensity between 0 and 1
prob_high = 99.9
im_ms_adj = rescale_image_intensity(im_ms_res, im_cloud_res, prob_high, plot_bool)
im_pan_adj = rescale_image_intensity(im_pan, im_cloud_res, prob_high, plot_bool)

# 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

im_ms_ps = pansharpen(im_ms_adj[:,:,[0,1,2]], im_pan_adj, im_cloud_res)
# Add resized bands for NIR and SWIR
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()

im_ndwi_nir = normalized_difference(im_ms_ps[:,:,3], im_ms_ps[:,:,1], im_cloud_res, plot_bool)
vec_ndwi_nir = im_ndwi_nir.reshape(size_pan[0] * size_pan[1])

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)

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()