VedgeSat_DriverTemplate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e30042e3",
   "metadata": {},
   "source": [
    "# Notebook to drive Vegetation Edge extraction from Satellite Images\n",
    "The programming language we are using is called Python. The code has all been written and this notebook will guide you through modifying the analysis for your own area of interest, and executing the analysis.\n",
    "\n",
    "**To run a code block, click in a cell, hold down shift, and press enter.** An asterisk in square brackets `In [*]:` will appear while the code is being executed, and this will change to a number `In [1]:` when the code is finished. *The order in which you execute the code blocks matters, they must be run in sequence.* Some cells are optional however! Please read the descriptions to decide if you need to run each cell.\n",
    "\n",
    "If the heading above a cell is <font color='red'>IN RED</font>, this means you should edit this cell with your own requirements and settings.\n",
    "\n",
    "Inside blocks of python code there are comments indicated by lines that start with `#`. These lines are not computer code but rather comments providing information about what the code is doing to help you follow along and troubleshoot. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24582257",
   "metadata": {},
   "source": [
    "### Import Python packages/modules\n",
    "Load in the packages we'll need for running the tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad15b040",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import pickle\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import matplotlib\n",
    "matplotlib.use('Qt5Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "plt.ion()\n",
    "from datetime import datetime\n",
    "from Toolshed import Download, Toolbox, VegetationLine, Plotting, PlottingSeaborn, Transects\n",
    "import ee\n",
    "import geopandas as gpd\n",
    "import geemap\n",
    "\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc712e56",
   "metadata": {},
   "source": [
    "### <font color='red'>User requirements</font>\n",
    "Edit these variables to set up your site of interest, including desired coordinate projections, site name and bounding box, satellite platforms (`L5` is Landsat 5, `L8` is Landsat 8, and `S2` is Sentinel 2) and date range of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1384b8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define AOI using coordinates of a rectangle\n",
    "# The points represent the corners of a bounding box that go around your site\n",
    "sitename = 'SITENAME'\n",
    "\n",
    "# Date range\n",
    "dates = ['2021-05-01', '2021-07-02']\n",
    "\n",
    "# Satellite missions\n",
    "# Input a list of containing any/all of 'L5', 'L7', 'L8', 'L9', 'S2', 'PSScene4Band'\n",
    "# L5: 1984-2013; L7: 1999-2017 (SLC error from 2003); L8: 2013-present; S2: 2014-present; L9: 2021-present\n",
    "sat_list = ['L5','L8','S2']\n",
    "\n",
    "# Cloud threshold for screening out cloudy imagery (0.5 or 50% recommended)\n",
    "cloud_thresh = 0.5\n",
    "\n",
    "# Extract shoreline (wet-dry boundary) as well as veg edge\n",
    "wetdry = True\n",
    "\n",
    "# Reference shoreline/veg line shapefile name (should be stored in a folder called referenceLines in Data)\n",
    "# Line should be ONE CONTINUOUS linestring along the shore, stored as a shapefile in WGS84 coord system\n",
    "referenceLineShp = 'SITENAME_refLine.shp'\n",
    "# Maximum amount in metres by which to buffer the reference line for capturing veg edges within\n",
    "max_dist_ref = 150"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "019b80d3",
   "metadata": {},
   "source": [
    "### Set Up Site Directory\n",
    "Make a series of directories for the new sitename."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bcb8b10-e7f2-4de3-9af9-84d5274823be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory where the data will be stored\n",
    "filepath = Toolbox.CreateFileStructure(sitename, sat_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12cf39f9-3336-47af-91d2-53fd8220968b",
   "metadata": {},
   "source": [
    "### Reference Shore Option 1: Define AOI from the bounding box of the <font color='red'>reference shoreline shapefile</font>\n",
    "You shouldn't need to change anything here; the `referenceLineShp` comes from the <font color='red'>User requirements</font> cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b60ea3d8-6713-4f22-a3be-74e4058fa3a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Return AOI from reference line bounding box and save AOI folium map HTML in sitename directory\n",
    "referenceLinePath = os.path.join(filepath, 'referenceLines', referenceLineShp)\n",
    "referenceLineDF = gpd.read_file(referenceLinePath)\n",
    "polygon, point, lonmin, lonmax, latmin, latmax = Toolbox.AOIfromLine(referenceLinePath, max_dist_ref, sitename)\n",
    "\n",
    "# It's recommended to convert the polygon to the smallest rectangle (sides parallel to coordinate axes)       \n",
    "polygon = Toolbox.smallest_rectangle(polygon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e025f9a-c3a2-40c3-b4a2-b1eb197881a0",
   "metadata": {},
   "source": [
    "### Reference Shore Option 2: <font color='red'>Draw desired AOI on map</font>\n",
    "Use the polygon or rectangle tool in the Leaflet map window (that will appear below the next cell) to draw out your area of interest. Defining only one feature is valid at present. Once you are happy with your AOI, run the next cell to save the feature as your AOI for gathering satellite metadata within."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c3b3769-393f-46f5-adf2-c0451f077b11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell to generate an interactive map\n",
    "Map = geemap.Map(center=[0,0],zoom=2)\n",
    "Map.add_basemap('HYBRID')\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a9a7d39-f93d-4233-be52-6aeffb3bb0e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell after defining your region of interest\n",
    "roi = Map.user_roi.geometries().getInfo()[0]['coordinates']\n",
    "polygon = [[roi[0][0],roi[0][3],roi[0][1],roi[0][2]]]\n",
    "point = ee.Geometry.Point(roi[0][0])\n",
    "\n",
    "polygon, point, BBoxGDF = Toolbox.AOIfromLeaflet(polygon, point, sitename)\n",
    "NewBBox = BBoxGDF.to_json()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "275f2c99-f464-4523-8ec9-2c89cb174a33",
   "metadata": {},
   "outputs": [],
   "source": [
    "NewBBox"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7078c4a-4a41-42ae-b104-06b9e8f7b4d8",
   "metadata": {},
   "source": [
    "### Reference Shore Option 3: Define AOI from <font color='red'>user defined coordinates</font> \n",
    "Make sure you get the latitude and longitude around the right way, and minimum and maximum correct (including any - signs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c412314",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDIT THESE to define an AOI bounding box (minimum and maxiumum latitudes and longitudes)\n",
    "lonmin, lonmax = -2.84869, -2.79878\n",
    "latmin, latmax = 56.32641, 56.39814\n",
    "\n",
    "# Return AOI after checking coords and saving folium map HTML in sitename directory\n",
    "polygon, point = Toolbox.AOI(lonmin, lonmax, latmin, latmax, sitename)\n",
    "\n",
    "# It's recommended to convert the polygon to the smallest rectangle (sides parallel to coordinate axes)       \n",
    "polygon = Toolbox.smallest_rectangle(polygon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2789d6c",
   "metadata": {},
   "source": [
    "### Compile Input Settings for Gathering Imagery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b179f10a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if len(dates)>2:\n",
    "    daterange='no'\n",
    "else:\n",
    "    daterange='yes'\n",
    "years = list(Toolbox.daterange(datetime.strptime(dates[0],'%Y-%m-%d'), datetime.strptime(dates[-1],'%Y-%m-%d')))\n",
    "\n",
    "# Put all the inputs into a dictionary\n",
    "inputs = {'polygon': polygon, 'dates': dates, 'daterange':daterange, 'sat_list': sat_list, 'sitename': sitename, 'filepath':filepath}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56855b8c",
   "metadata": {},
   "source": [
    "### Image Retrieval\n",
    "Before downloading the images, check how many images are available for your inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2408a975",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = Download.check_images_available(inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0363f53f",
   "metadata": {},
   "source": [
    "### Image Download\n",
    "Compile the metadata from the Google Earth Engine server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "196a0e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "Sat = Toolbox.image_retrieval(inputs)\n",
    "metadata = Toolbox.metadata_collection(inputs, Sat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a3a880f",
   "metadata": {},
   "source": [
    "### Local Image Retrieval (Planet)\n",
    "For the time being, Planet API is very slow at downloading files (depending on the size of your area). You may want to use the [Planet online data](https://developers.planet.com/docs/apps/explorer/) download portal and move the files to local folders that will have been created at the Toolbox.CreateFileStructure() step. \n",
    "- Move the image TIFFs into 'Data/YOURSITENAME/local_images/PlanetScope';\n",
    "- and the respective cloud masks into 'Data/YOURSITENAME/local_images/PlanetScope/cloudmasks'.\n",
    "- You can move any leftover extra files into 'Data/YOURSITENAME/AuxiliaryImages'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a99d4f52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you want to include Landsat 7 but DON'T want to include Scan Line Corrector affected images, set SLC=False\n",
    "Sat = Download.RetrieveImages(inputs, SLC=True)\n",
    "metadata = Download.CollectMetadata(inputs, Sat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b8a14c",
   "metadata": {},
   "source": [
    "### <font color='red'>Vegetation Edge Settings</font>\n",
    "ONLY EDIT THESE IF ADJUSTMENTS ARE NEEDED; for example, if you are getting no/disjointed edges, or  if you want to adjust the threshold used for detection in each image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1765c639",
   "metadata": {},
   "outputs": [],
   "source": [
    "BasePath = 'Data/' + sitename + '/lines'\n",
    "\n",
    "if os.path.isdir(BasePath) is False:\n",
    "    os.mkdir(BasePath)\n",
    "\n",
    "projection_epsg, _ = Toolbox.FindUTM(polygon[0][0][1],polygon[0][0][0])\n",
    "\n",
    "settings = {\n",
    "    # general parameters:\n",
    "    'cloud_thresh': cloud_thresh,        # threshold on maximum cloud cover\n",
    "    'output_epsg': projection_epsg,     # epsg code of spatial reference system desired for the output   \n",
    "    'wetdry': wetdry,              # extract wet-dry boundary as well as veg\n",
    "    # quality control:\n",
    "    'check_detection': True,    # if True, shows each shoreline detection to the user for validation\n",
    "    'adjust_detection': False,  # if True, allows user to adjust the postion of each shoreline by changing the threhold\n",
    "    'save_figure': True,        # if True, saves a figure showing the mapped shoreline for each image\n",
    "    # [ONLY FOR ADVANCED USERS] shoreline detection parameters:\n",
    "    'min_beach_area': 200,     # minimum area (in metres^2) for an object to be labelled as a beach\n",
    "    'buffer_size': 250,         # radius (in metres) for buffer around sandy pixels considered in the shoreline detection\n",
    "    'min_length_sl': 500,       # minimum length (in metres) of shoreline perimeter to be valid\n",
    "    'cloud_mask_issue': False,  # switch this parameter to True if sand pixels are masked (in black) on many images  \n",
    "    # add the inputs defined previously\n",
    "    'inputs': inputs,\n",
    "    'projection_epsg': projection_epsg,\n",
    "    'year_list': years\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25e2d531",
   "metadata": {},
   "source": [
    "### <font color='red'>Tidal Information</font>\n",
    "Compute tides from FES2014 or FES2022 data which is downloaded from the pyFES server. Only relevant for shoreline processing (as it is used to correct for the effects of tides on the cross-shore waterline position). \n",
    "\n",
    "(ONLY RUN IF YOU HAVE `pyfes` INSTALLED AND WANT TIDAL INFO SAVED. If you do, change `tidepath` to the path to your `aviso-fes` folder, see the README for details)\n",
    "\n",
    "Note: FES2022 is more accurate than FES2014 but takes several minutes longer to compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25e49142",
   "metadata": {},
   "outputs": [],
   "source": [
    "if wetdry is True:\n",
    "    tidepath = \"/path/to/your/aviso-fes/data/fes2014\"\n",
    "    daterange = dates\n",
    "    tidelatlon = [(latmin+latmax)/2, (lonmin+lonmax)/2] # centre of bounding box\n",
    "    Toolbox.ComputeTides(settings,tidepath,daterange,tidelatlon) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a970a990",
   "metadata": {},
   "source": [
    "### Vegetation Edge Reference Line Load-In\n",
    "Read in a shapefile representing the rough edge of vegetation that you want to investigate along. Does not neet to be accurate as it is used to create a coastal buffer for constraining the edge extraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98f7bc10",
   "metadata": {},
   "outputs": [],
   "source": [
    "referenceLinePath = os.path.join(inputs['filepath'], 'referenceLines', referenceLineShp)\n",
    "referenceLine, ref_epsg = Toolbox.ProcessRefline(referenceLinePath,settings)\n",
    "\n",
    "settings['reference_shoreline'] = referenceLine\n",
    "settings['ref_epsg'] = ref_epsg\n",
    "# Distance to buffer reference line by (this is in metres)\n",
    "settings['max_dist_ref'] = max_dist_ref"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ccb641",
   "metadata": {},
   "source": [
    "### <font color='red'>Reference Image for Coregistration</font>\n",
    "You can now coregister your satellite images using AROSICS. If you want to try coregistering your images to improve timeseries accuracy, <font color='red'>provide a filepath to a reference RGB image.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7734d9d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "settings['reference_coreg_im'] = None # leave as None if no coregistration is to be performed"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85492e2e",
   "metadata": {},
   "source": [
    "### Vegetation Line Extraction\n",
    "**OPTION 1**: Run the extraction tool and return output veg edges as a dictionary of lines with associated image info attached.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac527f86",
   "metadata": {},
   "outputs": [],
   "source": [
    "output, output_latlon, output_proj = VegetationLine.extract_veglines(metadata, settings, polygon, dates)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc12a784",
   "metadata": {},
   "source": [
    "### Vegetation Line Extraction Load-In\n",
    "**OPTION 2**: Load in pre-existing output file generated from a previous run of the vegetation edge extraction tool.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c545df8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "output, output_latlon, output_proj = Toolbox.ReadOutput(inputs)\n",
    "\n",
    "# Remove Duplicate Lines\n",
    "# For images taken on the same date by the same satellite, keep only the longest line\n",
    "output = Toolbox.RemoveDuplicates(output) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9621c82",
   "metadata": {},
   "source": [
    "### Save Veglines as Local Shapefiles\n",
    "You **DO NOT** need to run this again if you have already exported a previous run to a shapefile, especially not if you have then edited/tidied any unwanted lines for further analysis in an external GIS software (as rerunning this cell will overwrite your edits!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "040c8bb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save output veglines \n",
    "Toolbox.SaveConvShapefiles(output, BasePath, sitename, settings['output_epsg'])\n",
    "# Save output shorelines if they were generated\n",
    "if settings['wetdry'] == True:\n",
    "    Toolbox.SaveConvShapefiles_Water(output, BasePath, sitename, settings['output_epsg'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2431a28f",
   "metadata": {},
   "source": [
    "### <font color='red'>Define Settings for Cross-shore Transects</font>\n",
    "Edit me to define the characteristics of the cross-shore transects generated for intersecting with the extracted vegetation edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87874b42",
   "metadata": {},
   "outputs": [],
   "source": [
    "SmoothingWindowSize = 21 \n",
    "NoSmooths = 100\n",
    "TransectSpacing = 10\n",
    "DistanceInland = 100\n",
    "DistanceOffshore = 100\n",
    "\n",
    "# Provide average beach slope (tanBeta) for site, for calculating corrected beach widths\n",
    "# Set to 'None' if you want to use CoastSat.slope to calculate per-transect slopes for correcting with\n",
    "beachslope = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f0cbace",
   "metadata": {},
   "source": [
    "### Create Cross-shore Transects\n",
    "Generate transects based on the above settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83f6b995",
   "metadata": {},
   "outputs": [],
   "source": [
    "VegBasePath = 'Data/' + sitename + '/lines'\n",
    "VeglineShp = glob.glob(BasePath+'/*veglines.shp')\n",
    "VeglineGDF = gpd.read_file(VeglineShp[0])\n",
    "VeglineGDF = VeglineGDF.sort_values(by='dates') # sort GDF by dates to ensure transect intersects occur in chronological order\n",
    "VeglineGDF = VeglineGDF.reset_index(drop=True) # reset GDF index after date sorting\n",
    "if settings['wetdry'] == True:\n",
    "    WaterlineShp = glob.glob(BasePath+'/*waterlines.shp')\n",
    "    WaterlineGDF = gpd.read_file(WaterlineShp[0])\n",
    "    WaterlineGDF = WaterlineGDF.sort_values(by='dates') # as above with VeglineGDF date sorting\n",
    "    WaterlineGDF = WaterlineGDF.reset_index(drop=True)\n",
    "# Produces Transects for the reference line\n",
    "TransectSpec =  os.path.join(BasePath, sitename+'_Transects.shp')\n",
    "\n",
    "# If transects already exist, load them in\n",
    "if os.path.isfile(TransectSpec[:-3]+'pkl') is False:\n",
    "    TransectGDF = Transects.ProduceTransects(settings, SmoothingWindowSize, NoSmooths, TransectSpacing, DistanceInland, DistanceOffshore, VegBasePath, referenceLineShp)\n",
    "else:\n",
    "    print('Transects already exist and were loaded')\n",
    "    with open(TransectSpec[:-3]+'pkl', 'rb') as Tfile: \n",
    "        TransectGDF = pickle.load(Tfile)\n",
    "    \n",
    "# make new transect intersections folder\n",
    "if os.path.isdir(os.path.join(filepath, sitename, 'intersections')) is False:\n",
    "    os.mkdir(os.path.join(filepath, sitename, 'intersections'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5426d3a8",
   "metadata": {},
   "source": [
    "### Transect-Veg Intersections\n",
    "Create (or load existing) a GeoDataFrame holding intersections with veg edges along each transect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b13143be",
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.isfile(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_intersects.pkl')):\n",
    "    print('Transect Intersect GDF exists and was loaded')\n",
    "    with open(os.path.join\n",
    "              (filepath , sitename, 'intersections', sitename + '_transect_intersects.pkl'), 'rb') as f:\n",
    "        TransectInterGDF = pickle.load(f)\n",
    "else:\n",
    "    # Get intersections\n",
    "    TransectInterGDF = Transects.GetIntersections(BasePath, TransectGDF, VeglineGDF)\n",
    "    # Save newly intersected transects as shapefile\n",
    "    TransectInterGDF = Transects.SaveIntersections(TransectInterGDF, VeglineGDF, BasePath, sitename)\n",
    "    # Repopulate dict with intersection distances along transects normalised to transect midpoints\n",
    "    TransectInterGDF = Transects.CalculateChanges(TransectInterGDF)\n",
    "    \n",
    "    with open(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(TransectInterGDF, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27ffb709",
   "metadata": {},
   "source": [
    "### Transect-Water Intersections\n",
    "Create (or load existing) a GeoDataFrame holding intersections with waterlines along each transect. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26c15322",
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.isfile(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_water_intersects.pkl')):\n",
    "    print('Transect Intersect + Water GDF exists and was loaded')\n",
    "    with open(os.path.join\n",
    "              (filepath , sitename, 'intersections', sitename + '_transect_water_intersects.pkl'), 'rb') as f:\n",
    "        TransectInterGDFWater = pickle.load(f)\n",
    "else:        \n",
    "    if settings['wetdry'] == True:\n",
    "        TransectInterGDFWater = Transects.GetWaterIntersections(BasePath, TransectGDF, TransectInterGDF, WaterlineGDF, settings, output)  \n",
    "    \n",
    "    with open(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_water_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(TransectInterGDFWater, f)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18a5bbcd",
   "metadata": {},
   "source": [
    "### Transect-Waves Intersections\n",
    "Create (or load existing) a GeoDataFrame holding intersections with topographic data along each transect. This is for comparing veg edge positions with nearshore wave conditions at the time the image was taken. NOTE: this requires you to have a Copernicus Marine Service (CMEMS) account with access to their hindcast model, as you will be asked for a username and password. This should also be run before the tidal corrections and beach with steps if you want to include runup in your corrections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aba62430",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if os.path.isfile(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_wave_intersects.pkl')):\n",
    "    print('Transect Intersect + Wave GDF exists and was loaded')\n",
    "    with open(os.path.join\n",
    "              (filepath , sitename, 'intersections', sitename + '_transect_wave_intersects.pkl'), 'rb') as f:\n",
    "        TransectInterGDFWave = pickle.load(f)\n",
    "else:\n",
    "    TransectInterGDFWave = Transects.WavesIntersect(settings, TransectInterGDF, BasePath, output, lonmin, lonmax, latmin, latmax)\n",
    "    \n",
    "    with open(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_wave_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(TransectInterGDFWave, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b47cc8",
   "metadata": {},
   "source": [
    "### Additional wave-based WL metrics\n",
    "This is for comparing shoreline change with vegetation change, and for quantifying the beach width between the two for each image. If you would like to include runup in your waterline corrections, add `TransectInterGDFWave` to `GetWaterIntersections()`:\n",
    "\n",
    "```TransectInterGDFWater = Transects.GetWaterIntersections(BasePath, TransectGDF, TransectInterGDF, WaterlineGDF, settings, output, TransectInterGDFWave, beachslope)```\n",
    "\n",
    "If you want to include runup AND calculate slopes using CoastSat.slope (recommended), exclude the `beachslope` variable:\n",
    "\n",
    "`TransectInterGDFWater = Transects.GetWaterIntersections(BasePath, TransectGDF, TransectInterGDF, WaterlineGDF, settings, output, TransectInterGDFWave)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "783dd34d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if 'wlcorrdist' not in TransectInterGDFWater.columns:\n",
    "    # Tidal correction to get corrected distances along transects\n",
    "    TransectInterGDFWater = Transects.WLCorrections(settings, output, TransectInterGDFWater, TransectInterGDFWave)     \n",
    "    # Calculate width between VE and corrected WL\n",
    "    TransectInterGDFWater = Transects.CalcBeachWidth(settings, TransectGDF, TransectInterGDFWater)\n",
    "    # Calculate rates of change on corrected WL and save as Transects shapefile\n",
    "    TransectInterGDFWater = Transects.SaveWaterIntersections(TransectInterGDFWater, WaterlineGDF,  BasePath, sitename)\n",
    "    with open(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_water_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(TransectInterGDFWater, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1baed9ef",
   "metadata": {},
   "source": [
    "### Transect-Topo Intersections\n",
    "Create (or load existing) a GeoDataFrame holding intersections with topographic data along each transect. This is for comparing veg edge positions with dune slopes. <font color='red'>Edit the TIF filename</font> if you want to run this intersection process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cb2ca9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDIT ME: Path to slope raster for extracting slope values\n",
    "TIF = '/path/to/Slope_Raster.tif'\n",
    "\n",
    "if os.path.isfile(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_topo_intersects.pkl')):\n",
    "    print('Transect Intersect + Topo GDF exists and was loaded')\n",
    "    with open(os.path.join\n",
    "              (filepath , sitename, 'intersections', sitename + '_transect_topo_intersects.pkl'), 'rb') as f:\n",
    "        TransectInterGDFTopo = pickle.load(f)\n",
    "else:\n",
    "    # Update Transects with Transition Zone widths and slope if available\n",
    "    TransectInterGDFTopo = Transects.TZIntersect(settings, TransectInterGDF, VeglineGDF, BasePath)\n",
    "    TransectInterGDFTopo = Transects.SlopeIntersect(settings, TransectInterGDFTopo, VeglineGDF, BasePath, TIF)\n",
    "    \n",
    "    with open(os.path.join(filepath, sitename, 'intersections', sitename + '_transect_topo_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(TransectInterGDFTopo, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05dd13b8",
   "metadata": {},
   "source": [
    "### <font color='red'>Plotting</font>\n",
    "Edit the desired Transect IDs to plot timeseries of veg change and beach width change across "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e43780b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Timeseries Plotting\n",
    "\n",
    "# EDIT ME: Select transect ID to plot\n",
    "# You can plot subplots within a list of plot IDs, e.g. [[sub1, sub2], plot2]\n",
    "# You can also comment Line 1 out and uncomment Line 2 to create plots for ALL Transect IDs\n",
    "# NOTE: If you want to plot ALL transects, it's recommended you switch ShowPlot=False\n",
    "\n",
    "TransectIDs = [[25,30,35],50,75] # Line 1\n",
    "# TransectIDs = list(TransectInterGDF['TransectID']) # Line 2\n",
    "\n",
    "for TransectID in TransectIDs:\n",
    "    # Plot timeseries of cross-shore veg position\n",
    "    Plotting.VegTimeseries(sitename, TransectInterGDF, TransectID, Hemisphere='N', ShowPlot=True)\n",
    "    # If plotting veg and water lines together\n",
    "    if settings['wetdry']:\n",
    "        Plotting.VegWaterTimeseries(sitename, TransectInterGDFWater, TransectID, Hemisphere='N', ShowPlot=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9263e1d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beach Width Plotting\n",
    "\n",
    "# Select transect ID to plot\n",
    "TransectIDs = [[25,30,35],50,75]\n",
    "for TransectID in TransectIDs:\n",
    "    # Plot timeseries of cross-shore width between water edge and veg edge \n",
    "    Plotting.WidthTimeseries(sitename, TransectInterGDFWater, TransectID, Hemisphere='N')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3e9a24d",
   "metadata": {},
   "source": [
    "## OPTIONAL: Validation\n",
    "### <font color='red'>Validation Settings</font>\n",
    "Most likely you won't need to validate your lines, but if you do, edit these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09fedccb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Name of date column in validation edges shapefile (case sensitive!) \n",
    "DatesCol = 'Date'\n",
    "ValidationShp = './Validation/StAndrews_Veg_Edge_combined_2007_2022_singlepart.shp'\n",
    "# EDIT ME: List of transect ID tuples (startID, finishID)\n",
    "TransectIDList = [(595,711),(726,889),(972,1140),(1141,1297)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2df6249",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Satellite Edges Validation\n",
    "validpath = os.path.join(os.getcwd(), 'Data', sitename, 'validation')\n",
    "\n",
    "if os.path.isfile(os.path.join(validpath, sitename + '_valid_intersects.pkl')):\n",
    "    print('ValidDict exists and was loaded')\n",
    "    with open(os.path.join(validpath, sitename + '_valid_intersects.pkl'), 'rb') as f:\n",
    "        ValidInterGDF = pickle.load(f)\n",
    "else:\n",
    "    ValidInterGDF = Transects.ValidateSatIntersects(sitename, ValidationShp, DatesCol, TransectGDF, TransectInterGDF)\n",
    "    with open(os.path.join(validpath, sitename + '_valid_intersects.pkl'), 'wb') as f:\n",
    "        pickle.dump(ValidInterGDF, f)\n",
    "\n",
    "\n",
    "# Quantify errors between validation and satellite derived lines\n",
    "for TransectIDs in TransectIDList:\n",
    "    Toolbox.QuantifyErrors(sitename, VeglineGDF,'dates',ValidInterGDF,TransectIDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27154581",
   "metadata": {},
   "source": [
    "### <font color='red'>Validation Plots</font>\n",
    "Plot violin plots, probability density functions and regression lines of the cross-shore distance along each transect between each satellite-derived veg edge and validation veg edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9b119ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EDIT ME: List of transect ID tuples (startID, finishID)\n",
    "TransectIDList = [(0,1741)]\n",
    "\n",
    "for TransectIDs in TransectIDList:\n",
    "    PlotTitle = 'Accuracy of Transects ' + str(TransectIDs[0]) + ' to ' + str(TransectIDs[1])\n",
    "    PlottingSeaborn.SatViolin(sitename,VeglineGDF,'dates',ValidInterGDF,TransectIDs, PlotTitle)\n",
    "    PlottingSeaborn.SatPDF(sitename,VeglineGDF,'dates',ValidInterGDF,TransectIDs, PlotTitle)\n",
    "    Plotting.SatRegress(sitename,VeglineGDF,'dates',ValidInterGDF,TransectIDs, PlotTitle)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}