From 647ef2192e5f1a70dc432925708c508dc470f0e4 Mon Sep 17 00:00:00 2001 From: kvos Date: Fri, 16 Mar 2018 14:27:58 +1100 Subject: [PATCH] initial shoreline mapping Loads images from the Google Earth Engine database and maps shorelines using the marching squares algorithm --- geocheck.py | 221 ++++++++++++++++ sds.py | 519 ++++++++++++++++++++++++++++++++++++++ sds_extract.py | 200 +++++++++++++++ sds_extract_collection.py | 76 ++++++ shoreline_extraction.py | 359 ++++++++++++++++++++++++++ utils.py | 54 ++++ 6 files changed, 1429 insertions(+) create mode 100644 geocheck.py create mode 100644 sds.py create mode 100644 sds_extract.py create mode 100644 sds_extract_collection.py create mode 100644 shoreline_extraction.py create mode 100644 utils.py diff --git a/geocheck.py b/geocheck.py new file mode 100644 index 0000000..b39d3d8 --- /dev/null +++ b/geocheck.py @@ -0,0 +1,221 @@ +# -*- 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 +import skimage.measure as measure + +from shapely.geometry import Polygon + +from osgeo import gdal +from osgeo import osr +import tempfile +import urllib +from urllib.request import urlretrieve +import zipfile + +# my modules +from utils import * +# from sds import * + +np.seterr(all='ignore') # raise/ignore divisions by 0 and nans +ee.Initialize() +plot_bool = True # if you want the plots + +def download_tif(image, bandsId): + """downloads tif image (region and bands) from the ee server and stores it in a temp file""" + url = ee.data.makeDownloadUrl(ee.data.getDownloadId({ + 'image': image.serialize(), + 'bands': bandsId, + 'filePerBand': 'false', + 'name': 'data', + })) + local_zip, headers = urlretrieve(url) + with zipfile.ZipFile(local_zip) as local_zipfile: + return local_zipfile.extract('data.tif', tempfile.mkdtemp()) + +def load_image(image, bandsId): + """loads an ee.Image() as a np.array. e.Image() is retrieved from the EE database.""" + local_tif_filename = download_tif(image, bandsId) + dataset = gdal.Open(local_tif_filename, gdal.GA_ReadOnly) + bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)] + return np.stack(bands, 2), dataset + + + +im = ee.Image('LANDSAT/LC08/C01/T1_RT_TOA/LC08_089083_20130411') + +lon = [151.2820816040039, 151.3425064086914] +lat = [-33.68206818063878, -33.74775138989556] +polygon = [[lon[0], lat[0]], [lon[1], lat[0]], [lon[1], lat[1]], [lon[0], lat[1]]]; + +# get image metadata into dictionnary +im_dic = im.getInfo() +im_bands = im_dic.get('bands') +# delete dimensions key from dictionnary, otherwise the entire image is extracted +#for i in range(len(im_bands)): del im_bands[i]['dimensions'] +pan_band = [im_bands[7]] +ms_bands = [im_bands[1], im_bands[2], im_bands[3]] +im_full, dataset_full = load_image(im, ms_bands) + +plt.figure() +plt.imshow(np.clip(im_full[:,:,[2,1,0]] * 3, 0, 1)) +plt.show() +#%% + +def download_tif(image, polygon, bandsId): + """downloads tif image (region and bands) from the ee server and stores it in a temp file""" + url = ee.data.makeDownloadUrl(ee.data.getDownloadId({ + 'image': image.serialize(), + 'region': polygon, + 'bands': bandsId, + 'filePerBand': 'false', + 'name': 'data', + })) + local_zip, headers = urlretrieve(url) + with zipfile.ZipFile(local_zip) as local_zipfile: + return local_zipfile.extract('data.tif', tempfile.mkdtemp()) + +def load_image(image, polygon, bandsId): + """ + Loads an ee.Image() as a np.array. e.Image() is retrieved from the EE database. + The geographic area and bands to select can be specified + + KV WRL 2018 + + Arguments: + ----------- + image: ee.Image() + image objec from the EE database + polygon: list + coordinates of the points creating a polygon. Each point is a list with 2 values + bandsId: list + bands to select, each band is a dictionnary in the list containing the following keys: + crs, crs_transform, data_type and id. NOTE: you have to remove the key dimensions, otherwise + the entire image is retrieved. + + Returns: + ----------- + image_array : np.ndarray + An array containing the image (2D if one band, otherwise 3D) + """ + + local_tif_filename = download_tif(image, polygon, bandsId) + dataset = gdal.Open(local_tif_filename, gdal.GA_ReadOnly) + bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)] + return np.stack(bands, 2), dataset + + +for i in range(len(im_bands)): del im_bands[i]['dimensions'] +ms_bands = [im_bands[1], im_bands[2], im_bands[3]] + +im_cropped, dataset_cropped = load_image(im, polygon, ms_bands) + +plt.figure() +plt.imshow(np.clip(im_cropped[:,:,[2,1,0]] * 3, 0, 1)) +plt.show() + +#%% +crs_full = dataset_full.GetGeoTransform() +crs_cropped = dataset_cropped.GetGeoTransform() +scale = crs_full[1] +ul_full = np.array([crs_full[0], crs_full[3]]) +ul_cropped = np.array([crs_cropped[0], crs_cropped[3]]) + +delta = np.abs(ul_full - ul_cropped)/scale + +u0 = delta[0].astype('int') +v0 = delta[1].astype('int') + +im_full[v0,u0,:] +im_cropped[0,0,:] + +lrx = ul_cropped[0] + (dataset_cropped.RasterXSize * scale) +lry = ul_cropped[1] + (dataset_cropped.RasterYSize * (-scale)) + +lr_cropped = np.array([lrx, lry]) + +delta = np.abs(ul_full - lr_cropped)/scale +u1 = delta[0].astype('int') +v1 = delta[1].astype('int') + +im_cropped2 = im_full[v0:v1,u0:u1,:] + +#%% +crs_full = dataset_full.GetGeoTransform() +source = osr.SpatialReference() +source.ImportFromWkt(dataset_full.GetProjection()) + +target = osr.SpatialReference() +target.ImportFromEPSG(4326) + +transform = osr.CoordinateTransformation(source, target) + +transform.TransformPoint(ulx, uly) + +#%% +crs_cropped = dataset_cropped.GetGeoTransform() +ulx = crs_cropped[0] +uly = crs_cropped[3] +source = osr.SpatialReference() +source.ImportFromWkt(dataset_cropped.GetProjection()) + +target = osr.SpatialReference() +target.ImportFromEPSG(4326) + +transform = osr.CoordinateTransformation(source, target) + +transform.TransformPoint(lrx, lry) + + + +#%% +source = osr.SpatialReference() +source.ImportFromEPSG(4326) + +target = osr.SpatialReference() +target.ImportFromEPSG(32656) + +coords = transform.TransformPoint(151.2820816040039, -33.68206818063878) +coords[0] - ulx +coords[1] - uly +#%% +x_ul_full = ms_bands[0]['crs_transform'][2] +y_ul_full = ms_bands[0]['crs_transform'][5] +scale = ms_bands[0]['crs_transform'][0] + +x_ul_cropped = np.array([340756.105840223, 346357.851288875, 346474.839525944, 340877.362938763]) +y_ul_cropped = np.array([-3728229.45372866, -3728137.91775723, -3735421.58347927, -3735513.20696522]) + +dx = abs(x_ul_full - x_ul_cropped) +dy = abs(y_ul_full - y_ul_cropped) + +u_coord = np.round(dx/scale).astype('int') +v_coord = np.round(dy/scale).astype('int') + +im_cropped2 = im_full[np.min(v_coord):np.max(v_coord), np.min(u_coord):np.max(u_coord),:] + +plt.figure() +plt.imshow(np.clip(im_cropped2[:,:,[2,1,0]] * 3, 0, 1), cmap='gray') +plt.show() + +sum(sum(sum(np.equal(im_cropped,im_cropped2).astype('int')-1))) + diff --git a/sds.py b/sds.py new file mode 100644 index 0000000..799e409 --- /dev/null +++ b/sds.py @@ -0,0 +1,519 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Mar 1 11:20:35 2018 + +@author: z5030440 +""" + +"""This script contains the functions needed for satellite derived shoreline (SDS) extraction""" + +# Initial settings +import numpy as np +import matplotlib.pyplot as plt +import pdb +import ee + +# other modules +from osgeo import gdal +import tempfile +from urllib.request import urlretrieve +import zipfile + +# 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.measure as measure + + +# import own modules +from utils import * + + +# Download from ee server function + +def download_tif(image, polygon, bandsId): + """downloads tif image (region and bands) from the ee server and stores it in a temp file""" + url = ee.data.makeDownloadUrl(ee.data.getDownloadId({ + 'image': image.serialize(), + 'region': polygon, + 'bands': bandsId, + 'filePerBand': 'false', + 'name': 'data', + })) + local_zip, headers = urlretrieve(url) + with zipfile.ZipFile(local_zip) as local_zipfile: + return local_zipfile.extract('data.tif', tempfile.mkdtemp()) + +def load_image(image, polygon, bandsId): + """ + Loads an ee.Image() as a np.array. e.Image() is retrieved from the EE database. + The geographic area and bands to select can be specified + + KV WRL 2018 + + Arguments: + ----------- + image: ee.Image() + image objec from the EE database + polygon: list + coordinates of the points creating a polygon. Each point is a list with 2 values + bandsId: list + bands to select, each band is a dictionnary in the list containing the following keys: + crs, crs_transform, data_type and id. NOTE: you have to remove the key dimensions, otherwise + the entire image is retrieved. + + Returns: + ----------- + image_array : np.ndarray + An array containing the image (2D if one band, otherwise 3D) + georef : np.ndarray + 6 element vector containing the crs_parameters + [X_ul_corner Xscale Xshear Y_ul_corner Yshear Yscale] + """ + + local_tif_filename = download_tif(image, polygon, bandsId) + dataset = gdal.Open(local_tif_filename, gdal.GA_ReadOnly) + georef = np.array(dataset.GetGeoTransform()) + bands = [dataset.GetRasterBand(i + 1).ReadAsArray() for i in range(dataset.RasterCount)] + return np.stack(bands, 2), georef + +def create_cloud_mask(im_qa): + """ + Creates a cloud mask from the image containing the QA band information + + KV WRL 2018 + + Arguments: + ----------- + im_qa: np.ndarray + Image containing the QA band + + Returns: + ----------- + cloud_mask : np.ndarray of booleans + A boolean array with True where the cloud are present + """ + + # convert QA bits + cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908] + cloud_mask = np.isin(im_qa, cloud_values) + + #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): + """ + 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) + + KV WRL 2018 + + Arguments: + ----------- + im: ee.Image() + Image to read from the Google Earth Engine database + plot_bool: boolean + True if plot is wanted + + Returns: + ----------- + im_pan: np.ndarray (2D) + The panchromatic band (15m) + im_ms: np.ndarray (3D) + The multispectral bands interpolated at 15m + im_cloud: np.ndarray (2D) + The cloud mask at 15m + crs_params: list + EPSG code and affine transformation parameters + """ + + im_dic = im.getInfo() + im_bands = im_dic.get('bands') + # delete dimensions key from dictionnary, otherwise the entire image is extracted + for i in range(len(im_bands)): del im_bands[i]['dimensions'] + # load panchromatic band + pan_band = [im_bands[7]] + im_pan, crs_pan = load_image(im, polygon, pan_band) + im_pan = im_pan[:,:,0] + # load the multispectral 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_30m, crs_ms = load_image(im, polygon, ms_bands) + # create cloud mask + 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 = transform.resize(im_cloud, (im_pan.shape[0], im_pan.shape[1]), + order=0, preserve_range=True, mode='constant').astype('bool_') + if plot_bool: + + plt.figure() + + plt.subplot(221) + plt.imshow(im_pan, cmap='gray') + plt.title('PANCHROMATIC') + + plt.subplot(222) + plt.imshow(im_ms_30m[:,:,[2,1,0]]) + plt.title('RGB') + + plt.subplot(223) + plt.imshow(im_ms_30m[:,:,3], cmap='gray') + plt.title('NIR') + + plt.subplot(224) + plt.imshow(im_ms_30m[:,:,4], cmap='gray') + plt.title('SWIR') + + plt.show() + + # resize the image using bilinear interpolation (order 1) + im_ms = transform.resize(im_ms_30m,(im_pan.shape[0], im_pan.shape[1]), + order=1, preserve_range=True, mode='constant') + + # get the crs parameters for the image at 15m and 30m resolution + crs = {'crs_15m':crs_pan, 'crs_30m':crs_ms, 'epsg_code':pan_band[0]['crs']} + + return im_pan, im_ms, im_cloud, crs + + +def rescale_image_intensity(im, cloud_mask, prob_high, plot_bool): + """ + Rescales the intensity of an image (multispectral or single band) by applying + a cloud mask and clipping the prob_high upper percentile. This functions allows + to stretch the contrast of an image. + + KV WRL 2018 + + Arguments: + ----------- + im: np.ndarray + Image to rescale, can be 3D (multispectral) or 2D (single band) + cloud_mask: np.ndarray + 2D cloud mask with True where cloud pixels are + prob_high: float + probability of exceedence used to calculate the upper percentile + plot_bool: boolean + True if plot is wanted + + Returns: + ----------- + im_adj: np.ndarray + The rescaled image + """ + prc_low = 0 # lower percentile + vec_mask = cloud_mask.reshape(im.shape[0] * im.shape[1]) + + if plot_bool: + plt.figure() + + if len(im.shape) > 2: + vec = im.reshape(im.shape[0] * im.shape[1], im.shape[2]) + vec_adj = np.ones((len(vec_mask), im.shape[2])) * np.nan + + for i in range(im.shape[2]): + prc_high = np.percentile(vec[~vec_mask, i], prob_high) + vec_rescaled = exposure.rescale_intensity(vec[~vec_mask, i], in_range=(prc_low, prc_high)) + vec_adj[~vec_mask,i] = vec_rescaled + + if plot_bool: + plt.subplot(np.floor(im.shape[2]/2) + 1, np.floor(im.shape[2]/2), i+1) + plt.hist(vec[~vec_mask, i], bins=200, label='original') + plt.hist(vec_rescaled, bins=200, alpha=0.5, label='rescaled') + plt.legend() + plt.title('Band' + str(i+1)) + plt.show() + + im_adj = vec_adj.reshape(im.shape[0], im.shape[1], im.shape[2]) + + if plot_bool: + plt.figure() + ax1 = plt.subplot(121) + plt.imshow(im[:,:,[2,1,0]]) + plt.axis('off') + plt.title('Original') + ax2 = plt.subplot(122, sharex=ax1, sharey=ax1) + plt.imshow(im_adj[:,:,[2,1,0]]) + plt.axis('off') + plt.title('Rescaled') + plt.show() + + else: + vec = im.reshape(im.shape[0] * im.shape[1]) + vec_adj = np.ones(len(vec_mask)) * np.nan + prc_high = np.percentile(vec[~vec_mask], prob_high) + vec_rescaled = exposure.rescale_intensity(vec[~vec_mask], in_range=(prc_low, prc_high)) + vec_adj[~vec_mask] = vec_rescaled + + if plot_bool: + plt.hist(vec[~vec_mask], bins=200, label='original') + plt.hist(vec_rescaled, bins=200, alpha=0.5, label='rescaled') + plt.legend() + plt.title('Single band') + plt.show() + + im_adj = vec_adj.reshape(im.shape[0], im.shape[1]) + + if plot_bool: + plt.figure() + ax1 = plt.subplot(121) + plt.imshow(im, cmap='gray') + plt.axis('off') + plt.title('Original') + ax2 = plt.subplot(122, sharex=ax1, sharey=ax1) + plt.imshow(im_adj, cmap='gray') + plt.axis('off') + plt.title('Rescaled') + plt.show() + + return im_adj + + +def hist_match(source, template): + """ + Adjust the pixel values of a grayscale image such that its histogram + matches that of a target image + + Arguments: + ----------- + source: np.ndarray + Image to transform; the histogram is computed over the flattened + array + template: np.ndarray + Template image; can have different dimensions to source + Returns: + ----------- + matched: np.ndarray + The transformed output image + """ + + oldshape = source.shape + source = source.ravel() + template = template.ravel() + + # get the set of unique pixel values and their corresponding indices and + # counts + s_values, bin_idx, s_counts = np.unique(source, return_inverse=True, + return_counts=True) + t_values, t_counts = np.unique(template, return_counts=True) + + # take the cumsum of the counts and normalize by the number of pixels to + # get the empirical cumulative distribution functions for the source and + # template images (maps pixel value --> quantile) + s_quantiles = np.cumsum(s_counts).astype(np.float64) + s_quantiles /= s_quantiles[-1] + t_quantiles = np.cumsum(t_counts).astype(np.float64) + t_quantiles /= t_quantiles[-1] + + # interpolate linearly to find the pixel values in the template image + # that correspond most closely to the quantiles in the source image + interp_t_values = np.interp(s_quantiles, t_quantiles, t_values) + + return interp_t_values[bin_idx].reshape(oldshape) + +def pansharpen(im_ms, im_pan, cloud_mask, plot_bool): + """ + Pansharpens a multispectral image (3D), using the panchromatic band (2D) + and a cloud mask + + KV WRL 2018 + + Arguments: + ----------- + im_ms: np.ndarray + Multispectral image to pansharpen (3D) + im_pan: np.ndarray + Panchromatic band (2D) + cloud_mask: np.ndarray + 2D cloud mask with True where cloud pixels are + plot_bool: boolean + True if plot is wanted + + Returns: + ----------- + im_ms_ps: np.ndarray + Pansharpened multisoectral image (3D) + """ + + # reshape image into vector and apply cloud mask + vec = im_ms.reshape(im_ms.shape[0] * im_ms.shape[1], im_ms.shape[2]) + vec_mask = cloud_mask.reshape(im_ms.shape[0] * im_ms.shape[1]) + vec = vec[~vec_mask, :] + # apply PCA to RGB bands + pca = decomposition.PCA() + vec_pcs = pca.fit_transform(vec) + # replace 1st PC with pan band (after matching histograms) + vec_pan = im_pan.reshape(im_pan.shape[0] * im_pan.shape[1]) + vec_pan = vec_pan[~vec_mask] + vec_pcs[:,0] = hist_match(vec_pan, 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])) + # reshape vector into image + vec_ms_ps_full = np.ones((len(vec_mask), im_ms.shape[2])) * np.nan + vec_ms_ps_full[~vec_mask,:] = vec_ms_ps + im_ms_ps = vec_ms_ps_full.reshape(im_ms.shape[0], im_ms.shape[1], im_ms.shape[2]) + + if plot_bool: + plt.figure() + ax1 = plt.subplot(121) + plt.imshow(im_ms[:,:,[2,1,0]]) + plt.axis('off') + plt.title('Original') + ax2 = plt.subplot(122, sharex=ax1, sharey=ax1) + plt.imshow(im_ms_ps[:,:,[2,1,0]]) + plt.axis('off') + plt.title('Pansharpened') + plt.show() + + return im_ms_ps + +def nd_index(im1, im2, cloud_mask, plot_bool): + """ + Computes normalised difference index on 2 images (2D), given a cloud mask (2D) + + KV WRL 2018 + + Arguments: + ----------- + im1, im2: np.ndarray + Images (2D) with which to calculate the ND index + cloud_mask: np.ndarray + 2D cloud mask with True where cloud pixels are + plot_bool: boolean + True if plot is wanted + + Returns: ----------- + im_nd: np.ndarray + + Image (2D) containing the ND index + """ + vec_mask = cloud_mask.reshape(im1.shape[0] * im1.shape[1]) + vec_nd = np.ones(len(vec_mask)) * np.nan + vec1 = im1.reshape(im1.shape[0] * im1.shape[1]) + vec2 = im2.reshape(im2.shape[0] * im2.shape[1]) + temp = np.divide(vec1[~vec_mask] - vec2[~vec_mask], + vec1[~vec_mask] + vec2[~vec_mask]) + vec_nd[~vec_mask] = temp + im_nd = vec_nd.reshape(im1.shape[0], im1.shape[1]) + + if plot_bool: + plt.figure() + plt.imshow(im_nd, cmap='seismic') + plt.colorbar() + plt.title('Normalised index') + plt.show() + + return im_nd + +def find_wl_contours(im_ndwi, cloud_mask, min_contour_points, plot_bool): + """ + Computes normalised difference index on 2 images (2D), given a cloud mask (2D) + + KV WRL 2018 + + Arguments: + ----------- + im_ndwi: np.ndarray + Image (2D) with the NDWI (water index) + cloud_mask: np.ndarray + 2D cloud mask with True where cloud pixels are + min_contour_points: int + minimum number of points in each contour line + plot_bool: boolean + True if plot is wanted + + Returns: ----------- + contours_wl: list of np.arrays + contains the (row,column) coordinates of the contour lines + + """ + + # reshape image to vector + vec_ndwi = im_ndwi.reshape(im_ndwi.shape[0] * im_ndwi.shape[1]) + vec_mask = cloud_mask.reshape(cloud_mask.shape[0] * cloud_mask.shape[1]) + vec = vec_ndwi[~vec_mask] + # apply otsu's threshold + t_otsu = filters.threshold_otsu(vec) + # use Marching Squares algorithm to detect contours on ndwi image + contours = measure.find_contours(im_ndwi, t_otsu) + # filter water lines + contours_wl = [] + for i, contour in enumerate(contours): + # remove contour points that are around clouds (nan values) + if np.any(np.isnan(contour)): + index_nan = np.where(np.isnan(contour))[0] + contour = np.delete(contour, index_nan, axis=0) + # remove contours that have only few points (less than min_contour_points) + if contour.shape[0] > min_contour_points: + contours_wl.append(contour) + + if plot_bool: + # plot otsu's histogram segmentation + plt.figure() + vals = plt.hist(vec, bins=200) + plt.plot([t_otsu, t_otsu],[0, np.max(vals[0])], 'r-', label='Otsu threshold') + plt.legend() + plt.show() + + # plot the water line contours on top of water index + plt.figure() + plt.imshow(im_ndwi, cmap='seismic') + plt.colorbar() + for i,contour in enumerate(contours_wl): plt.plot(contour[:, 1], contour[:, 0], linewidth=3, color='k') + plt.axis('image') + plt.title('Detected water lines') + plt.show() + + return contours_wl + +def convert_pix2world(points, crs_vec): + """ + Converts pixel coordinates (row,columns) to world projected coordinates + performing an affine transformation + + KV WRL 2018 + + Arguments: + ----------- + points: np.ndarray or list of np.ndarray + array with 2 columns (rows first and columns second) + crs_vec: np.ndarray + vector of 6 elements [Xtr, Xscale, Xshear, Ytr, Yshear, Yscale] + + Returns: ----------- + points_converted: np.ndarray or list of np.ndarray + converted coordinates, first columns with X and second column with Y + + """ + + # make affine transformation matrix + aff_mat = np.array([[crs_vec[1], crs_vec[2], crs_vec[0]], + [crs_vec[4], crs_vec[5], crs_vec[3]], + [0, 0, 1]]) + # create affine transformation + tform = transform.AffineTransform(aff_mat) + + if type(points) is list: + points_converted = [] + # iterate over the list + for i, arr in enumerate(points): + tmp = arr[:,[1,0]] + points_converted.append(tform(tmp)) + + elif type(points) is np.ndarray: + tmp = points[:,[1,0]] + points_converted = tform(tmp) + + else: + print('invalid input type') + raise + + return points_converted diff --git a/sds_extract.py b/sds_extract.py new file mode 100644 index 0000000..95c1bfc --- /dev/null +++ b/sds_extract.py @@ -0,0 +1,200 @@ +# -*- 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() + + + + diff --git a/sds_extract_collection.py b/sds_extract_collection.py new file mode 100644 index 0000000..e89292e --- /dev/null +++ b/sds_extract_collection.py @@ -0,0 +1,76 @@ +# -*- 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 +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 +import skimage.measure as measure + +# my modules +from utils import * +from sds import * + +np.seterr(all='ignore') # raise/ignore divisions by 0 and nans +ee.Initialize() + +# 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 + +# select collection +input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_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]]]; +# Dates +start_date = '2016-01-01' +end_date = '2016-12-01' +# filter by location +flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra)).filterDate(start_date, end_date) + +n_img = flt_col.size().getInfo() +print('Number of images covering Narrabeen:', n_img) +im_all = flt_col.getInfo().get('features') +output = [] +# 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 = read_eeimage(im, rect_narra, plot_bool) + # rescale intensities + im_ms = rescale_image_intensity(im_ms, im_cloud, prob_high, plot_bool) + im_pan = rescale_image_intensity(im_pan, im_cloud, prob_high, plot_bool) + # pansharpen rgb image + im_ms_ps = 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 = nd_index(im_ms_ps[:,:,3], im_ms_ps[:,:,1], im_cloud, plot_bool) + # edge detection + wl_pix = find_wl_contours(im_ndwi, im_cloud, min_contour_points, True) + # convert from pixels to world coordinates + wl_coords = convert_pix2world(wl_pix, crs['crs_15m']) + output.append(wl_coords) + diff --git a/shoreline_extraction.py b/shoreline_extraction.py new file mode 100644 index 0000000..9ef6659 --- /dev/null +++ b/shoreline_extraction.py @@ -0,0 +1,359 @@ +# -*- 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 diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..c809a18 --- /dev/null +++ b/utils.py @@ -0,0 +1,54 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Mar 1 11:30:31 2018 + +@author: z5030440 + +Contains all the utilities, convenience functions and small functions that do simple things +""" + +import matplotlib.pyplot as plt +import numpy as np + + + +def ecdf(x): + """convenience function for computing the empirical CDF""" + vals, counts = np.unique(x, return_counts=True) + ecdf = np.cumsum(counts).astype(np.float64) + ecdf /= ecdf[-1] + return vals, ecdf + +def intensity_histogram(image): + """plots histogram and cumulative distribution of the pixel intensities in an image""" + imSize = image.shape + if len(imSize) == 2: + im = image[:,:].reshape(imSize[0] * imSize[1]) + im = im[~np.isnan(im)] + fig, (ax1, ax2) = plt.subplots(2,1, sharex=True, figsize = (8,6)) + ax1.hist(im, bins=300) + ax1.set_title('Probability density function') + ax2.hist(im, bins=300, cumulative=True, histtype='step') + ax2.set_title('Cumulative distribution') + plt.show() + + else: + for i in range(imSize[2]): + im = image[:,:,i].reshape(imSize[0] * imSize[1]) + im = im[~np.isnan(im)] + fig, (ax1, ax2) = plt.subplots(2,1, sharex=True, figsize = (8,6)) + ax1.hist(im, bins=300) + ax1.set_title('Probability density function') + ax2.hist(im, bins=300, cumulative=True, histtype='step') + ax2.set_title('Cumulative distribution') + plt.show() + +def compare_images(im1, im2): + """plots 2 images next to each other, sharing the axis""" + plt.figure() + ax1 = plt.subplot(121) + plt.imshow(im1, cmap='gray') + ax2 = plt.subplot(122, sharex=ax1, sharey=ax1) + plt.imshow(im2, cmap='gray') + plt.show() +