diff --git a/README.md b/README.md index 00eeb39e..2ecbafbd 100644 --- a/README.md +++ b/README.md @@ -7,35 +7,42 @@ CoastSat is an open-source software toolkit written in Python that enables users to obtain time-series of shoreline position at any coastline worldwide from 30+ years (and growing) of publicly available satellite imagery. -Visit the [CoastSat website](http://coastsat.wrl.unsw.edu.au/) to explore and download regional-scale datasets of satellite-derived shorelines and beach slopes. - ![Alt text](https://github.com/kvos/CoastSat/blob/master/doc/example.gif) -The underlying approach of the CoastSat toolkit is described in detail in the following publications: +:point_right: Relevant publications: + +- Shoreline detection algorithm: https://doi.org/10.1016/j.envsoft.2019.104528 (Open Access) +- Accuracy assessment and applications: https://doi.org/10.1016/j.coastaleng.2019.04.004 +- Beach slope estimation: https://doi.org/10.1029/2020GL088365 (preprint [here](https://www.essoar.org/doi/10.1002/essoar.10502903.2)) +- Satellite-derived shorelines along meso-macrotidal beaches: https://doi.org/10.1016/j.geomorph.2021.107707 +- Beach-face slope dataset for Australia: https://doi.org/10.5194/essd-14-1345-2022 + +:point_right: Other repositories and addons related to this toolbox: +- [CoastSat.slope](https://github.com/kvos/CoastSat.slope): estimates the beach-face slope from the satellite-derived shorelines obtained with CoastSat. +- [CoastSat.PlanetScope](https://github.com/ydoherty/CoastSat.PlanetScope): shoreline extraction for PlanetScope Dove imagery (near-daily since 2017 at 3m resolution). +- [InletTracker](https://github.com/VHeimhuber/InletTracker): monitoring of intermittent open/close estuary entrances. +- [CoastSat.islands](https://github.com/mcuttler/CoastSat.islands): 2D planform measurements for small reef islands. +- [CoastSeg](https://github.com/dbuscombe-usgs/CoastSeg): image segmentation, deep learning, doodler. +- [CoastSat.Maxar](https://github.com/kvos/CoastSat.Maxar): shoreline extraction on Maxar World-View images (in progress) -1. Shoreline detection algorithm: https://doi.org/10.1016/j.envsoft.2019.104528 (Open Access) -2. Accuracy assessment and applications: https://doi.org/10.1016/j.coastaleng.2019.04.004 -3. Beach slope estimation: https://doi.org/10.1029/2020GL088365 (preprint [here](https://www.essoar.org/doi/10.1002/essoar.10502903.2)) -4. Satellite-derived shorelines along meso-macrotidal beaches: https://doi.org/10.1016/j.geomorph.2021.107707 +:point_right: Visit the [CoastSat website](http://coastsat.wrl.unsw.edu.au/) to explore and download regional-scale datasets of satellite-derived shorelines and beach slopes generated with CoastSat. -Extensions to this toolbox: -1. [CoastSat.slope](https://github.com/kvos/CoastSat.slope): estimates the beach-face slope from the satellite-derived shorelines obtained with CoastSat. -2. [CoastSat.islands](https://github.com/mcuttler/CoastSat.islands): 2D planform measurements for small reef islands. -3. [CoastSat.PlanetScope](https://github.com/ydoherty/CoastSat.PlanetScope): shoreline extraction for PlanetScope Dove imagery (near-daily since 2017 at 3m resolution). -4. [InletTracker](https://github.com/VHeimhuber/InletTracker): monitoring of intermittent open/close estuary entrances. +:star: **If you like the repo put a star on it!** :star: -### Description +### Latest updates +:arrow_forward: *(2022/05/02)* +Compatibility with Landsat 9 and Landsat Collection 2 -Satellite remote sensing can provide low-cost long-term shoreline data capable of resolving the temporal scales of interest to coastal scientists and engineers at sites where no in-situ field measurements are available. CoastSat enables the non-expert user to extract shorelines from Landsat 5, Landsat 7, Landsat 8 and Sentinel-2 images. -The shoreline detection algorithm implemented in CoastSat is optimised for sandy beach coastlines. It combines a sub-pixel border segmentation and an image classification component, which refines the segmentation into four distinct categories such that the shoreline detection is specific to the sand/water interface. +### Project description -The toolbox has three main functionalities: -- assisted retrieval from Google Earth Engine of all available satellite images spanning the user-defined region of interest and time period -- automated extraction of shorelines from all the selected images using a sub-pixel resolution technique -- intersection of the 2D shorelines with user-defined shore-normal transects -- tidal correction using measured water levels and an estimate of the beach slope +Satellite remote sensing can provide low-cost long-term shoreline data capable of resolving the temporal scales of interest to coastal scientists and engineers at sites where no in-situ field measurements are available. CoastSat enables the non-expert user to extract shorelines from Landsat 5, Landsat 7, Landsat 8, Landsat 9 and Sentinel-2 images. +The shoreline detection algorithm implemented in CoastSat is optimised for sandy beach coastlines. It combines a sub-pixel border segmentation and an image classification component, which refines the segmentation into four distinct categories such that the shoreline detection is specific to the sand/water interface. -**If you like the repo put a star on it!** +The toolbox has four main functionalities: +1. assisted retrieval from Google Earth Engine of all available satellite images spanning the user-defined region of interest and time period +2. automated extraction of shorelines from all the selected images using a sub-pixel resolution technique +3. intersection of the 2D shorelines with user-defined shore-normal transects +4. tidal correction using measured water levels and an estimate of the beach slope ## 1. Installation @@ -45,13 +52,19 @@ To run the toolbox you first need to install the required Python packages in an Once you have it installed on your PC, open the Anaconda prompt (in Mac and Linux, open a terminal window) and use the `cd` command (change directory) to go the folder where you have downloaded this repository. -Create a new environment named `coastsat` with all the required packages: +Create a new environment named `coastsat` with all the required packages by entering these commands in succession: ``` -conda env create -f environment.yml -n coastsat +conda create -n coastsat python=3.8 +conda activate coastsat +conda install -c conda-forge earthengine-api=0.1.236 +conda install gdal geopandas +conda install scikit-image +conda install -c anaconda astropy +conda install spyder notebook ``` -All the required packages have now been installed in an environment called `coastsat`. Now, activate the new environment: +All the required packages have now been installed in an environment called `coastsat`. Always make sure that the environment is activated with: ``` conda activate coastsat @@ -59,7 +72,12 @@ conda activate coastsat To confirm that you have successfully activated CoastSat, your terminal command line prompt should now start with (coastsat). -**In case errors are raised:**: open the **Anaconda Navigator**, in the *Environments* tab click on *Import* and select the `environment.yml` file. For more details, the following [link](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#creating-an-environment-with-commands) shows how to create and manage an environment with Anaconda. +:warning: **In case errors are raised** :warning:: clean things up with the following command (better to have the Anaconda Prompt open as administrator) before attempting to install `coastsat` again: +``` +conda clean --all +``` + +You can also install the packages with the **Anaconda Navigator**, in the *Environments* tab. For more details, the following [link](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html#creating-an-environment-with-commands) shows how to create and manage an environment with Anaconda. ### 1.2 Activate Google Earth Engine Python API @@ -93,22 +111,24 @@ If using `example.py` on **Spyder**, make sure that the Graphics Backend is set A Jupyter Notebook combines formatted text and code. To run the code, place your cursor inside one of the code sections and click on the `run cell` button (or press `Shift` + `Enter`) and progress forward. -![run_cell](https://user-images.githubusercontent.com/7217258/60766570-c2100080-a0ee-11e9-9675-e2aeba87e4a7.png) +![image](https://user-images.githubusercontent.com/7217258/165960239-e8870f7e-0dab-416e-bbdd-089b136b7d20.png) + ### 2.1 Retrieval of the satellite images To retrieve from the GEE server the available satellite images cropped around the user-defined region of coastline for the particular time period of interest, the following variables are required: - `polygon`: the coordinates of the region of interest (longitude/latitude pairs in WGS84) - `dates`: dates over which the images will be retrieved (e.g., `dates = ['2017-12-01', '2018-01-01']`) -- `sat_list`: satellite missions to consider (e.g., `sat_list = ['L5', 'L7', 'L8', 'S2']` for Landsat 5, 7, 8 and Sentinel-2 collections) +- `sat_list`: satellite missions to consider (e.g., `sat_list = ['L5', 'L7', 'L8', 'L9', 'S2']` for Landsat 5, 7, 8, 9 and Sentinel-2 collections) - `sitename`: name of the site (this is the name of the subfolder where the images and other accompanying files will be stored) - `filepath`: filepath to the directory where the data will be stored +- :new: `landsat_collection`: whether to use Collection 1 (`C01`) or Collection 2 (`C02`). Note that after 2022/01/01, Landsat images are only available in Collection 2. Landsat 9 is therefore only available as Collection 2. So if the user has selected `C01`, images prior to 2022/01/01 will be downloaded from Collection 1, while images captured after that date will be automatically taken from `C02`. Also note that at the time of writing `C02` is not complete in Google Earth Engine and still being uploaded. The call `metadata = SDS_download.retrieve_images(inputs)` will launch the retrieval of the images and store them as .TIF files (under */filepath/sitename*). The metadata contains the exact time of acquisition (in UTC time) of each image, its projection and its geometric accuracy. If the images have already been downloaded previously and the user only wants to run the shoreline detection, the metadata can be loaded directly by running `metadata = SDS_download.get_metadata(inputs)`. The screenshot below shows an example of inputs that will retrieve all the images of Collaroy-Narrabeen (Australia) acquired by Sentinel-2 in December 2017. -![doc1](https://user-images.githubusercontent.com/7217258/56278746-20f65700-614a-11e9-8715-ba5b8f938063.PNG) +![doc1](https://user-images.githubusercontent.com/7217258/166197244-9f41de17-f387-40a6-945e-8a78b581c4b1.png) **Note:** The area of the polygon should not exceed 100 km2, so for very long beaches split it into multiple smaller polygons. @@ -240,4 +260,6 @@ A fork is a copy on which you can make your changes. 4. Castelle B., Masselink G., Scott T., Stokes C., Konstantinou A., Marieu V., Bujan S. (2021). Satellite-derived shoreline detection at a high-energy meso-macrotidal beach. *Geomorphology*. volume 383, 107707. https://doi.org/10.1016/j.geomorph.2021.107707 -5. Training dataset used for pixel-wise classification in CoastSat: https://doi.org/10.5281/zenodo.3334147 +5. Vos, K. and Deng, W. and Harley, M. D. and Turner, I. L. and Splinter, K. D. M. (2022). Beach-face slope dataset for Australia. *Earth System Science Data*. volume 14, 3, p. 1345--1357. https://doi.org/10.5194/essd-14-1345-2022 + +6. Training dataset used for pixel-wise classification in CoastSat: https://doi.org/10.5281/zenodo.3334147 diff --git a/coastsat/SDS_download.py b/coastsat/SDS_download.py index abe9bebe..cdaf04de 100644 --- a/coastsat/SDS_download.py +++ b/coastsat/SDS_download.py @@ -5,6 +5,7 @@ Author: Kilian Vos, Water Research Laboratory, University of New South Wales """ + # load basic modules import os import numpy as np @@ -119,7 +120,7 @@ def retrieve_images(inputs): im_epsg.append(int(im_meta['bands'][0]['crs'][5:])) # get geometric accuracy - if satname in ['L5','L7','L8']: + if satname in ['L5','L7','L8','L9']: if 'GEOMETRIC_RMSE_MODEL' in im_meta['properties'].keys(): acc_georef = im_meta['properties']['GEOMETRIC_RMSE_MODEL'] else: @@ -174,7 +175,7 @@ def retrieve_images(inputs): 'epsg':im_epsg[i]} # Landsat 7 and 8 download - elif satname in ['L7', 'L8']: + elif satname in ['L7', 'L8', 'L9']: if satname == 'L7': bands['pan'] = [im_bands[8]] # panchromatic band bands['ms'] = [im_bands[0], im_bands[1], im_bands[2], im_bands[3], @@ -279,9 +280,9 @@ def retrieve_images(inputs): # once all images have been downloaded, load metadata from .txt files metadata = get_metadata(inputs) - # merge overlapping images (necessary only if the polygon is at the boundary of an image) if 'S2' in metadata.keys(): + print("\n Called merge_overlapping_images\n") try: metadata = merge_overlapping_images(metadata,inputs) except: @@ -323,7 +324,7 @@ def get_metadata(inputs): # initialize metadata dict metadata = dict([]) # loop through the satellite missions - for satname in ['L5','L7','L8','S2']: + for satname in ['L5','L7','L8','L9','S2']: # if a folder has been created for the given satellite mission if satname in os.listdir(filepath): # update the metadata dict @@ -356,13 +357,16 @@ def get_metadata(inputs): pickle.dump(metadata, f) return metadata + + ################################################################################################### # AUXILIARY FUNCTIONS ################################################################################################### def check_images_available(inputs): """ - Create the structure of subfolders for each satellite mission + Scan the GEE collections to see how many images are available for each + satellite mission (L5,L7,L8,L9,S2), collection (C01,C02) and tier (T1,T2). KV WRL 2018 @@ -379,81 +383,125 @@ def check_images_available(inputs): list of images in Tier 2 (Landsat only) """ - # check if dates are in correct order dates = [datetime.strptime(_,'%Y-%m-%d') for _ in inputs['dates']] + dates_str = inputs['dates'] + polygon = inputs['polygon'] + + # check if dates are in chronological order if dates[1] <= dates[0]: - raise Exception('Verify that your dates are in the correct order') + raise Exception('Verify that your dates are in the correct chronological order') # check if EE was initialised or not try: ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA') except: ee.Initialize() - - print('Images available between %s and %s:'%(inputs['dates'][0],inputs['dates'][1]), end='\n') - # check how many images are available in Tier 1 and Sentinel Level-1C - col_names_T1 = {'L5':'LANDSAT/LT05/C01/T1_TOA', - 'L7':'LANDSAT/LE07/C01/T1_TOA', - 'L8':'LANDSAT/LC08/C01/T1_TOA', - 'S2':'COPERNICUS/S2'} - + + print('Number of images available between %s and %s:'%(dates_str[0],dates_str[1]), end='\n') + + # get images in Landsat Tier 1 as well as Sentinel Level-1C print('- In Landsat Tier 1 & Sentinel-2 Level-1C:') + col_names_T1 = {'L5':'LANDSAT/LT05/%s/T1_TOA'%inputs['landsat_collection'], + 'L7':'LANDSAT/LE07/%s/T1_TOA'%inputs['landsat_collection'], + 'L8':'LANDSAT/LC08/%s/T1_TOA'%inputs['landsat_collection'], + 'L9':'LANDSAT/LC09/C02/T1_TOA', # only C02 for Landsat 9 + 'S2':'COPERNICUS/S2'} im_dict_T1 = dict([]) sum_img = 0 for satname in inputs['sat_list']: - - # get list of images in EE collection - while True: - try: - ee_col = ee.ImageCollection(col_names_T1[satname]) - col = ee_col.filterBounds(ee.Geometry.Polygon(inputs['polygon']))\ - .filterDate(inputs['dates'][0],inputs['dates'][1]) - im_list = col.getInfo().get('features') - break - except: - continue - # remove very cloudy images (>95% cloud cover) - im_list_upt = remove_cloudy_images(im_list, satname) - sum_img = sum_img + len(im_list_upt) - print(' %s: %d images'%(satname,len(im_list_upt))) - im_dict_T1[satname] = im_list_upt - - print(' Total: %d images'%sum_img) - - # in only S2 is in sat_list, stop here + im_list = get_image_info(col_names_T1[satname],satname,polygon,dates_str) + sum_img = sum_img + len(im_list) + print(' %s: %d images'%(satname,len(im_list))) + im_dict_T1[satname] = im_list + + # if using C01 (only goes to the end of 2021), complete with C02 for L7 and L8 + if dates[1] > datetime(2022,1,1) and inputs['landsat_collection'] == 'C01': + print(' -> completing Tier 1 with C02 after %s...'%'2022-01-01') + col_names_C02 = {'L7':'LANDSAT/LE07/C02/T1_TOA', + 'L8':'LANDSAT/LC08/C02/T1_TOA'} + dates_C02 = ['2022-01-01',dates_str[1]] + for satname in inputs['sat_list']: + if satname not in ['L7','L8']: continue # only L7 and L8 + im_list = get_image_info(col_names_C02[satname],satname,polygon,dates_C02) + sum_img = sum_img + len(im_list) + print(' %s: %d images'%(satname,len(im_list))) + im_dict_T1[satname] += im_list + + print(' Total to download: %d images'%sum_img) + + # if only S2 is in sat_list, stop here as no Tier 2 for Sentinel if len(inputs['sat_list']) == 1 and inputs['sat_list'][0] == 'S2': return im_dict_T1, [] - # otherwise check how many images are available in Landsat Tier 2 - col_names_T2 = {'L5':'LANDSAT/LT05/C01/T2_TOA', - 'L7':'LANDSAT/LE07/C01/T2_TOA', - 'L8':'LANDSAT/LC08/C01/T2_TOA'} - print('- In Landsat Tier 2:', end='\n') + # if user also requires Tier 2 images, check the T2 collections as well + col_names_T2 = {'L5':'LANDSAT/LT05/%s/T2_TOA'%inputs['landsat_collection'], + 'L7':'LANDSAT/LE07/%s/T2_TOA'%inputs['landsat_collection'], + 'L8':'LANDSAT/LC08/%s/T2_TOA'%inputs['landsat_collection']} + print('- In Landsat Tier 2 (not suitable for time-series analysis):', end='\n') im_dict_T2 = dict([]) sum_img = 0 for satname in inputs['sat_list']: - if satname == 'S2': continue - # get list of images in EE collection - while True: - try: - ee_col = ee.ImageCollection(col_names_T2[satname]) - col = ee_col.filterBounds(ee.Geometry.Polygon(inputs['polygon']))\ - .filterDate(inputs['dates'][0],inputs['dates'][1]) - im_list = col.getInfo().get('features') - break - except: - continue - # remove very cloudy images (>95% cloud cover) - im_list_upt = remove_cloudy_images(im_list, satname) - sum_img = sum_img + len(im_list_upt) - print(' %s: %d images'%(satname,len(im_list_upt))) - im_dict_T2[satname] = im_list_upt - - print(' Total: %d images'%sum_img) + if satname in ['L9','S2']: continue # no Tier 2 for Sentinel-2 and Landsat 9 + im_list = get_image_info(col_names_T2[satname],satname,polygon,dates_str) + sum_img = sum_img + len(im_list) + print(' %s: %d images'%(satname,len(im_list))) + im_dict_T2[satname] = im_list + + # also complete with C02 for L7 and L8 after 2022 + if dates[1] > datetime(2022,1,1) and inputs['landsat_collection'] == 'C01': + print(' -> completing Tier 2 with C02 after %s...'%'2022-01-01') + col_names_C02 = {'L7':'LANDSAT/LE07/C02/T2_TOA', + 'L8':'LANDSAT/LC08/C02/T2_TOA'} + dates_C02 = ['2022-01-01',dates_str[1]] + for satname in inputs['sat_list']: + if satname not in ['L7','L8']: continue # only L7 and L8 + im_list = get_image_info(col_names_C02[satname],satname,polygon,dates_C02) + sum_img = sum_img + len(im_list) + print(' %s: %d images'%(satname,len(im_list))) + im_dict_T2[satname] += im_list + + print(' Total Tier 2: %d images'%sum_img) return im_dict_T1, im_dict_T2 +def get_image_info(collection,satname,polygon,dates): + """ + Reads info about EE images for the specified collection, satellite and dates + + KV WRL 2022 + + Arguments: + ----------- + collection: str + name of the collection (e.g. 'LANDSAT/LC08/C02/T1_TOA') + satname: str + name of the satellite mission + polygon: list + coordinates of the polygon in lat/lon + dates: list of str + start and end dates (e.g. '2022-01-01') + + Returns: + ----------- + im_list: list of ee.Image objects + list with the info for the images + """ + while True: + try: + # get info about images + ee_col = ee.ImageCollection(collection) + col = ee_col.filterBounds(ee.Geometry.Polygon(polygon))\ + .filterDate(dates[0],dates[1]) + im_list = col.getInfo().get('features') + break + except: + continue + # remove very cloudy images (>95% cloud cover) + im_list = remove_cloudy_images(im_list, satname) + return im_list + + def download_tif(image, polygon, bandsId, filepath): """ Downloads a .TIF image from the ee server. The image is downloaded as a @@ -555,7 +603,7 @@ def create_folder_structure(im_folder, satname): # subfolders depending on satellite mission if satname == 'L5': filepaths.append(os.path.join(im_folder, satname, '30m')) - elif satname in ['L7','L8']: + elif satname in ['L7','L8','L9']: filepaths.append(os.path.join(im_folder, satname, 'pan')) filepaths.append(os.path.join(im_folder, satname, 'ms')) elif satname in ['S2']: @@ -591,7 +639,7 @@ def remove_cloudy_images(im_list, satname, prc_cloud_cover=95): """ # remove very cloudy images from the collection (>95% cloud) - if satname in ['L5','L7','L8']: + if satname in ['L5','L7','L8','L9']: cloud_property = 'CLOUD_COVER' elif satname in ['S2']: cloud_property = 'CLOUDY_PIXEL_PERCENTAGE' @@ -604,6 +652,9 @@ def remove_cloudy_images(im_list, satname, prc_cloud_cover=95): return im_list_upt +################################################################################################### +# Sentinel-2 ONLY +################################################################################################### def filter_S2_collection(im_list): """ @@ -660,7 +711,6 @@ def filter_S2_collection(im_list): return im_list_flt - def merge_overlapping_images(metadata,inputs): """ Merge simultaneous overlapping images that cover the area of interest. @@ -717,10 +767,14 @@ def duplicates_dict(lst): "return duplicates and indices" def duplicates(lst, item): return [i for i, x in enumerate(lst) if x == item] + return dict((x, duplicates(lst, x)) for x in set(lst) if lst.count(x) > 1) - + # first pass on images that have the exact same timestamp duplicates = duplicates_dict([_.split('_')[0] for _ in filenames]) + # {"S2-2029-2020": [0,1,2,3]} + # {"duplicate_filename": [indices of duplicated files]"} + total_removed_step1 = 0 if len(duplicates) > 0: # loop through each pair of duplicates and merge them @@ -734,9 +788,16 @@ def duplicates(lst, item): os.path.join(filepath, 'S2', '20m', filenames[idx_dup[index]].replace('10m','20m')), os.path.join(filepath, 'S2', '60m', filenames[idx_dup[index]].replace('10m','60m')), os.path.join(filepath, 'S2', 'meta', filenames[idx_dup[index]].replace('_10m','').replace('.tif','.txt'))]) - # bounding polygons - polygons.append(SDS_tools.get_image_bounds(fn_im[index][0])) - im_epsg.append(metadata[sat]['epsg'][idx_dup[index]]) + try: + # bounding polygons + polygons.append(SDS_tools.get_image_bounds(fn_im[index][0])) + im_epsg.append(metadata[sat]['epsg'][idx_dup[index]]) + except AttributeError: + print("\n Error getting the TIF. Skipping this iteration of the loop") + continue + except FileNotFoundError: + print(f"\n The file {fn_im[index][0]} did not exist") + continue # check if epsg are the same, print a warning message if len(np.unique(im_epsg)) > 1: print('WARNING: there was an error as two S2 images do not have the same epsg,'+ @@ -796,19 +857,22 @@ def duplicates(lst, item): pair_first = [_[0] for _ in pairs] idx_remove_pair = [] for idx in np.unique(pair_first): - # quadruplicate if trying to merge 3 times the same image with a successive image - if sum(pair_first == idx) == 3: - # remove the last image: 3 .tif files + the .txt file - idx_last = [pairs[_] for _ in np.where(pair_first == idx)[0]][-1][-1] - fn_im = [os.path.join(filepath, 'S2', '10m', filenames[idx_last]), - os.path.join(filepath, 'S2', '20m', filenames[idx_last].replace('10m','20m')), - os.path.join(filepath, 'S2', '60m', filenames[idx_last].replace('10m','60m')), - os.path.join(filepath, 'S2', 'meta', filenames[idx_last].replace('_10m','').replace('.tif','.txt'))] - for k in range(4): - os.chmod(fn_im[k], 0o777) - os.remove(fn_im[k]) - # store the index of the pair to remove it outside the loop - idx_remove_pair.append(np.where(pair_first == idx)[0][-1]) + # calculate the number of duplicates + n_duplicates = sum(pair_first == idx) + # if more than 3 duplicates, delete the other images so that a max of 3 duplicates are handled + if n_duplicates > 2: + for i in range(2,n_duplicates): + # remove the last image: 3 .tif files + the .txt file + idx_last = [pairs[_] for _ in np.where(pair_first == idx)[0]][i][-1] + fn_im = [os.path.join(filepath, 'S2', '10m', filenames[idx_last]), + os.path.join(filepath, 'S2', '20m', filenames[idx_last].replace('10m','20m')), + os.path.join(filepath, 'S2', '60m', filenames[idx_last].replace('10m','60m')), + os.path.join(filepath, 'S2', 'meta', filenames[idx_last].replace('_10m','').replace('.tif','.txt'))] + for k in range(4): + os.chmod(fn_im[k], 0o777) + os.remove(fn_im[k]) + # store the index of the pair to remove it outside the loop + idx_remove_pair.append(np.where(pair_first == idx)[0][i]) # remove quadruplicates from list of pairs pairs = [i for j, i in enumerate(pairs) if j not in idx_remove_pair] @@ -823,11 +887,25 @@ def duplicates(lst, item): os.path.join(filepath, 'S2', '60m', filenames[pair[index]].replace('10m','60m')), os.path.join(filepath, 'S2', 'meta', filenames[pair[index]].replace('_10m','').replace('.tif','.txt'))]) # get polygon for first image - polygon0 = SDS_tools.get_image_bounds(fn_im[0][0]) - im_epsg0 = metadata[sat]['epsg'][pair[0]] + try: + polygon0 = SDS_tools.get_image_bounds(fn_im[0][0]) + im_epsg0 = metadata[sat]['epsg'][pair[0]] + except AttributeError: + print("\n Error getting the TIF. Skipping this iteration of the loop") + continue + except FileNotFoundError: + print(f"\n The file {fn_im[index][0]} did not exist") + continue # get polygon for second image - polygon1 = SDS_tools.get_image_bounds(fn_im[1][0]) - im_epsg1 = metadata[sat]['epsg'][pair[1]] + try: + polygon1 = SDS_tools.get_image_bounds(fn_im[1][0]) + im_epsg1 = metadata[sat]['epsg'][pair[1]] + except AttributeError: + print("\n Error getting the TIF. Skipping this iteration of the loop") + continue + except FileNotFoundError: + print(f"\n The file {fn_im[index][0]} did not exist") + continue # check if epsg are the same if not im_epsg0 == im_epsg1: print('WARNING: there was an error as two S2 images do not have the same epsg,'+ diff --git a/coastsat/SDS_preprocess.py b/coastsat/SDS_preprocess.py index 142dc1a9..32fa0d9a 100644 --- a/coastsat/SDS_preprocess.py +++ b/coastsat/SDS_preprocess.py @@ -30,12 +30,12 @@ np.seterr(all='ignore') # raise/ignore divisions by 0 and nans -# Main function to preprocess a satellite image (L5,L7,L8 or S2) +# Main function to preprocess a satellite image (L5, L7, L8, L9 or S2) def preprocess_single(fn, satname, cloud_mask_issue): """ Reads the image and outputs the pansharpened/down-sampled multispectral bands, the georeferencing vector of the image (coordinates of the upper left pixel), - the cloud mask, the QA band and a no_data image. + the cloud mask, the QA band and a no_data image. For Landsat 7-8 it also outputs the panchromatic band and for Sentinel-2 it also outputs the 20m SWIR band. @@ -44,7 +44,7 @@ def preprocess_single(fn, satname, cloud_mask_issue): Arguments: ----------- fn: str or list of str - filename of the .TIF file containing the image. For L7, L8 and S2 this + filename of the .TIF file containing the image. For L7, L8 and S2 this is a list of filenames, one filename for each band at different resolution (30m and 15m for Landsat 7-8, 10m, 20m, 60m for Sentinel-2) satname: str @@ -104,8 +104,8 @@ def preprocess_single(fn, satname, cloud_mask_issue): georef[5] = -15 georef[0] = georef[0] + 7.5 georef[3] = georef[3] - 7.5 - - # check if -inf or nan values on any band and eventually add those pixels to cloud mask + + # check if -inf or nan values on any band and eventually add those pixels to cloud mask im_nodata = np.zeros(cloud_mask.shape).astype(bool) for k in range(im_ms.shape[2]): im_inf = np.isin(im_ms[:,:,k], -np.inf) @@ -117,10 +117,10 @@ def preprocess_single(fn, satname, cloud_mask_issue): for k in [1,3,4]: # loop through the Green, NIR and SWIR bands im_zeros = np.logical_and(np.isin(im_ms[:,:,k],0), im_zeros) # add zeros to im nodata - im_nodata = np.logical_or(im_zeros, im_nodata) + im_nodata = np.logical_or(im_zeros, im_nodata) # update cloud mask with all the nodata pixels cloud_mask = np.logical_or(cloud_mask, im_nodata) - + # no extra image for Landsat 5 (they are all 30 m bands) im_extra = [] @@ -157,7 +157,7 @@ def preprocess_single(fn, satname, cloud_mask_issue): # resize the image using nearest neighbour interpolation (order 0) cloud_mask = transform.resize(cloud_mask, (nrows, ncols), order=0, preserve_range=True, mode='constant').astype('bool_') - # check if -inf or nan values on any band and eventually add those pixels to cloud mask + # check if -inf or nan values on any band and eventually add those pixels to cloud mask im_nodata = np.zeros(cloud_mask.shape).astype(bool) for k in range(im_ms.shape[2]): im_inf = np.isin(im_ms[:,:,k], -np.inf) @@ -169,7 +169,7 @@ def preprocess_single(fn, satname, cloud_mask_issue): for k in [1,3,4]: # loop through the Green, NIR and SWIR bands im_zeros = np.logical_and(np.isin(im_ms[:,:,k],0), im_zeros) # add zeros to im nodata - im_nodata = np.logical_or(im_zeros, im_nodata) + im_nodata = np.logical_or(im_zeros, im_nodata) # update cloud mask with all the nodata pixels cloud_mask = np.logical_or(cloud_mask, im_nodata) @@ -187,10 +187,10 @@ def preprocess_single(fn, satname, cloud_mask_issue): im_extra = im_pan #=============================================================================================# - # L8 images + # L8 and L9 images #=============================================================================================# - elif satname == 'L8': - + elif satname in ['L8','L9']: + # read pan image fn_pan = fn[0] data = gdal.Open(fn_pan, gdal.GA_ReadOnly) @@ -219,7 +219,7 @@ def preprocess_single(fn, satname, cloud_mask_issue): # resize the image using nearest neighbour interpolation (order 0) cloud_mask = transform.resize(cloud_mask, (nrows, ncols), order=0, preserve_range=True, mode='constant').astype('bool_') - # check if -inf or nan values on any band and eventually add those pixels to cloud mask + # check if -inf or nan values on any band and eventually add those pixels to cloud mask im_nodata = np.zeros(cloud_mask.shape).astype(bool) for k in range(im_ms.shape[2]): im_inf = np.isin(im_ms[:,:,k], -np.inf) @@ -231,10 +231,10 @@ def preprocess_single(fn, satname, cloud_mask_issue): for k in [1,3,4]: # loop through the Green, NIR and SWIR bands im_zeros = np.logical_and(np.isin(im_ms[:,:,k],0), im_zeros) # add zeros to im nodata - im_nodata = np.logical_or(im_zeros, im_nodata) + im_nodata = np.logical_or(im_zeros, im_nodata) # update cloud mask with all the nodata pixels cloud_mask = np.logical_or(cloud_mask, im_nodata) - + # pansharpen Blue, Green, Red (where there is overlapping with pan band in L8) try: im_ms_ps = pansharpen(im_ms[:,:,[0,1,2]], im_pan, cloud_mask) @@ -317,7 +317,7 @@ def preprocess_single(fn, satname, cloud_mask_issue): # dilate if image was merged as there could be issues at the edges if 'merged' in fn10: im_nodata = morphology.dilation(im_nodata,morphology.square(5)) - + # update cloud mask with all the nodata pixels cloud_mask = np.logical_or(cloud_mask, im_nodata) @@ -326,7 +326,6 @@ def preprocess_single(fn, satname, cloud_mask_issue): return im_ms, georef, cloud_mask, im_extra, im_QA, im_nodata - ################################################################################################### # AUXILIARY FUNCTIONS ################################################################################################### @@ -351,11 +350,11 @@ def create_cloud_mask(im_QA, satname, cloud_mask_issue): ----------- cloud_mask : np.array boolean array with True if a pixel is cloudy and False otherwise - + """ # convert QA bits (the bits allocated to cloud cover vary depending on the satellite mission) - if satname == 'L8': + if satname in ['L8','L9']: cloud_values = [2800, 2804, 2808, 2812, 6896, 6900, 6904, 6908] elif satname == 'L7' or satname == 'L5' or satname == 'L4': cloud_values = [752, 756, 760, 764] @@ -390,12 +389,12 @@ def hist_match(source, template): array template: np.array Template image; can have different dimensions to source - + Returns: ----------- matched: np.array The transformed output image - + """ oldshape = source.shape @@ -425,9 +424,9 @@ def hist_match(source, template): def pansharpen(im_ms, im_pan, cloud_mask): """ Pansharpens a multispectral image, using the panchromatic band and a cloud mask. - A PCA is applied to the image, then the 1st PC is replaced, after histogram + A PCA is applied to the image, then the 1st PC is replaced, after histogram matching with the panchromatic band. Note that it is essential to match the - histrograms of the 1st PC and the panchromatic band before replacing and + histrograms of the 1st PC and the panchromatic band before replacing and inverting the PCA. KV WRL 2018 @@ -445,7 +444,7 @@ def pansharpen(im_ms, im_pan, cloud_mask): ----------- im_ms_ps: np.ndarray Pansharpened multispectral image (3D) - + """ # reshape image into vector and apply cloud mask @@ -469,7 +468,6 @@ def pansharpen(im_ms, im_pan, cloud_mask): return im_ms_ps - def rescale_image_intensity(im, cloud_mask, prob_high): """ Rescales the intensity of an image (multispectral or single band) by applying @@ -530,7 +528,7 @@ def rescale_image_intensity(im, cloud_mask, prob_high): def create_jpg(im_ms, cloud_mask, date, satname, filepath): """ Saves a .jpg file with the RGB image as well as the NIR and SWIR1 grayscale images. - This functions can be modified to obtain different visualisations of the + This functions can be modified to obtain different visualisations of the multispectral images. KV WRL 2018 @@ -591,7 +589,6 @@ def create_jpg(im_ms, cloud_mask, date, satname, filepath): fig.savefig(os.path.join(filepath, date + '_' + satname + '.jpg'), dpi=150) plt.close() - def save_jpg(metadata, settings, **kwargs): """ Saves a .jpg image for all the images contained in metadata. @@ -606,18 +603,18 @@ def save_jpg(metadata, settings, **kwargs): 'inputs': dict input parameters (sitename, filepath, polygon, dates, sat_list) 'cloud_thresh': float - value between 0 and 1 indicating the maximum cloud fraction in + value between 0 and 1 indicating the maximum cloud fraction in the cropped image that is accepted 'cloud_mask_issue': boolean True if there is an issue with the cloud mask and sand pixels are erroneously being masked on the images - + Returns: ----------- Stores the images as .jpg in a folder named /preprocessed - + """ - + sitename = settings['inputs']['sitename'] cloud_thresh = settings['cloud_thresh'] filepath_data = settings['inputs']['filepath'] @@ -666,7 +663,7 @@ def save_jpg(metadata, settings, **kwargs): def get_reference_sl(metadata, settings): """ Allows the user to manually digitize a reference shoreline that is used seed - the shoreline detection algorithm. The reference shoreline helps to detect + the shoreline detection algorithm. The reference shoreline helps to detect the outliers, making the shoreline detection more robust. KV WRL 2018 @@ -679,7 +676,7 @@ def get_reference_sl(metadata, settings): 'inputs': dict input parameters (sitename, filepath, polygon, dates, sat_list) 'cloud_thresh': float - value between 0 and 1 indicating the maximum cloud fraction in + value between 0 and 1 indicating the maximum cloud fraction in the cropped image that is accepted 'cloud_mask_issue': boolean True if there is an issue with the cloud mask and sand pixels @@ -690,7 +687,7 @@ def get_reference_sl(metadata, settings): Returns: ----------- reference_shoreline: np.array - coordinates of the reference shoreline that was manually digitized. + coordinates of the reference shoreline that was manually digitized. This is also saved as a .pkl and .geojson file. """ @@ -707,214 +704,207 @@ def get_reference_sl(metadata, settings): with open(os.path.join(filepath, sitename + '_reference_shoreline.pkl'), 'rb') as f: refsl = pickle.load(f) return refsl - - # otherwise get the user to manually digitise a shoreline on S2, L8 or L5 images (no L7 because of scan line error) - else: - # first try to use S2 images (10m res for manually digitizing the reference shoreline) - if 'S2' in metadata.keys(): - satname = 'S2' - filepath = SDS_tools.get_filepath(settings['inputs'],satname) - filenames = metadata[satname]['filenames'] - # if no S2 images, try L8 (15m res in the RGB with pansharpening) - elif not 'S2' in metadata.keys() and 'L8' in metadata.keys(): - satname = 'L8' - filepath = SDS_tools.get_filepath(settings['inputs'],satname) - filenames = metadata[satname]['filenames'] - # if no S2 images and no L8, use L5 images (L7 images have black diagonal bands making it - # hard to manually digitize a shoreline) - elif not 'S2' in metadata.keys() and not 'L8' in metadata.keys() and 'L5' in metadata.keys(): - satname = 'L5' - filepath = SDS_tools.get_filepath(settings['inputs'],satname) - filenames = metadata[satname]['filenames'] - else: - raise Exception('You cannot digitize the shoreline on L7 images (because of gaps in the images), add another L8, S2 or L5 to your dataset.') - - # create figure - fig, ax = plt.subplots(1,1, figsize=[18,9], tight_layout=True) - mng = plt.get_current_fig_manager() - mng.window.showMaximized() - # loop trhough the images - for i in range(len(filenames)): - # read image - fn = SDS_tools.get_filenames(filenames[i],filepath, satname) - im_ms, georef, cloud_mask, im_extra, im_QA, im_nodata = preprocess_single(fn, satname, settings['cloud_mask_issue']) + # otherwise get the user to manually digitise a shoreline on + # S2, L8, L9 or L5 images (no L7 because of scan line error) + # first try to use S2 images (10m res for manually digitizing the reference shoreline) + if 'S2' in metadata.keys(): satname = 'S2' + # if no S2 images, use L8 or L9 (15m res in the RGB with pansharpening) + elif 'L8' in metadata.keys(): satname = 'L8' + elif 'L9' in metadata.keys(): satname = 'L9' + # if no S2, L8 or L9 use L5 (30m res) + elif 'L9' in metadata.keys(): satname = 'L9' + # if only L7 images, ask user to download other images + else: + raise Exception('You cannot digitize the shoreline on L7 images (because of gaps in the images), add another L8, S2 or L5 to your dataset.') + filepath = SDS_tools.get_filepath(settings['inputs'],satname) + filenames = metadata[satname]['filenames'] + # create figure + fig, ax = plt.subplots(1,1, figsize=[18,9], tight_layout=True) + mng = plt.get_current_fig_manager() + mng.window.showMaximized() + # loop trhough the images + for i in range(len(filenames)): + + # read image + fn = SDS_tools.get_filenames(filenames[i],filepath, satname) + im_ms, georef, cloud_mask, im_extra, im_QA, im_nodata = preprocess_single(fn, satname, settings['cloud_mask_issue']) + + # compute cloud_cover percentage (with no data pixels) + cloud_cover_combined = np.divide(sum(sum(cloud_mask.astype(int))), + (cloud_mask.shape[0]*cloud_mask.shape[1])) + if cloud_cover_combined > 0.99: # if 99% of cloudy pixels in image skip + continue + + # remove no data pixels from the cloud mask (for example L7 bands of no data should not be accounted for) + cloud_mask_adv = np.logical_xor(cloud_mask, im_nodata) + # compute updated cloud cover percentage (without no data pixels) + cloud_cover = np.divide(sum(sum(cloud_mask_adv.astype(int))), + (sum(sum((~im_nodata).astype(int))))) + + # skip image if cloud cover is above threshold + if cloud_cover > settings['cloud_thresh']: + continue + + # rescale image intensity for display purposes + im_RGB = rescale_image_intensity(im_ms[:,:,[2,1,0]], cloud_mask, 99.9) + + # plot the image RGB on a figure + ax.axis('off') + ax.imshow(im_RGB) + + # decide if the image if good enough for digitizing the shoreline + ax.set_title('Press if image is clear enough to digitize the shoreline.\n' + + 'If the image is cloudy press to get another image', fontsize=14) + # set a key event to accept/reject the detections (see https://stackoverflow.com/a/15033071) + # this variable needs to be immuatable so we can access it after the keypress event + skip_image = False + key_event = {} + def press(event): + # store what key was pressed in the dictionary + key_event['pressed'] = event.key + # let the user press a key, right arrow to keep the image, left arrow to skip it + # to break the loop the user can press 'escape' + while True: + btn_keep = plt.text(1.1, 0.9, 'keep ⇨', size=12, ha="right", va="top", + transform=ax.transAxes, + bbox=dict(boxstyle="square", ec='k',fc='w')) + btn_skip = plt.text(-0.1, 0.9, '⇦ skip', size=12, ha="left", va="top", + transform=ax.transAxes, + bbox=dict(boxstyle="square", ec='k',fc='w')) + btn_esc = plt.text(0.5, 0, ' to quit', size=12, ha="center", va="top", + transform=ax.transAxes, + bbox=dict(boxstyle="square", ec='k',fc='w')) + plt.draw() + fig.canvas.mpl_connect('key_press_event', press) + plt.waitforbuttonpress() + # after button is pressed, remove the buttons + btn_skip.remove() + btn_keep.remove() + btn_esc.remove() + # keep/skip image according to the pressed key, 'escape' to break the loop + if key_event.get('pressed') == 'right': + skip_image = False + break + elif key_event.get('pressed') == 'left': + skip_image = True + break + elif key_event.get('pressed') == 'escape': + plt.close() + raise StopIteration('User cancelled checking shoreline detection') + else: + plt.waitforbuttonpress() - # compute cloud_cover percentage (with no data pixels) - cloud_cover_combined = np.divide(sum(sum(cloud_mask.astype(int))), - (cloud_mask.shape[0]*cloud_mask.shape[1])) - if cloud_cover_combined > 0.99: # if 99% of cloudy pixels in image skip - continue + if skip_image: + ax.clear() + continue + else: + # create two new buttons + add_button = plt.text(0, 0.9, 'add', size=16, ha="left", va="top", + transform=plt.gca().transAxes, + bbox=dict(boxstyle="square", ec='k',fc='w')) + end_button = plt.text(1, 0.9, 'end', size=16, ha="right", va="top", + transform=plt.gca().transAxes, + bbox=dict(boxstyle="square", ec='k',fc='w')) + # add multiple reference shorelines (until user clicks on button) + pts_sl = np.expand_dims(np.array([np.nan, np.nan]),axis=0) + geoms = [] + while 1: + add_button.set_visible(False) + end_button.set_visible(False) + # update title (instructions) + ax.set_title('Click points along the shoreline (enough points to capture the beach curvature).\n' + + 'Start at one end of the beach.\n' + 'When finished digitizing, click ', + fontsize=14) + plt.draw() - # remove no data pixels from the cloud mask (for example L7 bands of no data should not be accounted for) - cloud_mask_adv = np.logical_xor(cloud_mask, im_nodata) - # compute updated cloud cover percentage (without no data pixels) - cloud_cover = np.divide(sum(sum(cloud_mask_adv.astype(int))), - (sum(sum((~im_nodata).astype(int))))) + # let user click on the shoreline + pts = ginput(n=50000, timeout=-1, show_clicks=True) + pts_pix = np.array(pts) + # convert pixel coordinates to world coordinates + pts_world = SDS_tools.convert_pix2world(pts_pix[:,[1,0]], georef) + + # interpolate between points clicked by the user (1m resolution) + pts_world_interp = np.expand_dims(np.array([np.nan, np.nan]),axis=0) + for k in range(len(pts_world)-1): + pt_dist = np.linalg.norm(pts_world[k,:]-pts_world[k+1,:]) + xvals = np.arange(0,pt_dist) + yvals = np.zeros(len(xvals)) + pt_coords = np.zeros((len(xvals),2)) + pt_coords[:,0] = xvals + pt_coords[:,1] = yvals + phi = 0 + deltax = pts_world[k+1,0] - pts_world[k,0] + deltay = pts_world[k+1,1] - pts_world[k,1] + phi = np.pi/2 - np.math.atan2(deltax, deltay) + tf = transform.EuclideanTransform(rotation=phi, translation=pts_world[k,:]) + pts_world_interp = np.append(pts_world_interp,tf(pt_coords), axis=0) + pts_world_interp = np.delete(pts_world_interp,0,axis=0) + + # save as geometry (to create .geojson file later) + geoms.append(geometry.LineString(pts_world_interp)) + + # convert to pixel coordinates and plot + pts_pix_interp = SDS_tools.convert_world2pix(pts_world_interp, georef) + pts_sl = np.append(pts_sl, pts_world_interp, axis=0) + ax.plot(pts_pix_interp[:,0], pts_pix_interp[:,1], 'r--') + ax.plot(pts_pix_interp[0,0], pts_pix_interp[0,1],'ko') + ax.plot(pts_pix_interp[-1,0], pts_pix_interp[-1,1],'ko') + + # update title and buttons + add_button.set_visible(True) + end_button.set_visible(True) + ax.set_title('click on to digitize another shoreline or on to finish and save the shoreline(s)', + fontsize=14) + plt.draw() - # skip image if cloud cover is above threshold - if cloud_cover > settings['cloud_thresh']: - continue + # let the user click again ( another shoreline or ) + pt_input = ginput(n=1, timeout=-1, show_clicks=False) + pt_input = np.array(pt_input) - # rescale image intensity for display purposes - im_RGB = rescale_image_intensity(im_ms[:,:,[2,1,0]], cloud_mask, 99.9) - - # plot the image RGB on a figure - ax.axis('off') - ax.imshow(im_RGB) - - # decide if the image if good enough for digitizing the shoreline - ax.set_title('Press if image is clear enough to digitize the shoreline.\n' + - 'If the image is cloudy press to get another image', fontsize=14) - # set a key event to accept/reject the detections (see https://stackoverflow.com/a/15033071) - # this variable needs to be immuatable so we can access it after the keypress event - skip_image = False - key_event = {} - def press(event): - # store what key was pressed in the dictionary - key_event['pressed'] = event.key - # let the user press a key, right arrow to keep the image, left arrow to skip it - # to break the loop the user can press 'escape' - while True: - btn_keep = plt.text(1.1, 0.9, 'keep ⇨', size=12, ha="right", va="top", - transform=ax.transAxes, - bbox=dict(boxstyle="square", ec='k',fc='w')) - btn_skip = plt.text(-0.1, 0.9, '⇦ skip', size=12, ha="left", va="top", - transform=ax.transAxes, - bbox=dict(boxstyle="square", ec='k',fc='w')) - btn_esc = plt.text(0.5, 0, ' to quit', size=12, ha="center", va="top", - transform=ax.transAxes, - bbox=dict(boxstyle="square", ec='k',fc='w')) - plt.draw() - fig.canvas.mpl_connect('key_press_event', press) - plt.waitforbuttonpress() - # after button is pressed, remove the buttons - btn_skip.remove() - btn_keep.remove() - btn_esc.remove() - # keep/skip image according to the pressed key, 'escape' to break the loop - if key_event.get('pressed') == 'right': - skip_image = False - break - elif key_event.get('pressed') == 'left': - skip_image = True - break - elif key_event.get('pressed') == 'escape': - plt.close() - raise StopIteration('User cancelled checking shoreline detection') - else: - plt.waitforbuttonpress() - - if skip_image: - ax.clear() - continue - else: - # create two new buttons - add_button = plt.text(0, 0.9, 'add', size=16, ha="left", va="top", - transform=plt.gca().transAxes, - bbox=dict(boxstyle="square", ec='k',fc='w')) - end_button = plt.text(1, 0.9, 'end', size=16, ha="right", va="top", - transform=plt.gca().transAxes, - bbox=dict(boxstyle="square", ec='k',fc='w')) - # add multiple reference shorelines (until user clicks on button) - pts_sl = np.expand_dims(np.array([np.nan, np.nan]),axis=0) - geoms = [] - while 1: + # if user clicks on , save the points and break the loop + if pt_input[0][0] > im_ms.shape[1]/2: add_button.set_visible(False) end_button.set_visible(False) - # update title (instructions) - ax.set_title('Click points along the shoreline (enough points to capture the beach curvature).\n' + - 'Start at one end of the beach.\n' + 'When finished digitizing, click ', - fontsize=14) + plt.title('Reference shoreline saved as ' + sitename + '_reference_shoreline.pkl and ' + sitename + '_reference_shoreline.geojson') plt.draw() + ginput(n=1, timeout=3, show_clicks=False) + plt.close() + break - # let user click on the shoreline - pts = ginput(n=50000, timeout=-1, show_clicks=True) - pts_pix = np.array(pts) - # convert pixel coordinates to world coordinates - pts_world = SDS_tools.convert_pix2world(pts_pix[:,[1,0]], georef) - - # interpolate between points clicked by the user (1m resolution) - pts_world_interp = np.expand_dims(np.array([np.nan, np.nan]),axis=0) - for k in range(len(pts_world)-1): - pt_dist = np.linalg.norm(pts_world[k,:]-pts_world[k+1,:]) - xvals = np.arange(0,pt_dist) - yvals = np.zeros(len(xvals)) - pt_coords = np.zeros((len(xvals),2)) - pt_coords[:,0] = xvals - pt_coords[:,1] = yvals - phi = 0 - deltax = pts_world[k+1,0] - pts_world[k,0] - deltay = pts_world[k+1,1] - pts_world[k,1] - phi = np.pi/2 - np.math.atan2(deltax, deltay) - tf = transform.EuclideanTransform(rotation=phi, translation=pts_world[k,:]) - pts_world_interp = np.append(pts_world_interp,tf(pt_coords), axis=0) - pts_world_interp = np.delete(pts_world_interp,0,axis=0) - - # save as geometry (to create .geojson file later) - geoms.append(geometry.LineString(pts_world_interp)) - - # convert to pixel coordinates and plot - pts_pix_interp = SDS_tools.convert_world2pix(pts_world_interp, georef) - pts_sl = np.append(pts_sl, pts_world_interp, axis=0) - ax.plot(pts_pix_interp[:,0], pts_pix_interp[:,1], 'r--') - ax.plot(pts_pix_interp[0,0], pts_pix_interp[0,1],'ko') - ax.plot(pts_pix_interp[-1,0], pts_pix_interp[-1,1],'ko') - - # update title and buttons - add_button.set_visible(True) - end_button.set_visible(True) - ax.set_title('click on to digitize another shoreline or on to finish and save the shoreline(s)', - fontsize=14) - plt.draw() + pts_sl = np.delete(pts_sl,0,axis=0) + # convert world image coordinates to user-defined coordinate system + image_epsg = metadata[satname]['epsg'][i] + pts_coords = SDS_tools.convert_epsg(pts_sl, image_epsg, settings['output_epsg']) + + # save the reference shoreline as .pkl + filepath = os.path.join(filepath_data, sitename) + with open(os.path.join(filepath, sitename + '_reference_shoreline.pkl'), 'wb') as f: + pickle.dump(pts_coords, f) + + # also store as .geojson in case user wants to drag-and-drop on GIS for verification + for k,line in enumerate(geoms): + gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries(line)) + gdf.index = [k] + gdf.loc[k,'name'] = 'reference shoreline ' + str(k+1) + # store into geodataframe + if k == 0: + gdf_all = gdf + else: + gdf_all = gdf_all.append(gdf) + gdf_all.crs = {'init':'epsg:'+str(image_epsg)} + # convert from image_epsg to user-defined coordinate system + gdf_all = gdf_all.to_crs({'init': 'epsg:'+str(settings['output_epsg'])}) + # save as geojson + gdf_all.to_file(os.path.join(filepath, sitename + '_reference_shoreline.geojson'), + driver='GeoJSON', encoding='utf-8') + + print('Reference shoreline has been saved in ' + filepath) + break - # let the user click again ( another shoreline or ) - pt_input = ginput(n=1, timeout=-1, show_clicks=False) - pt_input = np.array(pt_input) - - # if user clicks on , save the points and break the loop - if pt_input[0][0] > im_ms.shape[1]/2: - add_button.set_visible(False) - end_button.set_visible(False) - plt.title('Reference shoreline saved as ' + sitename + '_reference_shoreline.pkl and ' + sitename + '_reference_shoreline.geojson') - plt.draw() - ginput(n=1, timeout=3, show_clicks=False) - plt.close() - break - - pts_sl = np.delete(pts_sl,0,axis=0) - # convert world image coordinates to user-defined coordinate system - image_epsg = metadata[satname]['epsg'][i] - pts_coords = SDS_tools.convert_epsg(pts_sl, image_epsg, settings['output_epsg']) - - # save the reference shoreline as .pkl - filepath = os.path.join(filepath_data, sitename) - with open(os.path.join(filepath, sitename + '_reference_shoreline.pkl'), 'wb') as f: - pickle.dump(pts_coords, f) - - # also store as .geojson in case user wants to drag-and-drop on GIS for verification - for k,line in enumerate(geoms): - gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries(line)) - gdf.index = [k] - gdf.loc[k,'name'] = 'reference shoreline ' + str(k+1) - # store into geodataframe - if k == 0: - gdf_all = gdf - else: - gdf_all = gdf_all.append(gdf) - gdf_all.crs = {'init':'epsg:'+str(image_epsg)} - # convert from image_epsg to user-defined coordinate system - gdf_all = gdf_all.to_crs({'init': 'epsg:'+str(settings['output_epsg'])}) - # save as geojson - gdf_all.to_file(os.path.join(filepath, sitename + '_reference_shoreline.geojson'), - driver='GeoJSON', encoding='utf-8') - - print('Reference shoreline has been saved in ' + filepath) - break - # check if a shoreline was digitised if len(pts_coords) == 0: raise Exception('No cloud free images are available to digitise the reference shoreline,'+ - 'download more images and try again') - + 'download more images and try again') + return pts_coords diff --git a/coastsat/SDS_tools.py b/coastsat/SDS_tools.py index afef2945..96144fc6 100644 --- a/coastsat/SDS_tools.py +++ b/coastsat/SDS_tools.py @@ -17,8 +17,6 @@ import skimage.transform as transform from astropy.convolution import convolve -np.seterr(all='ignore') # raise/ignore divisions by 0 and nans - ################################################################################################### # COORDINATES CONVERSION FUNCTIONS ################################################################################################### @@ -114,7 +112,6 @@ def convert_world2pix(points, georef): return points_converted - def convert_epsg(points, epsg_in, epsg_out): """ Converts from one spatial reference to another using the epsg codes @@ -268,7 +265,6 @@ def mask_raster(fn, mask): # close dataset and flush cache raster = None - ################################################################################################### # UTILITIES ################################################################################################### @@ -301,7 +297,7 @@ def get_filepath(inputs,satname): 'sat_list': list of str list that contains the names of the satellite missions to include: ``` - sat_list = ['L5', 'L7', 'L8', 'S2'] + sat_list = ['L5', 'L7', 'L8', 'L9', 'S2'] ``` 'filepath_data': str filepath to the directory where the images are downloaded @@ -321,15 +317,10 @@ def get_filepath(inputs,satname): if satname == 'L5': # access downloaded Landsat 5 images filepath = os.path.join(filepath_data, sitename, satname, '30m') - elif satname == 'L7': + elif satname in ['L7','L8','L9']: # access downloaded Landsat 7 images - filepath_pan = os.path.join(filepath_data, sitename, 'L7', 'pan') - filepath_ms = os.path.join(filepath_data, sitename, 'L7', 'ms') - filepath = [filepath_pan, filepath_ms] - elif satname == 'L8': - # access downloaded Landsat 8 images - filepath_pan = os.path.join(filepath_data, sitename, 'L8', 'pan') - filepath_ms = os.path.join(filepath_data, sitename, 'L8', 'ms') + filepath_pan = os.path.join(filepath_data, sitename, satname, 'pan') + filepath_ms = os.path.join(filepath_data, sitename, satname, 'ms') filepath = [filepath_pan, filepath_ms] elif satname == 'S2': # access downloaded Sentinel 2 images @@ -364,7 +355,7 @@ def get_filenames(filename, filepath, satname): if satname == 'L5': fn = os.path.join(filepath, filename) - if satname == 'L7' or satname == 'L8': + if satname in ['L7','L8','L9']: filename_ms = filename.replace('pan','ms') fn = [os.path.join(filepath[0], filename), os.path.join(filepath[1], filename_ms)] @@ -716,11 +707,19 @@ def GetExtent(gt,cols,rows): return ext # load .tif file and get bounds + if not os.path.exists(fn): + raise FileNotFoundError(f"{fn}") data = gdal.Open(fn, gdal.GA_ReadOnly) - gt = data.GetGeoTransform() - cols = data.RasterXSize - rows = data.RasterYSize - ext = GetExtent(gt,cols,rows) + # Check if data is null meaning the open failed + if data is None: + print("TIF file: ",fn, "cannot be opened" ) + os.remove(fn) + raise AttributeError + else: + gt = data.GetGeoTransform() + cols = data.RasterXSize + rows = data.RasterYSize + ext = GetExtent(gt,cols,rows) return geometry.Polygon(ext) diff --git a/environment.yml b/environment.yml index 0db055a6..6207efde 100644 --- a/environment.yml +++ b/environment.yml @@ -2,20 +2,12 @@ channels: - defaults - conda-forge dependencies: - - python=3.7 + - python=3.8 - earthengine-api=0.1.236 - - google-api-python-client==1.12.8 - - gdal=2.3.3 - - pandas=0.24.2 - - geopandas=0.4.1 - - numpy + - gdal + - geopandas - matplotlib - - pillow - - pytz - scikit-image - - scikit-learn - - shapely - - scipy - astropy - spyder - notebook diff --git a/example.py b/example.py index b2c7155c..e37c1cdb 100644 --- a/example.py +++ b/example.py @@ -36,7 +36,7 @@ # satellite missions sat_list = ['S2'] - +collection = 'C01' # choose Landsat collection 'C01' or 'C02' # name of the site sitename = 'NARRA' @@ -49,7 +49,8 @@ 'dates': dates, 'sat_list': sat_list, 'sitename': sitename, - 'filepath': filepath_data + 'filepath': filepath_data, + 'landsat_collection': collection } # before downloading the images, check how many images are available for your inputs diff --git a/example_jupyter.ipynb b/example_jupyter.ipynb index ea2b250c..d14025a6 100644 --- a/example_jupyter.ipynb +++ b/example_jupyter.ipynb @@ -12,7 +12,7 @@ "- Beach slope estimation: https://doi.org/10.1029/2020GL088365\n", "\n", "It enables the users to extract time-series of shoreline change over the last 30+ years at their site of interest.\n", - "There are three main steps:\n", + "There are four main steps:\n", "1. Retrieval of the satellite images of the region of interest from Google Earth Engine\n", "2. Shoreline extraction at sub-pixel resolution\n", "3. Intersection of the shorelines with cross-shore transects\n", @@ -29,6 +29,8 @@ "metadata": {}, "outputs": [], "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", "import os\n", "import numpy as np\n", "import pickle\n", @@ -54,7 +56,9 @@ "\n", "Make sure the area of your ROI is smaller than 100 km2 (if larger split it into smaller ROIs).\n", "\n", - "The function `SDS_download.check_images_available(inputs)` will print the number of images available for your inputs. The Landsat images are divided in Tier 1 and Tier 2, only Tier 1 images can be used for time-series analysis. " + "The function `SDS_download.check_images_available(inputs)` will print the number of images available for your inputs. The Landsat images are divided in Tier 1 and Tier 2, only Tier 1 images can be used for time-series analysis.\n", + "\n", + "For Landsat, users can also choose between Collection 1 and Collection 2 with the `collection` variable. Note that at the time of writing (05/02/2022), Collection 2 is still being uploaded in the Google Earth Engine servers and may not be complete. Also only the satellite-derived shorelines extracted from Collection 1 were validated in previous work. Landsat 9 is only available from Collection 2, and any Landsat image after 01/01/2022 is also only available in Collection 2." ] }, { @@ -84,14 +88,17 @@ "polygon = SDS_tools.smallest_rectangle(polygon)\n", "# date range\n", "dates = ['2017-12-01', '2018-01-01']\n", - "# satellite missions\n", + "# satellite missions ['L5','L7','L8','L9','S2']\n", "sat_list = ['S2']\n", + "# choose Landsat collection 'C01' or 'C02'\n", + "collection = 'C01'\n", "# name of the site\n", "sitename = 'NARRA'\n", "# directory where the data will be stored\n", "filepath = os.path.join(os.getcwd(), 'data')\n", "# put all the inputs into a dictionnary\n", - "inputs = {'polygon': polygon, 'dates': dates, 'sat_list': sat_list, 'sitename': sitename, 'filepath':filepath}\n", + "inputs = {'polygon': polygon, 'dates': dates, 'sat_list': sat_list, 'sitename': sitename, 'filepath':filepath,\n", + " 'landsat_collection': collection}\n", "\n", "# before downloading the images, check how many images are available for your inputs\n", "SDS_download.check_images_available(inputs);"