{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "<a href=\"https://githubtocolab.com/gee-community/geemap/blob/master/examples/notebooks/water_app.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "# Using Earth Engine and geemap for mapping surface water dynamics\n",
    "\n",
    "**Steps to create Landsat timeseries:**\n",
    "\n",
    "1. Pan and zoom to your area of interest, and click on the map to select a polygon.\n",
    "2. Adjust the parameters (e.g., band combination, threshold, color) if needed.\n",
    "3. Click the `Submit` button to create timeseries of Landsat imagery and normalized difference indices.\n",
    "\n",
    "**Web Apps:** https://gishub.org/water-app, https://gishub.org/water-ngrok\n",
    "\n",
    "**Contact:** Dr. Qiusheng Wu ([Website](https://wetlands.io/), [LinkedIn](https://www.linkedin.com/in/qiushengwu), [Twitter](https://twitter.com/giswqs), [YouTube](https://youtube.com/@giswqs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check geemap installation\n",
    "import subprocess\n",
    "\n",
    "try:\n",
    "    import geemap\n",
    "except ImportError:\n",
    "    print(\"geemap package is not installed. Installing ...\")\n",
    "    subprocess.check_call([\"python\", \"-m\", \"pip\", \"install\", \"geemap\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import os\n",
    "import ee\n",
    "import geemap\n",
    "import ipywidgets as widgets\n",
    "from bqplot import pyplot as plt\n",
    "from ipyleaflet import WidgetControl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an interactive map\n",
    "Map = geemap.Map(center=[37.71, 105.47], zoom=4, add_google_map=False)\n",
    "Map.add_basemap(\"HYBRID\")\n",
    "Map.add_basemap(\"ROADMAP\")\n",
    "\n",
    "# Add Earth Engine data\n",
    "fc = ee.FeatureCollection(\"users/giswqs/public/chn_admin_level2\")\n",
    "Map.addLayer(fc, {}, \"二级行政区\")\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Designe interactive widgets\n",
    "\n",
    "style = {\"description_width\": \"initial\"}\n",
    "\n",
    "output_widget = widgets.Output(layout={\"border\": \"1px solid black\"})\n",
    "output_control = WidgetControl(widget=output_widget, position=\"bottomright\")\n",
    "Map.add_control(output_control)\n",
    "\n",
    "admin1_widget = widgets.Text(\n",
    "    description=\"一级行政区:\", value=\"广东省\", width=200, style=style\n",
    ")\n",
    "\n",
    "admin2_widget = widgets.Text(\n",
    "    description=\"二级行政区:\", value=\"广州市\", width=300, style=style\n",
    ")\n",
    "\n",
    "band_combo = widgets.Dropdown(\n",
    "    description=\"显示影像波段组合:\",\n",
    "    options=[\n",
    "        \"Red/Green/Blue\",\n",
    "        \"NIR/Red/Green\",\n",
    "        \"SWIR2/SWIR1/NIR\",\n",
    "        \"NIR/SWIR1/Red\",\n",
    "        \"SWIR2/NIR/Red\",\n",
    "        \"SWIR2/SWIR1/Red\",\n",
    "        \"SWIR1/NIR/Blue\",\n",
    "        \"NIR/SWIR1/Blue\",\n",
    "        \"SWIR2/NIR/Green\",\n",
    "        \"SWIR1/NIR/Red\",\n",
    "    ],\n",
    "    value=\"NIR/Red/Green\",\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "year_widget = widgets.IntSlider(\n",
    "    min=1984,\n",
    "    max=2020,\n",
    "    value=2010,\n",
    "    description=\"显示指定年份数据:\",\n",
    "    width=400,\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "fmask_widget = widgets.Checkbox(\n",
    "    value=True, description=\"应用fmask(去除云, 阴影, 雪)\", style=style\n",
    ")\n",
    "\n",
    "\n",
    "# Normalized Satellite Indices: https://www.usna.edu/Users/oceano/pguth/md_help/html/norm_sat.htm\n",
    "\n",
    "nd_options = [\n",
    "    \"Vegetation Index (NDVI)\",\n",
    "    \"Water Index (NDWI)\",\n",
    "    \"Modified Water Index (MNDWI)\",\n",
    "    \"Snow Index (NDSI)\",\n",
    "    \"Soil Index (NDSI)\",\n",
    "    \"Burn Ratio (NBR)\",\n",
    "    \"Customized\",\n",
    "]\n",
    "nd_indices = widgets.Dropdown(\n",
    "    options=nd_options,\n",
    "    value=\"Modified Water Index (MNDWI)\",\n",
    "    description=\"归一化指数:\",\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "first_band = widgets.Dropdown(\n",
    "    description=\"波段1:\",\n",
    "    options=[\"Blue\", \"Green\", \"Red\", \"NIR\", \"SWIR1\", \"SWIR2\"],\n",
    "    value=\"Green\",\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "second_band = widgets.Dropdown(\n",
    "    description=\"波段2:\",\n",
    "    options=[\"Blue\", \"Green\", \"Red\", \"NIR\", \"SWIR1\", \"SWIR2\"],\n",
    "    value=\"SWIR1\",\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "nd_threshold = widgets.FloatSlider(\n",
    "    value=0,\n",
    "    min=-1,\n",
    "    max=1,\n",
    "    step=0.01,\n",
    "    description=\"阈值:\",\n",
    "    orientation=\"horizontal\",\n",
    "    style=style,\n",
    ")\n",
    "\n",
    "nd_color = widgets.ColorPicker(\n",
    "    concise=False, description=\"颜色:\", value=\"blue\", style=style\n",
    ")\n",
    "\n",
    "\n",
    "def nd_index_change(change):\n",
    "    if nd_indices.value == \"Vegetation Index (NDVI)\":\n",
    "        first_band.value = \"NIR\"\n",
    "        second_band.value = \"Red\"\n",
    "    elif nd_indices.value == \"Water Index (NDWI)\":\n",
    "        first_band.value = \"NIR\"\n",
    "        second_band.value = \"SWIR1\"\n",
    "    elif nd_indices.value == \"Modified Water Index (MNDWI)\":\n",
    "        first_band.value = \"Green\"\n",
    "        second_band.value = \"SWIR1\"\n",
    "    elif nd_indices.value == \"Snow Index (NDSI)\":\n",
    "        first_band.value = \"Green\"\n",
    "        second_band.value = \"SWIR1\"\n",
    "    elif nd_indices.value == \"Soil Index (NDSI)\":\n",
    "        first_band.value = \"SWIR1\"\n",
    "        second_band.value = \"NIR\"\n",
    "    elif nd_indices.value == \"Burn Ratio (NBR)\":\n",
    "        first_band.value = \"NIR\"\n",
    "        second_band.value = \"SWIR2\"\n",
    "    elif nd_indices.value == \"Customized\":\n",
    "        first_band.value = None\n",
    "        second_band.value = None\n",
    "\n",
    "\n",
    "nd_indices.observe(nd_index_change, names=\"value\")\n",
    "\n",
    "submit = widgets.Button(\n",
    "    description=\"Submit\", button_style=\"primary\", tooltip=\"Click me\", style=style\n",
    ")\n",
    "\n",
    "full_widget = widgets.VBox(\n",
    "    [\n",
    "        widgets.HBox([admin1_widget, admin2_widget]),\n",
    "        widgets.HBox([band_combo, year_widget, fmask_widget]),\n",
    "        widgets.HBox([nd_indices, first_band, second_band, nd_threshold, nd_color]),\n",
    "        submit,\n",
    "    ]\n",
    ")\n",
    "\n",
    "full_widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Capture user interaction with the map\n",
    "\n",
    "\n",
    "def handle_interaction(**kwargs):\n",
    "    latlon = kwargs.get(\"coordinates\")\n",
    "    if kwargs.get(\"type\") == \"click\":\n",
    "        Map.default_style = {\"cursor\": \"wait\"}\n",
    "        xy = ee.Geometry.Point(latlon[::-1])\n",
    "        selected_fc = fc.filterBounds(xy)\n",
    "\n",
    "        with output_widget:\n",
    "            output_widget.clear_output()\n",
    "\n",
    "            try:\n",
    "                admin1_id = selected_fc.first().get(\"ADM1_ZH\").getInfo()\n",
    "                admin2_id = selected_fc.first().get(\"ADM2_ZH\").getInfo()\n",
    "                admin1_widget.value = admin1_id\n",
    "                admin2_widget.value = admin2_id\n",
    "                Map.layers = Map.layers[:4]\n",
    "                geom = selected_fc.geometry()\n",
    "                layer_name = admin1_id + \"-\" + admin2_id\n",
    "                Map.addLayer(\n",
    "                    ee.Image().paint(geom, 0, 2), {\"palette\": \"red\"}, layer_name\n",
    "                )\n",
    "                print(layer_name)\n",
    "            except:\n",
    "                print(\"找不到相关行政区\")\n",
    "                Map.layers = Map.layers[:4]\n",
    "\n",
    "        Map.default_style = {\"cursor\": \"pointer\"}\n",
    "\n",
    "\n",
    "Map.on_interaction(handle_interaction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Click event handler\n",
    "\n",
    "\n",
    "def submit_clicked(b):\n",
    "    with output_widget:\n",
    "        output_widget.clear_output()\n",
    "        print(\"Computing...\")\n",
    "        Map.default_style = {\"cursor\": \"wait\"}\n",
    "\n",
    "        try:\n",
    "            admin1_id = admin1_widget.value\n",
    "            admin2_id = admin2_widget.value\n",
    "            band1 = first_band.value\n",
    "            band2 = second_band.value\n",
    "            selected_year = year_widget.value\n",
    "            threshold = nd_threshold.value\n",
    "            bands = band_combo.value.split(\"/\")\n",
    "            apply_fmask = fmask_widget.value\n",
    "            palette = nd_color.value\n",
    "\n",
    "            roi = fc.filter(\n",
    "                ee.Filter.And(\n",
    "                    ee.Filter.eq(\"ADM1_ZH\", admin1_id),\n",
    "                    ee.Filter.eq(\"ADM2_ZH\", admin2_id),\n",
    "                )\n",
    "            )\n",
    "            Map.layers = Map.layers[:4]\n",
    "            geom = roi.geometry()\n",
    "            layer_name = admin1_id + \"-\" + admin2_id\n",
    "            Map.addLayer(ee.Image().paint(geom, 0, 2), {\"palette\": \"red\"}, layer_name)\n",
    "\n",
    "            images = geemap.landsat_timeseries(\n",
    "                roi=roi,\n",
    "                start_year=1984,\n",
    "                end_year=2020,\n",
    "                start_date=\"01-01\",\n",
    "                end_date=\"12-31\",\n",
    "                apply_fmask=apply_fmask,\n",
    "            )\n",
    "            nd_images = images.map(lambda img: img.normalizedDifference([band1, band2]))\n",
    "            result_images = nd_images.map(lambda img: img.gt(threshold))\n",
    "\n",
    "            selected_image = ee.Image(\n",
    "                images.toList(images.size()).get(selected_year - 1984)\n",
    "            )\n",
    "            selected_result_image = ee.Image(\n",
    "                result_images.toList(result_images.size()).get(selected_year - 1984)\n",
    "            ).selfMask()\n",
    "\n",
    "            vis_params = {\"bands\": bands, \"min\": 0, \"max\": 3000}\n",
    "\n",
    "            Map.addLayer(selected_image, vis_params, \"Landsat \" + str(selected_year))\n",
    "            Map.addLayer(\n",
    "                selected_result_image,\n",
    "                {\"palette\": palette},\n",
    "                \"Result \" + str(selected_year),\n",
    "            )\n",
    "\n",
    "            def cal_area(img):\n",
    "                pixel_area = img.multiply(ee.Image.pixelArea()).divide(1e6)\n",
    "                img_area = pixel_area.reduceRegion(\n",
    "                    **{\n",
    "                        \"geometry\": roi.geometry(),\n",
    "                        \"reducer\": ee.Reducer.sum(),\n",
    "                        \"scale\": 1000,\n",
    "                        \"maxPixels\": 1e12,\n",
    "                        \"bestEffort\": True,\n",
    "                    }\n",
    "                )\n",
    "                return img.set({\"area\": img_area})\n",
    "\n",
    "            areas = result_images.map(cal_area)\n",
    "            stats = areas.aggregate_array(\"area\").getInfo()\n",
    "            x = list(range(1984, 2021))\n",
    "            y = [item.get(\"nd\") for item in stats]\n",
    "\n",
    "            fig = plt.figure(1)\n",
    "            fig.layout.height = \"270px\"\n",
    "            plt.clear()\n",
    "            plt.plot(x, y)\n",
    "            plt.title(\"Temporal trend (1984-2020)\")\n",
    "            plt.xlabel(\"Year\")\n",
    "            plt.ylabel(\"Area (km2)\")\n",
    "\n",
    "            output_widget.clear_output()\n",
    "\n",
    "            plt.show()\n",
    "\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            print(\"An error occurred during computation.\")\n",
    "\n",
    "        Map.default_style = {\"cursor\": \"default\"}\n",
    "\n",
    "\n",
    "submit.on_click(submit_clicked)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}