#6 add script to test extract shorelines

test_extract_shorelines.py

@@ -0,0 +1,166 @@
+# ==========================================================#
+# Shoreline extraction from satellite images
+# ==========================================================#
+
+# Kilian Vos WRL 2018
+
+# %% 1. Initial settings
+
+# load modules
+import os
+import numpy as np
+import pickle
+import warnings
+
+warnings.filterwarnings("ignore")
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
+
+plt.ion()
+import pandas as pd
+from scipy import interpolate
+from scipy import stats
+from datetime import datetime, timedelta
+import pytz
+from coastsat import (
+    SDS_download,
+    SDS_preprocess,
+    SDS_shoreline,
+    SDS_tools,
+    SDS_transects,
+)
+
+# region of interest (longitude, latitude in WGS84)
+polygon = [
+    [
+        [151.301454, -33.700754],
+        [151.311453, -33.702075],
+        [151.307237, -33.739761],
+        [151.294220, -33.736329],
+        [151.301454, -33.700754],
+    ]
+]
+# can also be loaded from a .kml polygon
+# kml_polygon = os.path.join(os.getcwd(), 'examples', 'NARRA_polygon.kml')
+# polygon = SDS_tools.polygon_from_kml(kml_polygon)
+# convert polygon to a smallest rectangle (sides parallel to coordinate axes)
+polygon = SDS_tools.smallest_rectangle(polygon)
+
+# date range
+# dates = ["1984-01-01", "2022-01-01"]
+# dates = ["1984-01-01", "1989-05-01"]
+dates = ["2022-01-01", "2022-12-01"]
+
+# satellite missions
+# sat_list = ["L5", "L7", "L8"]
+# sat_list = ["L5", "L7", "L8", "L9", "S2"]
+sat_list = ["L8", "L9", "S2"]
+collection = "C02"  # choose Landsat collection 'C01' or 'C02'
+# name of the site
+sitename = "extract_shorelines_test"
+
+# filepath where data will be stored
+filepath_data = os.path.join(os.getcwd(), "data")
+
+# put all the inputs into a dictionnary
+inputs = {
+    "polygon": polygon,
+    "dates": dates,
+    "sat_list": sat_list,
+    "sitename": sitename,
+    "filepath": filepath_data,
+    "landsat_collection": collection,
+}
+
+# before downloading the images, check how many images are available for your inputs
+# SDS_download.check_images_available(inputs)
+
+# %% 2. Retrieve images
+
+# only uncomment this line if you want Landsat Tier 2 images (not suitable for time-series analysis)
+inputs["include_T2"] = True
+
+# retrieve satellite images from GEE
+# metadata = SDS_download.retrieve_images(inputs)
+
+# if you have already downloaded the images, just load the metadata file
+metadata = SDS_download.get_metadata(inputs)
+
+# %% 3. Batch shoreline detection
+
+# settings for the shoreline extraction
+settings = {
+    # general parameters:
+    "cloud_thresh": 0.1,  # threshold on maximum cloud cover
+    "dist_clouds": 300,  # ditance around clouds where shoreline can't be mapped
+    "output_epsg": 28356,  # epsg code of spatial reference system desired for the output
+    # quality control:
+    "check_detection": False,  # if True, shows each shoreline detection to the user for validation
+    "adjust_detection": False,  # if True, allows user to adjust the postion of each shoreline by changing the threhold
+    "save_figure": True,  # if True, saves a figure showing the mapped shoreline for each image
+    # [ONLY FOR ADVANCED USERS] shoreline detection parameters:
+    "min_beach_area": 1000,  # minimum area (in metres^2) for an object to be labelled as a beach
+    "min_length_sl": 500,  # minimum length (in metres) of shoreline perimeter to be valid
+    "cloud_mask_issue": False,  # switch this parameter to True if sand pixels are masked (in black) on many images
+    "sand_color": "default",  # 'default', 'latest', 'dark' (for grey/black sand beaches) or 'bright' (for white sand beaches)
+    "pan_off": False,  # True to switch pansharpening off for Landsat 7/8/9 imagery
+    # add the inputs defined previously
+    "inputs": inputs,
+    "create_plot": True,  # True create a matplotlib plot of the image with the datetime as the title. False save as a standard JPG
+    "remove_cloud_mask": False,  # True, Remove cloud mask from image processing
+}
+
+# [OPTIONAL] preprocess images (cloud masking, pansharpening/down-sampling)
+# SDS_preprocess.save_jpg(metadata, settings)
+
+# [OPTIONAL] create a reference shoreline (helps to identify outliers and false detections)
+settings["reference_shoreline"] = SDS_preprocess.get_reference_sl(metadata, settings)
+# set the max distance (in meters) allowed from the reference shoreline for a detected shoreline to be valid
+settings["max_dist_ref"] = 300
+
+# extract shorelines from all images (also saves output.pkl and shorelines.kml)
+output = SDS_shoreline.extract_shorelines(metadata, settings)
+
+# remove duplicates (images taken on the same date by the same satellite)
+output = SDS_tools.remove_duplicates(output)
+# remove inaccurate georeferencing (set threshold to 10 m)
+output = SDS_tools.remove_inaccurate_georef(output, 10)
+
+# for GIS applications, save output into a GEOJSON layer
+geomtype = "points"  # choose 'points' or 'lines' for the layer geometry
+gdf = SDS_tools.output_to_gdf(output, geomtype)
+gdf = gdf.set_crs(settings["output_epsg"])
+gdf.to_crs("epsg:4326", inplace=True)
+
+if gdf is None:
+    raise Exception("output does not contain any mapped shorelines")
+# gdf.crs = {"init": "epsg:" + str(settings["output_epsg"])}  # set layer projection
+# save GEOJSON layer to file
+gdf.to_file(
+    os.path.join(
+        inputs["filepath"],
+        inputs["sitename"],
+        "%s_output_%s.geojson" % (sitename, geomtype),
+    ),
+    driver="GeoJSON",
+    encoding="utf-8",
+)
+
+geomtype = "lines"  # choose 'points' or 'lines' for the layer geometry
+gdf = SDS_tools.output_to_gdf(output, geomtype)
+gdf = gdf.set_crs(settings["output_epsg"])
+gdf.to_crs("epsg:4326", inplace=True)
+
+if gdf is None:
+    raise Exception("output does not contain any mapped shorelines")
+# gdf.crs = {"init": "epsg:" + str(settings["output_epsg"])}  # set layer projection
+# save GEOJSON layer to file
+gdf.to_file(
+    os.path.join(
+        inputs["filepath"],
+        inputs["sitename"],
+        "%s_output_%s.geojson" % (sitename, geomtype),
+    ),
+    driver="GeoJSON",
+    encoding="utf-8",
+)