update zonal statistics

esciencecenter-digital-skills · Oct 9, 2023 · 3a21221 · 3a21221
1 parent 0b0bbea
commit 3a21221
Showing 1 changed file with 50 additions and 11 deletions.
diff --git a/notebooks/11-zonal-statistics.ipynb b/notebooks/11-zonal-statistics.ipynb
@@ -28,6 +28,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Load assets vector\n",
     "import geopandas as gpd\n",
     "assets = gpd.read_file('assets.gpkg')\n",
     "assets"
@@ -97,6 +98,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# The image has 1 band dimension\n",
     "burned.shape"
    ]
   },
@@ -129,6 +131,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Rasterize the vector data.\n",
+    "# transform represents the projection from pixel space to the projected coordinate space. \n",
+    "# By default, the pixels that are not contained within a polygon in our shapefile will be filled with 0.\n",
     "from rasterio import features\n",
     "assets_rasterized = features.rasterize(geom, out_shape=burned_squeeze.shape, transform=burned.rio.transform())"
    ]
@@ -225,36 +230,70 @@
    "source": [
     "Exercise: Zonal stats over slope zones\n",
     "\n",
-    "In the previouse section, we calculated the slope over Rhodes Island. Can you perform zonal statistics over slope and non slope region? We define non slope area as `slope<0.3`.\n",
+    "Now let's look into the burned areas in relation to slope classes like [Zhai et al., 2023](https://www.mdpi.com/1999-4907/14/4/807) have done. To reproduce their analysis, we will:\n",
+    "1. Reclassify the slope classes into 5 categories  0%–5%, 5%–10%, 10%–15%, 15%–25%, and above 25%, then,\n",
+    "2. Perform the zonal statistics on the above categories.\n",
     "\n",
     "Hint:\n",
     "1. Load slope data from 'slope_dask.tif'.\n",
-    "2. Convert slope data to zones using [`xr.where`](https://docs.xarray.dev/en/stable/generated/xarray.where.html) function.\n",
-    "3. Note that the slope data and burn index data are in different resolution."
+    "2. The big chellenge will be how to compute the slope zones. Consider:\n",
+    "    - Convert slope data to zones using [`numpy.digitize`](https://numpy.org/doc/stable/reference/generated/numpy.digitize.html).\n",
+    "    - Use [`xarray.apply_ufunc`](https://docs.xarray.dev/en/stable/generated/xarray.apply_ufunc.html) to apply `numpy.digitize` to all the elements in an Xarray. \n",
+    "4. Note that the slope data and burn index data are in different resolution."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0c9e89d8",
+   "id": "6f2f7d87",
    "metadata": {},
    "outputs": [],
    "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "\n",
     "# Load data and remove redundant dimension\n",
-    "slope = rioxarray.open_rasterio('slope_dask.tif').squeeze()\n",
+    "# slope = rioxarray.open_rasterio('slope.tif').squeeze()\n",
+    "slope = rioxarray.open_rasterio('slope_dask.tif', masked=True).squeeze()\n",
     "\n",
-    "# Covert slope data to Zone data\n",
-    "import xarray as xr\n",
-    "slope2 = xr.where(slope.isnull(), 0, slope) # Set all NaN values to 0\n",
-    "slope_zones = xr.where(slope2 < 0.3, 0, 1)\n",
-    "slope_zones.rio.set_crs(slope.rio.crs)\n",
+    "# Defines the bins for pixel values\n",
+    "# Zone 0 will be values <=0.; Zone 5 will be values >0.25; Zone 6 will be NaN values(>1000)\n",
+    "slope_bins = (0., 0.05, 0.1, 0.15, 0.25, 1000)\n",
+    "\n",
+    "slope_zones = xr.apply_ufunc(\n",
+    "    np.digitize,\n",
+    "    slope,\n",
+    "    slope_bins\n",
+    ")\n",
     "\n",
     "# Reproject burned to slope\n",
+    "# Because slope zones has lower resolution\n",
     "burned = rioxarray.open_rasterio('burned.tif').squeeze()\n",
     "burned_match = burned.rio.reproject_match(slope_zones)\n",
     "\n",
     "# Compute Zonal stats\n",
-    "stats = zonal_stats(slope_zones, burned_match)"
+    "stats = zonal_stats(slope_zones, burned_match)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4cab5e96",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Visualize slope zones\n",
+    "slope_zones.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "95d2cc47",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stats"
    ]
   }
  ],