{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "from bqplot import pyplot as plt\n",
    "from bqplot import Bars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map = geemap.Map()\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "roi = ee.Geometry.Point([-99.09805297851562, 47.10016154728056])\n",
    "# set visualization parameters\n",
    "vis = {\n",
    "    \"nir\": {\"bands\": [\"N\", \"R\", \"G\"]},\n",
    "    \"rgb\": {\"bands\": [\"R\", \"G\", \"B\"]},\n",
    "    \"water\": {\"palette\": [\"0000FF\"]},\n",
    "    \"ndwi\": {\n",
    "        \"palette\": [\n",
    "            \"#ece7f2\",\n",
    "            \"#d0d1e6\",\n",
    "            \"#a6bddb\",\n",
    "            \"#74a9cf\",\n",
    "            \"#3690c0\",\n",
    "            \"#0570b0\",\n",
    "            \"#045a8d\",\n",
    "            \"#023858\",\n",
    "        ]\n",
    "    },\n",
    "    \"ndvi\": {\n",
    "        \"palette\": [\n",
    "            \"#d73027\",\n",
    "            \"#f46d43\",\n",
    "            \"#fdae61\",\n",
    "            \"#fee08b\",\n",
    "            \"#d9ef8b\",\n",
    "            \"#a6d96a\",\n",
    "            \"#66bd63\",\n",
    "            \"#1a9850\",\n",
    "        ]\n",
    "    },\n",
    "}\n",
    "\n",
    "# search NAIP imagery that has RGBN bands\n",
    "collection = (\n",
    "    ee.ImageCollection(\"USDA/NAIP/DOQQ\")\n",
    "    .filterBounds(roi)\n",
    "    .filterDate(\"2009-01-01\", \"2019-12-31\")\n",
    "    .filter(ee.Filter.listContains(\"system:band_names\", \"N\"))\n",
    ")\n",
    "# print(collection.getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = collection.first()\n",
    "Map.centerObject(image)\n",
    "Map.addLayer(image, vis[\"nir\"], \"NAIP\")\n",
    "polygon = ee.Geometry(image.geometry())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the histogram of the NIR band.  The mean and variance are only FYI.\n",
    "histogram = image.select(\"N\").reduceRegion(\n",
    "    **{\n",
    "        \"reducer\": ee.Reducer.histogram(255, 2),\n",
    "        \"geometry\": polygon,\n",
    "        \"scale\": 10,\n",
    "        \"bestEffort\": True,\n",
    "    }\n",
    ")\n",
    "hist_dict = histogram.getInfo()\n",
    "# print(hist_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = hist_dict[\"N\"][\"bucketMeans\"]\n",
    "y = hist_dict[\"N\"][\"histogram\"]\n",
    "plt.bar(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Return the DN that maximizes interclass variance in B5 (in the region).\n",
    "def otsu(histogram):\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",
    "\n",
    "    indices = ee.List.sequence(1, size)\n",
    "\n",
    "    # Compute between sum of squares, where each mean partitions the data.\n",
    "\n",
    "    def func_xxx(i):\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 = (\n",
    "            aMeans.multiply(aCounts)\n",
    "            .reduce(ee.Reducer.sum(), [0])\n",
    "            .get([0])\n",
    "            .divide(aCount)\n",
    "        )\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",
    "        )\n",
    "\n",
    "    bss = indices.map(func_xxx)\n",
    "\n",
    "    # Return the mean value corresponding to the maximum BSS.\n",
    "    return means.sort(bss).get([-1])\n",
    "\n",
    "\n",
    "threshold = otsu(histogram.get(\"N\"))\n",
    "print(\"threshold\", threshold.getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "classA = image.select(\"N\").lt(threshold)\n",
    "Map.addLayer(classA.mask(classA), {\"palette\": \"blue\"}, \"class A\")\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_water(image):\n",
    "    histogram = image.select(\"N\").reduceRegion(\n",
    "        **{\n",
    "            \"reducer\": ee.Reducer.histogram(255, 2),\n",
    "            \"geometry\": polygon,\n",
    "            \"scale\": 10,\n",
    "            \"bestEffort\": True,\n",
    "        }\n",
    "    )\n",
    "    threshold = otsu(histogram.get(\"N\"))\n",
    "    water = image.select(\"N\").lt(threshold).selfMask()\n",
    "    return water.set({\"threshold\": threshold})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "water_images = collection.map(extract_water)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(water_images.first(), {\"palette\": \"blue\"}, \"first water\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "water_images.aggregate_array(\"threshold\").getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}