v0.1

MS algorithm for edge detection and transformation to world coordinates

README.md

@@ -1 +1,6 @@
-This python-based toolbox allows to work with the satellite imagery provided by Google Earth Engine.
+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

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
@@ -132,16 +132,30 @@ def read_eeimage(im, polygon, plot_bool):
     """   
 
     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,8 +163,27 @@ 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)
@@ -161,6 +194,7 @@ def read_eeimage(im, polygon, plot_bool):
         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,14 +205,7 @@ def read_eeimage(im, polygon, plot_bool):
         
         plt.show()
     
-    # resize the image using bilinear interpolation (order 1)
+    return im_pan, im_ms, im_cloud, crs, meta
-    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):
@@ -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

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):
@@ -51,4 +52,3 @@ def compare_images(im1, im2):
     ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
     plt.imshow(im2, cmap='gray')
     plt.show()  
-    

sds_extract.py

@@ -9,191 +9,172 @@ 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
+import skimage.measure as measure
-from utils import *
+
-from sds1 import *
+# 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.axis('image')
-plt.show()
+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.subplot(121)
+plt.axis('equal')
-plt.imshow(im_ms_adj[:,:,[3,1,0]])
+plt.grid()
-plt.title('Original NIR-GB')
+cmap = cm.get_cmap('jet')
-plt.show()
+colours = cmap(np.linspace(0, 1, num=len(idx_best)))
-
+for i, idx in enumerate(idx_best):
-plt.subplot(122)
+    for j in range(len(data['contours'][i])):
-plt.imshow(im_ms_ps[:,:,[3,1,0]])
+        if np.any(np.linalg.norm(data['contours'][i][j][:,[0,1]] - point_narra, axis=1) < 200):
-plt.title('Pansharpened NIR-GB')
+            plt.plot(data['contours'][i][j][:,0], data['contours'][i][j][:,1],
-plt.show()
+                     label=str(data['date_acquired'][i]),
-
+                     linewidth=2, color=colours[i,:])
-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()
     
 
     

sds_extract_collection.py

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