example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "affb366d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "from IPython.display import Markdown\n",
    "from geeCenterline import *\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3b2f498",
   "metadata": {},
   "source": [
    "# Centerline extraction & River/Sediment bar migration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f01b5cdf",
   "metadata": {},
   "source": [
    "## Centerline extraction of the Tallahatchie River\n",
    "This notebook provides an example of how to extract the centerline of river from remote sensing images. In this example, we are going to use the surface reflectance images provided by Landsat and PlanetScope. To concentrate on the workflow of centerline extraction, we have already upload images we need to assets of Google Earth Engine. However, users of Google Earth Engine can obtain Landsat collection directly. Instructions to upload Landsat images through Google Earth Engine is in [here](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db7d24d7",
   "metadata": {},
   "source": [
    "### Import images\n",
    "AssetIDs of images going to be uploaded are required and images should includes at least four bands: blue, green, red, and near-infrared, and the short-wave near-infrared band is also recommended to be included if avaliable. Names of these bands should be uniformed as 'Red', 'Green', 'Blue', 'NIR', 'SWIR1'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "efffbb09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# assits ids of images and region of interest\n",
    "PS_id = 'users/luoyee1997/example/Tallahatchie_Planet'\n",
    "LS_id = 'users/luoyee1997/example/Tallahatchie_Landsat2'\n",
    "roi_id = 'projects/ee-alpha-luoyee1997/assets/LT'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9b6092d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# band names of PlanetScope\n",
    "bands_in_PS = ['b1', 'b2', 'b3', 'b4']\n",
    "# required names\n",
    "bandnames_PS = ['Blue', 'Green', 'Red', 'NIR']\n",
    "# band names of Landsat\n",
    "bands_in_LS = ['SR_B4', 'SR_B3', 'SR_B2', 'SR_B5', 'SR_B6']\n",
    "# required names\n",
    "bandnames_LS = ['Red', 'Green', 'Blue', 'NIR', 'SWIR1']\n",
    "roi = load_file(roi_id)\n",
    "imgLS = load_file(LS_id, bands_in_LS, bandnames_LS).clip(roi)\n",
    "imgPS = load_file(PS_id, bands_in_PS, bandnames_PS).clip(roi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b3f9989",
   "metadata": {},
   "source": [
    "### Obtain river mask from Landsat image\n",
    "Comparing with PlanetScope, Landsat has lower resolution, which means lower computation. We use the Landsat image and the river mask obtained from it to filter out most unwanted water surfaces from the region of interest. In this section, we are going to\n",
    "1. Identify water surfaces.\n",
    "2. Identify the river from water surfaces.\n",
    "3. Create a mask for the PlanetScope image to locate the river."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6c093d4",
   "metadata": {},
   "source": [
    "#### Classification\n",
    "**thrs** of the classification function are the thresholds to classify water, bare soil and vegetation from the NDVI grey image, where values less than thr2 are considered as water. The value between thr3 and thr4 are considered as bare soil, and the value larger than thr5 are considered as vegetation. One thing should be noticed is that we are not going to use these threshold to classify the image. Instead, this classification are going to be used to select training set for the classifer. Thus, these thresholds are not accurite values. The threshold can be modified if the classifer does not perform well, but the default values are robust basing on our tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d4b7e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# landsat classification\n",
    "classLS = classification(imgLS, thr2=0, thr3=0.2, thr4=0.3, thr5=0.6, numPts=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1021c04",
   "metadata": {},
   "source": [
    "#### River mask\n",
    "The river is identified and the mask is dilated for PlanetScope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b8c913e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "waterLS = classLS.eq(1).clip(roi)\n",
    "\n",
    "# connect the water masks divided by small gaps. The radius is the half scale of the gap to be filled.\n",
    "waterLS1 = close(waterLS, radius=1.5, kernelType='square')\n",
    "\n",
    "# identify the river from other water masks. The second input is the minimum area that is considered as the river.\n",
    "# The second input is the connectivity of pixels. 4 means 4-connectivity and 8 means 8-connectivity.\n",
    "riverLS = noise_removal(waterLS1, 500, 8)\n",
    "\n",
    "# Dilate the river mask\n",
    "riverMask = mask_hrs(riverLS, r=1.5, i=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27196c87",
   "metadata": {},
   "source": [
    "### Obtain river mask from the PlanetScope image\n",
    "This section demonstrates how to mask the clasified PlanetScope image and obtain a simply connected river mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "27640f83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# planet classification\n",
    "classPS = classification(imgPS, thr2=0, thr3=0.2, thr4=0.3, thr5=0.6, numPts=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fd160281",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain the river mask\n",
    "waterPS = classPS.eq(1).updateMask(riverMask)\n",
    "waterPS2 = close(waterPS, radius=3)\n",
    "riverPS = noise_removal(waterPS2, 1000, 4)\n",
    "\n",
    "# fill holes of river masks\n",
    "notRiver = riverPS.unmask().Not()\n",
    "riverPS1 = noise_removal(notRiver, 50, 4)\n",
    "riverPS2 = riverPS1.unmask().Not()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "547a9a83",
   "metadata": {},
   "source": [
    "### Centerline extraction"
   ]
  },
  {
   "cell_type": "raw",
   "id": "480f0593",
   "metadata": {},
   "source": [
    "# the location to output the centerline. PLEASE use our own space!\n",
    "out_loc = 'users/luoyee1997/example/river_mask'\n",
    "task = ee.batch.Export.image.toAsset(**{\n",
    "    'image': riverPS2,\n",
    "    'description': 'river_mask',\n",
    "    'assetId': out_loc,\n",
    "    'scale': 3,\n",
    "    'maxPixels': 6e8,\n",
    "  'region': roi\n",
    "})\n",
    "task.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f69485f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "riverMask = ee.Image('users/luoyee1997/example/river_mask').gt(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4a7f36c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate centerline, thre is the threshold to decide wether a value will be considered as centerline.\n",
    "cl = CalculateCenterline(riverPS2.unmask(), thre=0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c19d4c1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# trim unwanted branches, the second input is the longest distance that a line will be considered as branches.\n",
    "# The third input is the number of iteration to do the trimming algrithm.\n",
    "trimed = noise_removal(cl, 100, 8)\n",
    "trimed1 = CleanCenterline(trimed, 300, 3, iterate=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e26d5b3",
   "metadata": {},
   "source": [
    "### Out put\n",
    "The result is going to be output as an asset in Google Earth Engine. This process may take a long time. Please check the process of the test at the Code Editor of Google Earthe Engine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ca4f78c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the location to output the centerline. PLEASE use our own space!\n",
    "out_loc = 'users/luoyee1997/example/trimed2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b6fc883e",
   "metadata": {},
   "outputs": [],
   "source": [
    "scl = trimed1.projection().nominalScale().getInfo()\n",
    "task = ee.batch.Export.image.toAsset(**{\n",
    "    'image': trimed1,\n",
    "    'description': 'river_centerline_trimed',\n",
    "    'assetId': out_loc,\n",
    "    'scale': scl,\n",
    "    'maxPixels': 6e8,\n",
    "    'region': roi\n",
    "})\n",
    "task.start()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26ef1b08",
   "metadata": {},
   "source": [
    "### Visualize the image\n",
    "Visualizing image can help you to check whether the parameters you selected is suitable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9e80719d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map = geemap.Map()\n",
    "vis = {\n",
    "    'bands': ['Red', 'Green', 'Blue'],\n",
    "    'min': 0,\n",
    "    'max': 0.7042,\n",
    "    'gamma': 2\n",
    "}\n",
    "Map.centerObject(imgLS, 12)\n",
    "Map.addLayer(imgLS, vis, 'LS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "05d342b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# classified Landsat iamge\n",
    "Map.addLayer(classLS, {'min':1, 'max':3, 'palette': ['#0080FF', 'red', 'green']}, 'classLS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1451ffc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the river mask obtained from Landsat\n",
    "Map.addLayer(riverMask.selfMask(), {}, 'mask')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f2dc50c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# classified PlanetScope image\n",
    "Map.addLayer(classPS, {'min':1, 'max':3, 'palette': ['#0080FF', 'red', 'green']}, 'classPS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "7ea540f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# river mask to compute the centerline.\n",
    "Map.addLayer(riverPS, {}, 'riverPS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "af71c2f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the computed centerlien\n",
    "Map.addLayer(trimed, {'palette': 'white'}, 'cl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d6d2276a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# upload the trimed centerline and visualize it\n",
    "trimed = ee.Image(out_loc)\n",
    "Map.addLayer(trimed, {'palette': 'coral'}, 'trimed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "9979840e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "edcb7b3763c844dca7d50ef11146493b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[33.519505219659685, -90.29477822219337], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75be65cc",
   "metadata": {},
   "source": [
    "# Point bar imgration of the Arkansas River"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c28cd17",
   "metadata": {},
   "source": [
    "## Load images and study area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3bb4e15e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a .txt file with assit IDs of the Landsat images\n",
    "f = 'arkansas.txt'\n",
    "ar_f = open(f).read().split()\n",
    "# ID of the region of interest\n",
    "roi_id = \"users/luoyee1997/geeRiver/Arkansas/roi\"\n",
    "roi = load_file(roi_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9e69168e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# required names\n",
    "bandnames_LS = ['Red', 'Green', 'Blue', 'NIR', 'SWIR1']\n",
    "ar_cl = batch_load_LS(ar_f, bandnames_LS)\n",
    "ar_cl = ar_cl.map(ref_cr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d37fb57",
   "metadata": {},
   "source": [
    "## Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "88e212bb",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "classified = wrap_classification(ar_cl, thr2=0, thr3=0.2, thr4=0.3, thr5=0.6, numPts=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c407c565",
   "metadata": {},
   "source": [
    "## Generate .gif from image collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "fec49ad1",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "[gif]https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/videoThumbnails/f946380aa52149fc109c40bf3c3d3178-20ba835189413c50714e431745aa64fb:getPixels"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "video_args = {\n",
    "    'region': roi,\n",
    "    'framesPerSecond': 2,\n",
    "    'min': 1,\n",
    "    'max': 3,\n",
    "    'palette': ['skyblue', 'sandybrown', 'green']\n",
    "}\n",
    "\n",
    "# Get URL that will produce the animation when accessed.\n",
    "gif_url1 = classified.getVideoThumbURL(video_args)\n",
    "Markdown('[gif]'+gif_url1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0a3f0c73",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "[gif]https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/videoThumbnails/6331610e94e5fdb363a6ad77ab455000-82fff5f12cb512479b41022a33cb4cb8:getPixels"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "video_args = {\n",
    "  'region': roi,\n",
    "  'framesPerSecond': 2,\n",
    "  'max': 0.8,\n",
    "  'min': 0,\n",
    "  'gamma': [2, 2, 2.5]\n",
    "}\n",
    "\n",
    "# Get URL that will produce the animation when accessed.\n",
    "gif_url2 = ar_cl.select(['SWIR1', 'NIR', 'Green']).getVideoThumbURL(video_args)\n",
    "Markdown('[gif]'+gif_url2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ac24b4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}