notebooks/zonal_module.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "291d02ac-38a9-435e-b2c0-204a787dda2e",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/UM-RMRS/raster_tools/blob/main/notebooks/zonal_module.ipynb)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "69f0b07e-bdad-4954-8538-67549332a6da",
   "metadata": {},
   "source": [
    "# Raster Tools Zonal Module\n",
    "## This notebook demonstrates the functionality of the Raster Tools zonal module \n",
    "by John Hogland 4/27/2023"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "605b128b-73a9-4eb7-b1f9-e2efec2eb7c6",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "16ffc7e8-3215-418c-8db9-66ab27d7ee23",
   "metadata": {},
   "source": [
    "# Install software for Colab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a84740a8-7ff0-4c8b-a1c8-53be5ec2320a",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --upgrade gdown\n",
    "!pip install --upgrade numba\n",
    "!pip install --upgrade geopandas\n",
    "!pip install mapclassify\n",
    "!pip install --upgrade datascience\n",
    "!pip install --upgrade gym\n",
    "!pip install --upgrade folium\n",
    "!pip install raster_tools\n",
    "!pip install py3dep"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b8108a3a-6ee9-46e2-9151-08d071d7ae8a",
   "metadata": {},
   "source": [
    "# The Process\n",
    "In this notebook we create some raster datasets and use the zonal module to summarize and extract values\n",
    "\n",
    "## Steps\n",
    "- 1. Create random raster surfaces\n",
    "- 2. Create zonal Raster and summarize\n",
    "- 3. Create zonal Vector and summarize\n",
    "- 4. Create random locations\n",
    "- 5. Extract raster values from those locations"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1a931bc6-6cb7-4624-b496-cf4f5964fb6c",
   "metadata": {},
   "source": [
    "## Step 1: Create a random raster surface\n",
    "### Import various packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d27e667-e1eb-47cf-a8e9-b74d1369e3c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from raster_tools import Raster, general, zonal\n",
    "import py3dep\n",
    "import geopandas as gpd\n",
    "import numpy as np"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0a012218-5497-4c71-90dc-84b6871b38dd",
   "metadata": {},
   "source": [
    "### Create a 3 band raster surface (y=1500 by x=1000 cells) of random numbers between 0 and 255\n",
    "#### This raster will be the raster values that get summarized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a3dcc09-ba89-44ce-951a-22a409606764",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "rs = Raster(np.random.randint(0, 255, (3, 1500, 1000)))\n",
    "display(rs.xdata)\n",
    "rs.plot(col=\"band\", col_wrap=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "311c2cd1",
   "metadata": {},
   "source": [
    "## Step 2: Create zonal Raster and summarize\n",
    "This raster will be what we use to summarize values (dimensions = 1,1500,1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "509c333c",
   "metadata": {},
   "outputs": [],
   "source": [
    "vls = np.ones((1, 1500, 1000))\n",
    "vls[0, 750:1500, 0:500] = 2\n",
    "vls[0, 0:750, 0:500] = 3\n",
    "vls[0, 750:1500, 500:1000] = 4\n",
    "z_rs = Raster(vls).astype(int)\n",
    "display(z_rs.xdata)\n",
    "z_rs.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "28d742a4",
   "metadata": {},
   "source": [
    "### Summarize by zones (Raster)\n",
    "Note, the output is organized in \"wide\" format by default with the zone ID as the index.\n",
    "The columns are a multi-level index with the band label on top and the stat names below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92c062a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonal.zonal_stats(z_rs, rs, list(zonal.ZONAL_STAT_FUNCS)).compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f293bf75-0bad-4ede-b108-30dabb7b3fbe",
   "metadata": {},
   "source": [
    "Pass `wide_format=False` to get a \"long\" format data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9278af8c-faeb-47bd-969e-28b6704c6624",
   "metadata": {},
   "outputs": [],
   "source": [
    "zonal.zonal_stats(\n",
    "    z_rs, rs, list(zonal.ZONAL_STAT_FUNCS), wide_format=False\n",
    ").compute()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "31d1c3b9",
   "metadata": {},
   "source": [
    "### Step 3: Create a zonal Vector and summarize\n",
    "- Create a Vector object with a Albers project ('EPSG:5070'). Note Vector and Raster projection must be specified and match. To match projection use geopandas to_crs() function.\n",
    "- Overlapping areas are accounted for when using Vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c51b1ac4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely import geometry\n",
    "\n",
    "xmin, ymin, xmax, ymax = z_rs.bounds\n",
    "x, y = (xmin, ymin)\n",
    "lngx = (xmax - xmin) / 2\n",
    "lngy = (ymax - ymin) / 2\n",
    "geom_lst = []\n",
    "while y < ymax:\n",
    "    while x < xmax:\n",
    "        geom = geometry.Polygon(\n",
    "            [\n",
    "                (x, y),\n",
    "                (x, y + lngy),\n",
    "                (x + lngx, y + lngy),\n",
    "                (x + lngx, y),\n",
    "                (x, y),\n",
    "            ]\n",
    "        )\n",
    "        geom_lst.append(geom)\n",
    "        x += lngx\n",
    "\n",
    "    x = xmin\n",
    "    y += lngy\n",
    "\n",
    "vct = zonal.get_vector(\n",
    "    gpd.GeoDataFrame(\n",
    "        {\"label\": [5, 6, 7, 8], \"geometry\": geom_lst}, crs=\"EPSG:5070\"\n",
    "    )\n",
    ")  # specifying a arbitrary projection for the example\n",
    "\n",
    "# visualize the data\n",
    "dt = vct.data.compute()  # compute data to a geopandas dataframe\n",
    "display(dt)\n",
    "dt.plot(facecolor=\"none\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b5eeb957",
   "metadata": {},
   "source": [
    "### Summarize values using the zonal Vector layer\n",
    "- Outputs are organized by zone in the index and by band label and stat names in the columns.\n",
    "- By default, the zone IDs are the index of the Vector object plus 1.\n",
    "- Use `features_field` to select a zone ID field."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f54dd2fc-8ef6-475e-a4da-558e46883210",
   "metadata": {},
   "source": [
    "With the default zone IDs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985bc2aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "odf = zonal.zonal_stats(\n",
    "    vct, rs.set_crs(\"EPSG:5070\"), list(zonal.ZONAL_STAT_FUNCS)\n",
    ").compute()  # don't forget to set a projection for the raster dataset\n",
    "odf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c0ab13-db5a-4bfb-bc1e-c455a1cccc71",
   "metadata": {},
   "source": [
    "With a zone IDs selected from the Vector's data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca57d214-3643-435f-b18d-551f99ecfd59",
   "metadata": {},
   "outputs": [],
   "source": [
    "odf = zonal.zonal_stats(\n",
    "    vct,\n",
    "    rs.set_crs(\"EPSG:5070\"),\n",
    "    list(zonal.ZONAL_STAT_FUNCS),\n",
    "    features_field=\"label\",\n",
    ").compute()\n",
    "odf"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7da9ffae",
   "metadata": {},
   "source": [
    "## Step 4: Create 150 random points to extract raster values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d62c67",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 150\n",
    "xmin, ymin, xmax, ymax = z_rs.bounds\n",
    "xdif = xmax - xmin\n",
    "ydif = ymax - ymin\n",
    "pnts_lst = []\n",
    "while len(pnts_lst) < n:\n",
    "    x = (np.random.random() * xdif) + xmin\n",
    "    y = (np.random.random() * ydif) + ymin\n",
    "    pnt = geometry.Point([x, y])\n",
    "    pnts_lst.append(pnt)\n",
    "\n",
    "dic = {\"geometry\": pnts_lst}\n",
    "pnts = gpd.GeoDataFrame(dic)\n",
    "# Visualize the points\n",
    "display(pnts)\n",
    "pnts.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "696e6f70",
   "metadata": {},
   "source": [
    "## Step 5: Extract raster values using the point locations\n",
    "- By default extracted values are returned in a one to many dataframe. The dataframe contains 2 columns and a row for each point band combination (pnts*bands)\n",
    "- Alternatively, you can specify the output to be merged by column by setting the optional variable axis to 1\n",
    "- We will set the optional variable making a new dataset ready to join with the pnts dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4e34a35",
   "metadata": {},
   "outputs": [],
   "source": [
    "edf = zonal.extract_points_eager(\n",
    "    pnts, rs, column_name=\"test\", axis=1\n",
    ").compute()\n",
    "edf"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "eb11e94f",
   "metadata": {},
   "source": [
    "### Joint edf to pnts dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a645ce00",
   "metadata": {},
   "outputs": [],
   "source": [
    "pnts_edf = pnts.join(edf)\n",
    "display(pnts_edf)\n",
    "pnts_edf.plot(column=\"test_1\", legend=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "95a52ee7-9618-49b7-b912-5e13c75bd693",
   "metadata": {
    "tags": []
   },
   "source": [
    "# This ends the Raster Tools zonal module notebook\n",
    "## Check out the other notebooks:\n",
    "- https://github.com/UM-RMRS/raster_tools/blob/main/notebooks/README.md\n",
    "## References\n",
    "- Raster-Tools GitHub: https://github.com/UM-RMRS/raster_tools\n",
    "- Hogland's Spatial Solutions: https://sites.google.com/view/hoglandsspatialsolutions/home\n",
    "- Dask: https://dask.org/\n",
    "- Geopandas:https://geopandas.org/en/stable/\n",
    "- Xarray: https://docs.xarray.dev/en/stable/\n",
    "- Jupyter: https://jupyter.org/\n",
    "- Anaconda:https://www.anaconda.com/\n",
    "- VS Code: https://code.visualstudio.com/\n",
    "- ipywidgets: https://ipywidgets.readthedocs.io/en/latest/\n",
    "- numpy:https://numpy.org/\n",
    "- matplotlib:https://matplotlib.org/\n",
    "- folium: https://python-visualization.github.io/folium/\n",
    "- pandas: https://pandas.pydata.org/\n",
    "- sklearn: https://scikit-learn.org/stable/index.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RStools Python 3.10",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}