{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map = geemap.Map()\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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,
   "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,
   "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,
   "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,
   "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,
   "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,
   "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "water_images = collection.map(extract_water)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map.addLayer(water_images.first(), {\"palette\": \"blue\"}, \"first water\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "water_images.aggregate_array(\"threshold\").getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Map"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}