diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..d088f68
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+*.pyc
+*.mat
+*.tif
+*.png
+*.mp4
+*.gif
\ No newline at end of file
diff --git a/.gitingnore b/.gitingnore
deleted file mode 100644
index 7e99e36..0000000
--- a/.gitingnore
+++ /dev/null
@@ -1 +0,0 @@
-*.pyc
\ No newline at end of file
diff --git a/functions/sds.py b/functions/sds.py
index a4cea6c..ec62f0e 100644
--- a/functions/sds.py
+++ b/functions/sds.py
@@ -25,6 +25,7 @@ import skimage.exposure as exposure
 import skimage.transform as transform
 import sklearn.decomposition as decomposition
 import skimage.measure as measure
+import skimage.morphology as morphology
 
 
 # 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)]
     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
     
@@ -89,6 +90,10 @@ def create_cloud_mask(im_qa):
     -----------
         im_qa: np.ndarray
             Image containing the QA band
+        satname: string
+            short name for the satellite (L8, L7, S2)
+        plot_bool: boolean
+            True if plot is wanted
             
     Returns:
     -----------
@@ -97,15 +102,27 @@ def create_cloud_mask(im_qa):
     """
     
     # 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)
+    # 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_mask = np.isin(im_qa, cloud_shadow_values)
     
     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)
     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]]
     im_qa, crs_qa = load_image(im, polygon, qa_band)
     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]),
                                 order=0, preserve_range=True, mode='constant').astype('bool_')
     
diff --git a/functions/utils.py b/functions/utils.py
index 35aecca..019c0d4 100644
--- a/functions/utils.py
+++ b/functions/utils.py
@@ -9,7 +9,6 @@ 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
 
 
@@ -56,3 +55,7 @@ def compare_images(im1, im2):
 def find_indices(lst, condition):
     "imitation of MATLAB find function"
     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)]
diff --git a/sds_extract.py b/sds_extract.py
index e04675d..8de4f40 100644
--- a/sds_extract.py
+++ b/sds_extract.py
@@ -1,13 +1,10 @@
 # -*- coding: utf-8 -*-
-"""
-Created on Thu Mar  1 14:32:08 2018
 
-@author: z5030440
-
-Main code to extract shorelines from Landsat imagery
-"""
-# Preamble
+#==========================================================#
+# Extract shorelines from Landsat images
+#==========================================================#
 
+# Initial settings
 import ee
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
@@ -27,206 +24,78 @@ import sklearn.decomposition as decomposition
 import skimage.morphology as morphology
 import skimage.measure as measure
 
-# my functions
+# my modules
 import functions.utils as utils
 import functions.sds as sds
 
+# some settings
 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()
 
-#%% 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
+cloud_thresh = 0.5      # threshold for cloud cover
+plot_bool = True      # 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
 
 # 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)
-rect_narra = [[[151.3473129272461,-33.69035274454718],
-              [151.2820816040039,-33.68206818063878], 
-              [151.27281188964844,-33.74775138989556],
-              [151.3425064086914,-33.75231878701767],
-              [151.3473129272461,-33.69035274454718]]];
-               
-with open('data/narra_beach.pkl', 'rb') as f:
-    pts_beach = pickle.load(f)
-               
-#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)
+polygon = [[[151.3473129272461,-33.69035274454718],
+            [151.2820816040039,-33.68206818063878], 
+            [151.27281188964844,-33.74775138989556],
+            [151.3425064086914,-33.75231878701767],
+            [151.3473129272461,-33.69035274454718]]];
+
+# dates
+start_date = '2013-01-01'
+end_date = '2018-12-31'
+
+# filter by location and date
+flt_col = input_col.filterBounds(ee.Geometry.Polygon(polygon)).filterDate(start_date, end_date)
 
 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')
 
-#%% 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 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()
+i = 0 # first image
     
-    # 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
+# find image in ee database
 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_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('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.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.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)
diff --git a/sds_extract_allcode.py b/sds_extract_allcode.py
deleted file mode 100644
index 9ef6659..0000000
--- a/sds_extract_allcode.py
+++ /dev/null
@@ -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()
\ No newline at end of file