LakeTimeSeriesv1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "LakeTimeSeriesv1.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "49jEY_NDOVQo",
        "2X5UWWBhOhMw"
      ],
      "authorship_tag": "ABX9TyMxFVN1ejzL4ayzWEqcgxVH",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ericslevenson/arctic-surface-water/blob/main/LakeTimeSeriesv1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#**Description:** \n",
        "This script ingests a shapefile of individual waterbodies and exports a csv of water, cloudy, and total pixels for each day of valid observations. Observations through time of individual waterbodies are merged outside of Earth Engine (script under development).\n",
        "\n",
        "#**Directions:** \n",
        "- Run the 'Setup' chunk and to authenticate your EE account with google colab. Currently the google drive and colab user authentication are commented out, but for development purposes sometimes it's helpful to include those.\n",
        "- Set your input and export variables and run the chunk. \n",
        "- Run the 'Functions' chunk. Note: you can collapse chunks for ease of use.\n",
        "- Run the 'Main' chunk.\n",
        "- Export progress can be monitored at code.earthengine.google.com under the 'Tasks' panel.\n",
        "\n",
        "#**Potential Errors:**\n",
        "Certain aspects of this workflow may need some tweaking in order to adapt to different analyses and starting shapefiles:\n",
        "- Currently set to run within one UTM zone.\n",
        "- Some of the waterbody properties depend on attributes associated with the imported shapefile, such as centroid, count, and ID. The 'lakeProps' function that combines all of the waterbody property extraction functions will likely fail with some of these attributes missing. In the future I can put in some try/except statements to avoid errors, but for now we likely have to alter the property extraction functions for specific use cases.\n",
        "- Cloudy image filtering currently occurs based on the cloudiness of the entire ROI. This can easily be adjusted to just include cloudiness over the lakes, but filtering per lake will be difficult. One option is to just not do any filtering because cloudy pixels are counted in the exported time series.\n"
      ],
      "metadata": {
        "id": "LGrbMpSWrl0-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Setup"
      ],
      "metadata": {
        "id": "42fyaZMGOqrB"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "vwu9PTkpONFv",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 197
        },
        "outputId": "89d474c7-96b5-45f0-9a4e-62b3c8c4cbd5"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does not start automatically, please manually browse the URL below.\n",
            "\n",
            "    https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=cCA8jC5RFv3gaADI6ZVs1jeoNmloqH0Nr8Vt1wGs3iA&tc=tcZBzFRKJ2ooMFE-wfshAdO3pt1kzfIuwtTberMX7jg&cc=oh-ZzQmffLHBsDpK5-bQQb4v1rD6Vy7AvDVv7MduoNw\n",
            "\n",
            "The authorization workflow will generate a code, which you should paste in the box below. \n",
            "Enter verification code: 4/1AX4XfWjYD4gKQZC6j3EOgWMPKU2add5f3Hbu42_p1JpGvkDr5kbvvHmpldM\n",
            "\n",
            "Successfully saved authorization token.\n",
            "geemap package not installed. Installing ...\n"
          ]
        }
      ],
      "source": [
        "# Authenticate private account (only required for exporting to drive/gee/gcp)\n",
        "#from google.colab import auth \n",
        "#auth.authenticate_user()\n",
        "\n",
        "# Earth Engine setup\n",
        "import ee # Trigger the authentication flow.\n",
        "ee.Authenticate()\n",
        "ee.Initialize() # Initialize the library.\n",
        "\n",
        "# Google Drive setup (if needed)\n",
        "#from google.colab import drive\n",
        "#drive.mount('/content/drive')\n",
        "\n",
        "# Some common imports\n",
        "from IPython.display import Image\n",
        "import folium\n",
        "# Installs geemap package\n",
        "import subprocess\n",
        "\n",
        "try:\n",
        "    import geemap\n",
        "except ImportError:\n",
        "    print('geemap package not installed. Installing ...')\n",
        "    subprocess.check_call([\"python\", '-m', 'pip', 'install', 'geemap']) \n",
        "\n",
        "import geemap"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Input and Output Settings"
      ],
      "metadata": {
        "id": "49jEY_NDOVQo"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "## ***INPUTS***\n",
        "\n",
        "# Time Period Dates\n",
        "start = ee.Date('2020-06-01')\n",
        "finish = ee.Date('2020-09-30')\n",
        "# Region of Interest\n",
        "xmin = -150.2\n",
        "xmax = -150\n",
        "ymin =  69.2\n",
        "ymax =  69.23\n",
        "roi = ee.Geometry.Rectangle([xmin, ymin, xmax, ymax])\n",
        "# UTM zone\n",
        "epsg = 'EPSG:32606'\n",
        "# Lake Shapefile\n",
        "lakes = ee.FeatureCollection('projects/toolik/assets/ToolikLakes1')\n",
        "# Minimum Pixel Count for a lake \n",
        "minPix = 10 \n",
        "# Image scale\n",
        "pixScale = 10\n",
        "\n",
        "# ***EXPORTS***\n",
        "\n",
        "# Export Properties\n",
        "exportSelectors = ['id', 'waterPixels', 'allPixels', 'clearPixels', 'cloudPixels', 'count', 'centroid', 'ratio'] # These will likely need to be changed to reflect your input shapefile. Specifically, 'count', 'centroid', and 'ratio' are all attributes from the imported EE asset. The others can stay the same.\n",
        "# ROI Description\n",
        "roiLabel = 'Toolik_small'\n",
        "# Export Folder\n",
        "exportFolder = 'tooliktimeseries2'"
      ],
      "metadata": {
        "id": "fHHMfvC4OWAI"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Functions"
      ],
      "metadata": {
        "id": "2X5UWWBhOhMw"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "## ***IMAGE PRE-PROCESSING METHODS***\n",
        "\n",
        "# Mask clouds in Sentinel-2\n",
        "def maskS2clouds(image):\n",
        "  '''Takes an input and adds two bands: cloud mask and clear mask'''\n",
        "  qa = image.select('QA60')\n",
        "  cloudBitMask = 1 << 10\n",
        "  cirrusBitMask = 1 << 11\n",
        "  clear_mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(qa.bitwiseAnd(cirrusBitMask).eq(0)).rename('clear_mask')\n",
        "  cloud_mask = qa.bitwiseAnd(cloudBitMask).eq(1).And(qa.bitwiseAnd(cirrusBitMask).eq(1)).rename('cloud_mask')\n",
        "  return image.addBands([cloud_mask,clear_mask])\n",
        "\n",
        "# Clip image\n",
        "def clip_image(image):\n",
        "  '''Clips to the roi defined at the beginning of the script'''\n",
        "  return image.clip(roi)\n",
        "\n",
        "def clip2lakes(image):\n",
        "  '''Clips an image based on the lake boundaries'''\n",
        "  return image.clip(lakes)\n",
        "\n",
        "# Get percentile cover   \n",
        "def getCover(image):\n",
        "  '''calculates percentage of the roi covered by the clear mask. NOTE: this function\n",
        "  calls the global totPixels variable that needs to be calculated in the main script.'''\n",
        "  actPixels = ee.Number(image.updateMask(image.select('clear_mask')).reduceRegion(\n",
        "      reducer = ee.Reducer.count(),\n",
        "      scale = 100,\n",
        "      geometry = roi,\n",
        "      maxPixels=1e12,\n",
        "      ).values().get(0))\n",
        "  # calculate the perc of cover OF CLEAR PIXELS \n",
        "  percCover = actPixels.divide(totPixels).multiply(100).round()\n",
        "  # number as output\n",
        "  return image.set('percCover', percCover,'actPixels',actPixels)\n",
        "  \n",
        "# Mosaic images by date, orbit, - basically combines images together that were taken on the same day \n",
        "def mosaicBy(imcol):\n",
        "  '''Takes an image collection (imcol) and creates a mosaic for each day\n",
        "  Returns: An image collection of daily mosaics'''\n",
        "  #return the collection as a list of images (not an image collection)\n",
        "  imlist = imcol.toList(imcol.size())\n",
        "  # Get all the dates as list\n",
        "  def imdate(im):\n",
        "    date = ee.Image(im).date().format(\"YYYY-MM-dd\")\n",
        "    return date\n",
        "  all_dates = imlist.map(imdate)\n",
        "  # get all orbits as list\n",
        "  def orbitId(im):\n",
        "    orb = ee.Image(im).get('SENSING_ORBIT_NUMBER')\n",
        "    return orb\n",
        "  all_orbits = imlist.map(orbitId)\n",
        "  # get all spacecraft names as list\n",
        "  def spacecraft(im):\n",
        "    return ee.Image(im).get('SPACECRAFT_NAME')\n",
        "  all_spNames = imlist.map(spacecraft)\n",
        "  # this puts dates, orbits and names into a nested list\n",
        "  concat_all = all_dates.zip(all_orbits).zip(all_spNames);\n",
        "  # here we unnest the list with flatten, and then concatenate the list elements with \" \"\n",
        "  def concat(el):\n",
        "    return ee.List(el).flatten().join(\" \")\n",
        "  concat_all = concat_all.map(concat)\n",
        "  # here, just get distinct combintations of date, orbit and name\n",
        "  concat_unique = concat_all.distinct()\n",
        "  # mosaic\n",
        "  def mosaicIms(d):\n",
        "    d1 = ee.String(d).split(\" \")\n",
        "    date1 = ee.Date(d1.get(0))\n",
        "    orbit = ee.Number.parse(d1.get(1)).toInt()\n",
        "    spName = ee.String(d1.get(2))\n",
        "    im = imcol.filterDate(date1, date1.advance(1, \"day\")).filterMetadata('SPACECRAFT_NAME', 'equals', spName).filterMetadata('SENSING_ORBIT_NUMBER','equals', orbit).mosaic()\n",
        "    return im.set(\n",
        "        \"system:time_start\", date1.millis(),\n",
        "        \"system:date\", date1.format(\"YYYY-MM-dd\"),\n",
        "        \"system:id\", d1)\n",
        "  mosaic_imlist = concat_unique.map(mosaicIms)\n",
        "  return ee.ImageCollection(mosaic_imlist)\n",
        "\n",
        "def reprojectMosaic(image):\n",
        "  '''Reproject to UTM. A future function should take the image location and return\n",
        "  the UTM zone. For now, I'm manually entering the EPSG code.'''\n",
        "  image_projected = image.reproject(epsg)\n",
        "  return image_projected\n",
        "\n",
        "###########################################################################\n",
        "## ***WATER CLASSIFICATION METHODS***\n",
        "\n",
        "# Define NDWI image\n",
        "def ndwi(image):\n",
        "  '''Adds an NDWI band to the input image'''\n",
        "  return image.normalizedDifference(['B3', 'B8']).rename('NDWI').multiply(1000)\n",
        "\n",
        "# Basic ndwi classification  \n",
        "def ndwi_classify(image):\n",
        "  '''Creates a binary image based on an NDWI threshold of 0'''\n",
        "  ndwimask = image.select('NDWI')\n",
        "  water = ndwimask.gte(0)\n",
        "  land = ndwimask.lt(0)\n",
        "  return(water)\n",
        "\n",
        "# OTSU thresholding from histogram\n",
        "def otsu(histogram):\n",
        "  '''Returns the NDWI threshold for binary water classification'''\n",
        "  counts = ee.Array(ee.Dictionary(histogram).get('histogram'))\n",
        "  means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'))\n",
        "  size = means.length().get([0])\n",
        "  total = counts.reduce(ee.Reducer.sum(), [0]).get([0])\n",
        "  sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0])\n",
        "  mean = sum.divide(total)\n",
        "  indices = ee.List.sequence(1, size)\n",
        "  def func_xxx(i):\n",
        "    '''Compute between sum of squares, where each mean partitions the data.'''\n",
        "    aCounts = counts.slice(0, 0, i)\n",
        "    aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0])\n",
        "    aMeans = means.slice(0, 0, i)\n",
        "    aMean = aMeans.multiply(aCounts) \\\n",
        "        .reduce(ee.Reducer.sum(), [0]).get([0]) \\\n",
        "        .divide(aCount)\n",
        "    bCount = total.subtract(aCount)\n",
        "    bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount)\n",
        "    return aCount.multiply(aMean.subtract(mean).pow(2)).add(\n",
        "           bCount.multiply(bMean.subtract(mean).pow(2)))\n",
        "  bss = indices.map(func_xxx)\n",
        "  # Return the mean value corresponding to the maximum BSS.\n",
        "  return means.sort(bss).get([-1])\n",
        "\n",
        "# OTSU thresholding for an image\n",
        "def otsu_thresh(water_image):\n",
        "  '''Calculate NDWI and create histogram. Return the OTSU threshold.'''\n",
        "  NDWI = ndwi(water_image).select('NDWI').updateMask(water_image.select('clear_mask'))\n",
        "  histogram = ee.Dictionary(NDWI.reduceRegion(\n",
        "    geometry = roi,\n",
        "    reducer = ee.Reducer.histogram(255, 2).combine('mean', None, True).combine('variance', None, True),\n",
        "    scale = pixScale,\n",
        "    maxPixels = 1e12\n",
        "  ))\n",
        "  return otsu(histogram.get('NDWI_histogram'))\n",
        "\n",
        "# Classify an image using OTSU threshold.\n",
        "def otsu_classify(water_image):\n",
        "  '''(1) Calculate NDWI and create histogram. (2) Calculate NDWI threshold for \n",
        "  binary classification using OTSU method. (3) Classify image and add layer to input image.\n",
        "  '''\n",
        "  NDWI = ndwi(water_image).select('NDWI')\n",
        "  histogram = ee.Dictionary(NDWI.reduceRegion(\n",
        "    geometry = roi,\n",
        "    reducer = ee.Reducer.histogram(255, 2).combine('mean', None, True).combine('variance', None, True),\n",
        "    scale = pixScale,\n",
        "    maxPixels = 1e12\n",
        "  ))\n",
        "  threshold = otsu(histogram.get('NDWI_histogram'))\n",
        "  otsu_classed = NDWI.gt(ee.Number(threshold)).And(water_image.select('B8').lt(2000)).rename('otsu_classed')\n",
        "  return water_image.addBands([otsu_classed])\n",
        "\n",
        "def adaptive_thresholding(water_image):\n",
        "  '''Takes an image clipped to lakes and returns the water mask'''\n",
        "  NDWI = ndwi(water_image).select('NDWI')#.updateMask(water_image.select('clear_mask')) # get NDWI **TURNED OFF CLOUD MASK, SHOULD THIS STAY OFF?**\n",
        "  threshold = ee.Number(otsu_thresh(water_image)) \n",
        "  threshold = threshold.divide(10).round().multiply(10)\n",
        "  # get fixed histogram\n",
        "  histo = NDWI.reduceRegion(\n",
        "      geometry = roi,\n",
        "      reducer = ee.Reducer.fixedHistogram(-1000, 1000, 200),\n",
        "      scale = pixScale, # This was 30, keep at 10!?!?\n",
        "      maxPixels = 1e12\n",
        "  )\n",
        "  hist = ee.Array(histo.get('NDWI'))\n",
        "  counts = hist.cut([-1,1])\n",
        "  buckets = hist.cut([-1,0])\n",
        "  #find split points from otsu threshold\n",
        "  threshold = ee.Array([threshold]).toList()\n",
        "  buckets_list = buckets.toList()\n",
        "  split = buckets_list.indexOf(threshold)\n",
        "  # split into land and water slices\n",
        "  land_slice = counts.slice(0,0,split)\n",
        "  water_slice = counts.slice(0,split.add(1),-1)\n",
        "  # find max of land and water slices\n",
        "  land_max = land_slice.reduce(ee.Reducer.max(),[0])\n",
        "  water_max = water_slice.reduce(ee.Reducer.max(),[0])\n",
        "  land_max = land_max.toList().get(0)\n",
        "  water_max = water_max.toList().get(0)\n",
        "  land_max = ee.List(land_max).getNumber(0)\n",
        "  water_max = ee.List(water_max).getNumber(0)\n",
        "  #find difference between land, water and otsu val\n",
        "  counts_list = counts.toList()\n",
        "  otsu_val = ee.Number(counts_list.get(split))\n",
        "  otsu_val = ee.List(otsu_val).getNumber(0)\n",
        "  land_prom = ee.Number(land_max).subtract(otsu_val)\n",
        "  water_prom = ee.Number(water_max).subtract(otsu_val)\n",
        "  #find land and water buckets corresponding to 0.9 times the prominence\n",
        "  land_thresh = ee.Number(land_max).subtract((land_prom).multiply(ee.Number(0.9)))\n",
        "  water_thresh = ee.Number(water_max).subtract((water_prom).multiply(ee.Number(0.9)))\n",
        "  land_max_ind = land_slice.argmax().get(0)\n",
        "  water_max_ind = water_slice.argmax().get(0)\n",
        "  li = ee.Number(land_max_ind).subtract(1)\n",
        "  li = li.max(ee.Number(1))\n",
        "  wi = ee.Number(water_max_ind).add(1)\n",
        "  wi = wi.min(ee.Number(199))\n",
        "  land_slice2 = land_slice.slice(0,li,-1).subtract(land_thresh)\n",
        "  water_slice2 = water_slice.slice(0,0,wi).subtract(water_thresh)\n",
        "  land_slice2  = land_slice2.abs().multiply(-1)\n",
        "  water_slice2 = water_slice2.abs().multiply(-1)\n",
        "  land_index = ee.Number(land_slice2.argmax().get(0)).add(land_max_ind)\n",
        "  water_index = ee.Number(water_slice2.argmax().get(0)).add(split)\n",
        "  land_level = ee.Number(buckets_list.get(land_index))\n",
        "  water_level = ee.Number(buckets_list.get(water_index))\n",
        "  land_level = ee.Number(ee.List(land_level).get(0)).add(5)\n",
        "  water_level = ee.Number(ee.List(water_level).get(0)).add(5)\n",
        "  #calculate water fraction and classify\n",
        "  water_fraction = (NDWI.subtract(land_level)).divide(water_level.subtract(land_level)).multiply(100).rename('water_fraction')\n",
        "  #water_fraction = conditional(water_fraction) #sets values less than 0 to 0 and greater than 100 to 100\n",
        "  water_75 = water_fraction.gte(75).rename('water_75'); #note, this is a non-binary classification, so we use 75% water as \"water\"\n",
        "  all_mask = water_image.select('B2').gt(5).rename('all_mask')\n",
        "  cloud_mask_ed = water_image.select('clear_mask').neq(1).rename('cloud_mask_ed')\n",
        "  return water_image.addBands([water_fraction,water_75,NDWI,cloud_mask_ed])\n",
        "\n",
        "###############################################################################\n",
        "## ***PROPERTY EXTRACTION METHODS***\n",
        "def sumWater(lake):\n",
        "  '''sums the water pixels within a buffered lake polygon and adds the result to the feature'''\n",
        "  watersum = lakeIm.select('water_75').reduceRegion(\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      geometry = lake.geometry(),\n",
        "      scale = 10,\n",
        "      maxPixels=1e9\n",
        "  ).get('water_75')\n",
        "  return watersum\n",
        "\n",
        "def sumClear(lake):\n",
        "  '''sums the number of clear pixels within a buffered lake polygon'''\n",
        "  clearsum = lakeIm.select('clear_mask').reduceRegion(\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      geometry = lake.geometry(),\n",
        "      scale = 10,\n",
        "      maxPixels=1e9\n",
        "  ).get('clear_mask')\n",
        "  return clearsum\n",
        "\n",
        "def sumClouds(lake):\n",
        "  '''sums the number of clear pixels within a buffered lake polygon'''\n",
        "  cloudsum = lakeIm.select('cloud_mask_ed').reduceRegion(\n",
        "      reducer=ee.Reducer.sum(),\n",
        "      geometry = lake.geometry(),\n",
        "      scale = 10,\n",
        "      maxPixels=1e9\n",
        "  ).get('cloud_mask_ed')\n",
        "  return cloudsum\n",
        "\n",
        "def sumAll(lake):\n",
        "  '''sums the number of pixels within a buffered lake polygon'''\n",
        "  all_mask = lakeIm.select('B2').gt(5).rename('all_mask')\n",
        "  allsum = all_mask.reduceRegion(\n",
        "      reducer=ee.Reducer.count(),\n",
        "      geometry = lake.geometry(),\n",
        "      scale = 10,\n",
        "      maxPixels=1e9\n",
        "      ).get('all_mask')\n",
        "  return allsum\n",
        "\n",
        "def troid(lake):\n",
        "  center = ee.Array(lake.centroid().geometry().coordinates())\n",
        "  return center\n",
        "\n",
        "def getID(lake):\n",
        "  id = ee.Number(lake.id())\n",
        "  return id\n",
        "\n",
        "def getCount(lake):\n",
        "  count = ee.Number(lake.get('count'))\n",
        "  return count\n",
        "\n",
        "def getLabel(lake):\n",
        "  label = ee.String(lake.get('label'))\n",
        "  return label\n",
        "\n",
        "def lakeProps(lake):\n",
        "  water = sumWater(lake)\n",
        "  clear = sumClear(lake)\n",
        "  clouds = sumClouds(lake)\n",
        "  all = sumAll(lake)\n",
        "  centroid = troid(lake)\n",
        "  id = getID(lake)\n",
        "  count = getCount(lake)\n",
        "  label = getLabel(lake)\n",
        "  return ee.Feature(None, {'id': id, 'label':label, 'waterPixels': water, 'clearPixels': clear, 'cloudPixels': clouds, 'allPixels': all, 'centroid': centroid,'count': count})\n",
        "\n",
        "###############################################################################\n",
        "## ***EXPORT METHODS***\n",
        "def export_lakes(collection, description, fileNamePrefix, fileFormat, folder, selectors):\n",
        "  '''Export a feature collection of lake properties to google drive for a given day.'''\n",
        "  task = ee.batch.Export.table.toDrive(**{\n",
        "    'collection': collection,\n",
        "    'description': description, \n",
        "    'fileNamePrefix': fileNamePrefix,\n",
        "    'fileFormat': fileFormat,\n",
        "    'folder': folder,\n",
        "    'selectors': selectors\n",
        "  })\n",
        "  task.start()"
      ],
      "metadata": {
        "id": "IZErfZ2sOfyo"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Main"
      ],
      "metadata": {
        "id": "6AFgm4ybOknX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "###############################################################################\n",
        "## ***Data Pre-Processing***\n",
        "\n",
        "# Filter lakes\n",
        "lakes = lakes.filter(ee.Filter.gt('count', minPix))\n",
        "lakes = lakes.filterBounds(roi)\n",
        "\n",
        "# Get Images, mosaic, clip to ROI, filter by coverage, and clip to lake mask\n",
        "images = ee.ImageCollection('COPERNICUS/S2').filterBounds(roi).filterDate(start,finish).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',50)) # Get Images\n",
        "images_all = mosaicBy(images) # Mosaic images\n",
        "images_all = images_all.map(maskS2clouds) # Create cloud/clear masks\n",
        "# Filter by percentile cover\n",
        "image_mask = images.select('B2').mean().clip(roi).gte(0) #First, calculate total number of pixels\n",
        "totPixels = ee.Number(image_mask.reduceRegion(\n",
        "    reducer = ee.Reducer.count(),\n",
        "    scale = 100,\n",
        "    geometry = roi,\n",
        "    maxPixels = 1e12\n",
        "    ).values().get(0))\n",
        "images_all = images_all.map(clip_image) # Clip to roi\n",
        "images_all = images_all.map(getCover) # Add percent cover as an image property\n",
        "images = images_all.filterMetadata('percCover','greater_than',50) # remove images covering less than 50% of the ROI)\n",
        "lakeimages = images.map(clip2lakes) # Clip images to buffered lake mask\n",
        "dates = lakeimages.aggregate_array('system:date').getInfo() # get a non-EE list of dates so it is iterable.\n",
        "\n",
        "###############################################################################\n",
        "## ***WATER CLASSIFICATION***\n",
        "lakeimages = lakeimages.map(adaptive_thresholding)\n",
        "\n",
        "###############################################################################\n",
        "## ***ITERATE THROUGH DAYS AND EXPORT***\n",
        "for i, date in enumerate(dates):\n",
        "  # Get lake properties\n",
        "  eedate = ee.Date(date)\n",
        "  lakeIm = lakeimages.filterDate(eedate).first()\n",
        "  lakes2 = lakes.map(lakeProps)\n",
        "  # Export\n",
        "  exportDate = date.replace('-', '_')\n",
        "  description = roiLabel +'_'+ exportDate\n",
        "  fileformat = 'CSV'\n",
        "  export_lakes(lakes2, description, description, fileformat, exportFolder, exportSelectors)"
      ],
      "metadata": {
        "id": "yXCa2bkyOmpA"
      },
      "execution_count": 4,
      "outputs": []
    }
  ]
}