diff --git a/README.md b/README.md index 012fed7..a5ccfb4 100644 --- a/README.md +++ b/README.md @@ -1 +1,6 @@ -This python-based toolbox allows to work with the satellite imagery provided by Google Earth Engine. \ No newline at end of file +This toolbox uses the Google Earth Engine Python API to explore collections of publicly available satellite imagery (Landsat, Sentinel). + +It has .py routines to: + - read and preprocess satellite images (cloud masking, contrast stretching) + - pansharpen Landsat 8 images + - extract shorelines with the Marching Squares algorithm \ No newline at end of file diff --git a/data/quadbike_dates.mat b/data/quadbike_dates.mat new file mode 100644 index 0000000..a59ca57 Binary files /dev/null and b/data/quadbike_dates.mat differ diff --git a/data_2016.pkl b/data_2016.pkl new file mode 100644 index 0000000..a8fbc5c Binary files /dev/null and b/data_2016.pkl differ diff --git a/sds.py b/functions/sds.py similarity index 86% rename from sds.py rename to functions/sds.py index 799e409..a4cea6c 100644 --- a/sds.py +++ b/functions/sds.py @@ -14,7 +14,7 @@ import pdb import ee # other modules -from osgeo import gdal +from osgeo import gdal, ogr, osr import tempfile from urllib.request import urlretrieve import zipfile @@ -28,7 +28,7 @@ import skimage.measure as measure # import own modules -from utils import * +from functions.utils import * # Download from ee server function @@ -130,18 +130,32 @@ def read_eeimage(im, polygon, plot_bool): crs_params: list EPSG code and affine transformation parameters """ - + im_dic = im.getInfo() + # save metadata + im_meta = im_dic.get('properties') + meta = {'timestamp':im_meta['system:time_start'], + 'date_acquired':im_meta['DATE_ACQUIRED'], + 'geom_rmse_model':im_meta['GEOMETRIC_RMSE_MODEL'], + 'gcp_model':im_meta['GROUND_CONTROL_POINTS_MODEL'], + 'quality':im_meta['IMAGE_QUALITY_OLI'], + 'sun_azimuth':im_meta['SUN_AZIMUTH'], + 'sun_elevation':im_meta['SUN_ELEVATION']} + 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) @@ -149,18 +163,38 @@ def read_eeimage(im, polygon, plot_bool): 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_') + + # 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') + + # check if -inf values (means out of image) and add to cloud mask + im_inf = np.isin(im_ms[:,:,0], -np.inf) + im_nan = np.isnan(im_ms[:,:,0]) + im_cloud = np.logical_or(np.logical_or(im_cloud, im_inf), im_nan) + + # get the crs parameters for the image at 15m and 30m resolution + crs = {'crs_15m':crs_pan, 'crs_30m':crs_ms, 'epsg_code':int(pan_band[0]['crs'][5:])} + if plot_bool: + + # if there are -inf in the image, set them to 0 before plotting + if sum(sum(np.isin(im_ms_30m[:,:,0], -np.inf).astype(int))) > 0: + idx = np.isin(im_ms_30m[:,:,0], -np.inf) + im_ms_30m[idx,0] = 0; im_ms_30m[idx,1] = 0; im_ms_30m[idx,2] = 0; + im_ms_30m[idx,3] = 0; im_ms_30m[idx,4] = 0 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') @@ -171,15 +205,8 @@ def read_eeimage(im, polygon, plot_bool): 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 - + return im_pan, im_ms, im_cloud, crs, meta + def rescale_image_intensity(im, cloud_mask, prob_high, plot_bool): """ @@ -214,7 +241,7 @@ def rescale_image_intensity(im, cloud_mask, prob_high, plot_bool): 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)) @@ -517,3 +544,47 @@ def convert_pix2world(points, crs_vec): raise return points_converted + +def convert_epsg(points, epsg_in, epsg_out): + """ + Converts from one spatial reference to another using the epsg codes + + KV WRL 2018 + + Arguments: + ----------- + points: np.ndarray or list of np.ndarray + array with 2 columns (rows first and columns second) + epsg_in: int + epsg code of the spatial reference in which the input is + epsg_out: int + epsg code of the spatial reference in which the output will be + + Returns: ----------- + points_converted: np.ndarray or list of np.ndarray + converted coordinates + + """ + + # define input and output spatial references + inSpatialRef = osr.SpatialReference() + inSpatialRef.ImportFromEPSG(epsg_in) + outSpatialRef = osr.SpatialReference() + outSpatialRef.ImportFromEPSG(epsg_out) + # create a coordinates transform + coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef) + # transform points + if type(points) is list: + points_converted = [] + # iterate over the list + for i, arr in enumerate(points): + points_converted.append(np.array(coordTransform.TransformPoints(arr))) + + elif type(points) is np.ndarray: + points_converted = np.array(coordTransform.TransformPoints(points)) + + else: + print('invalid input type') + raise + + return points_converted diff --git a/utils.py b/functions/utils.py similarity index 97% rename from utils.py rename to functions/utils.py index c809a18..fea6cd5 100644 --- a/utils.py +++ b/functions/utils.py @@ -9,7 +9,8 @@ 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 def ecdf(x): @@ -50,5 +51,4 @@ def compare_images(im1, im2): plt.imshow(im1, cmap='gray') ax2 = plt.subplot(122, sharex=ax1, sharey=ax1) plt.imshow(im2, cmap='gray') - plt.show() - + plt.show() diff --git a/sds_extract.py b/sds_extract.py index 95c1bfc..b789091 100644 --- a/sds_extract.py +++ b/sds_extract.py @@ -9,192 +9,173 @@ Main code to extract shorelines from Landsat imagery # Preamble import ee -from IPython import display -import math import matplotlib.pyplot as plt +import matplotlib.cm as cm import numpy as np +import pandas as pd +from datetime import datetime +import pickle import pdb +import pytz + + # 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 * +import skimage.measure as measure + +# my functions +import functions.utils as utils +import functions.sds as sds np.seterr(all='ignore') # raise/ignore divisions by 0 and nans ee.Initialize() -plot_bool = True +#%% 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 + +# select collection input_col = ee.ImageCollection('LANDSAT/LC08/C01/T1_RT_TOA') -# filter collection on location (Narrabeen-Collaroy rectangle) + +# 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]]]; -flt_col = input_col.filterBounds(ee.Geometry.Polygon(rect_narra)) +#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) + n_img = flt_col.size().getInfo() print('Number of images covering Narrabeen:', 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 100% cloud + if sum(sum(im_cloud.astype(int)))/(im_cloud.shape[0]*im_cloud.shape[1]) > cloud_threshold: + skipped_images[i] = True + continue + # 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']) + # 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) + # 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) + + output_wl.append(wl) + 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_pan, im_ms, im_cloud, crs, meta = sds.read_eeimage(im, rect_narra, 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_ms_ps = sds.pansharpen(im_ms[:,:,[0,1,2]], im_pan, im_cloud, plot_bool) -# 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.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.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.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.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() - +plt.show() + + diff --git a/shoreline_extraction.py b/sds_extract_allcode.py similarity index 100% rename from shoreline_extraction.py rename to sds_extract_allcode.py diff --git a/sds_extract_collection.py b/sds_extract_collection.py deleted file mode 100644 index e89292e..0000000 --- a/sds_extract_collection.py +++ /dev/null @@ -1,76 +0,0 @@ -# -*- 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) -