notebooks/zonalstats/wradlib_zonalstats_classes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "This notebook is part of the $\\omega radlib$ documentation: http://wradlib.org/wradlib-docs.\n",
    "\n",
    "Copyright (c) 2018, $\\omega radlib$ developers.\n",
    "Distributed under the MIT License. See LICENSE.txt for more info."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Zonal Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [wradlib.zonalstats](http://wradlib.org/wradlib-docs/latest/zonalstats.html) module provides classes and functions for calculation of zonal statistics for data on arbitrary grids and projections.\n",
    "\n",
    "It provides classes for:\n",
    "\n",
    "- managing georeferenced data (grid points or grid polygons, zonal polygons),\n",
    "- calculation of geographic intersections and managing resulting vector data\n",
    "- calculation of zonal statistics and managing result data as vector attributes\n",
    "- output to vector and raster files available within ogr/gdal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import wradlib as wrl\n",
    "import matplotlib.pyplot as pl\n",
    "import matplotlib as mpl\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "try:\n",
    "    get_ipython().magic(\"matplotlib inline\")\n",
    "except:\n",
    "    pl.ion()\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DataSource"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [wradlib.zonalstats.DataSource](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.DataSource.html) class handles point or polygon vector data by wrapping ogr.DataSource with special functions.\n",
    "\n",
    "The following example shows how to create different DataSource objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import osr\n",
    "\n",
    "# create gk zone 2 projection osr object\n",
    "proj_gk2 = osr.SpatialReference()\n",
    "proj_gk2.ImportFromEPSG(31466)\n",
    "\n",
    "# Setting up DataSource\n",
    "box0 = np.array([[2600000., 5630000.],[2600000., 5640000.],\n",
    "                 [2610000., 5640000.],[2610000., 5630000.],\n",
    "                 [2600000., 5630000.]])\n",
    "box1 = np.array([[2610000., 5630000.],[2610000., 5640000.],\n",
    "                 [2620000., 5640000.],[2620000., 5630000.],\n",
    "                 [2610000., 5630000.]])\n",
    "box2 = np.array([[2600000., 5640000.],[2600000., 5650000.],\n",
    "                 [2610000., 5650000.],[2610000., 5640000.],\n",
    "                 [2600000., 5640000.]])\n",
    "box3 = np.array([[2610000., 5640000.],[2610000., 5650000.],\n",
    "                 [2620000., 5650000.],[2620000., 5640000.],\n",
    "                 [2610000., 5640000.]])\n",
    "\n",
    "point0 = np.array(wrl.georef.get_centroid(box0))\n",
    "point1 = np.array(wrl.georef.get_centroid(box1))\n",
    "point2 = np.array(wrl.georef.get_centroid(box2))\n",
    "point3 = np.array(wrl.georef.get_centroid(box3))\n",
    "\n",
    "# creates Polygons in Datasource\n",
    "poly = wrl.zonalstats.DataSource(np.array([box0, box1, box2, box3]), srs=proj_gk2, name='poly')\n",
    "\n",
    "# creates Points in Datasource\n",
    "point = wrl.zonalstats.DataSource(np.vstack((point0, point1, point2, point3)),\n",
    "                                  srs=proj_gk2, name='point')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the data, which will be exported as numpy arrays. The property ``data`` exports all available data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(poly.data)\n",
    "print(point.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently data can also be retrieved by:\n",
    "\n",
    "- index - [wradlib.zonalstats.DataSource.get_data_by_idx()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.DataSource.get_data_by_idx.html),\n",
    "- attribute - [wradlib.zonalstats.DataSource.get_data_by_att()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.DataSource.get_data_by_att.html) and\n",
    "- geometry - [wradlib.zonalstats.DataSource.get_data_by_geom()](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.DataSource.get_data_by_geom.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with the DataSource being created, we can add/set attribute data of the features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add attribute\n",
    "poly.set_attribute('mean', np.array([10.1, 20.2, 30.3, 40.4]))\n",
    "point.set_attribute('mean', np.array([10.1, 20.2, 30.3, 40.4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attributes associated with features can also be retrieved:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get attributes\n",
    "print(poly.get_attributes(['mean']))\n",
    "# get attributes filtered\n",
    "print(poly.get_attributes(['mean'], filt=('index', 2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can export the contained data to OGR/GDAL supported [vector](http://www.gdal.org/ogr_formats.html) and [raster](http://www.gdal.org/formats_list.html) files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dump as 'ESRI Shapefile', default\n",
    "poly.dump_vector('test_poly.shp')\n",
    "point.dump_vector('test_point.shp')\n",
    "# dump as 'GeoJSON'\n",
    "poly.dump_vector('test_poly.geojson', 'GeoJSON')\n",
    "point.dump_vector('test_point.geojson', 'GeoJSON')\n",
    "# dump as 'GTiff', default\n",
    "poly.dump_raster('test_poly_raster.tif', attr='mean', pixel_size=100.)\n",
    "# dump as 'netCDF'\n",
    "poly.dump_raster('test_poly_raster.nc', 'netCDF', attr='mean', pixel_size=100.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ZonalData"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ZonalData is usually available as georeferenced regular gridded data. Here the [wradlib.zonalstats.ZonalDataBase](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataBase.html) class manages the grid data, the zonal data (target polygons) and the intersection data of source grid and target polygons.\n",
    "Because the calculation of intersection is different for point grids and polygon grids, we have subclasses [wradlib.zonalstats.ZonalDataPoly](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataPoly.html) and [wradlib.zonalstats.ZonalDataPoint](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataPoint.html).\n",
    "\n",
    "Basically, [wradlib.zonalstats.ZonalDataBase](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataBase.html) encapsulates three [wradlib.zonalstats.DataSource](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.DataSource.html) objects:\n",
    "\n",
    "- source grid (points/polygons)\n",
    "- target polygons\n",
    "- destination (intersection) (points/polygons)\n",
    "\n",
    "The destination DataSource object is created from the provided source grid and target polygons at initialisation time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example the creation of a [wradlib.zonalstats.ZonalDataPoly](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataPoly.html) class instance is shown:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup test grid and catchment\n",
    "lon = 7.071664\n",
    "lat = 50.730521\n",
    "r = np.array(range(50, 100*1000 + 50 , 100))\n",
    "a = np.array(range(0, 90, 1))\n",
    "rays = a.shape[0]\n",
    "bins = r.shape[0]\n",
    "\n",
    "# setup OSR objects\n",
    "proj_gk = osr.SpatialReference()\n",
    "proj_gk.ImportFromEPSG(31466)\n",
    "\n",
    "# create polar grid polygon vertices in gk\n",
    "radar_gk = wrl.georef.spherical_to_polyvert(r, a, 0, (lon, lat), \n",
    "                                            proj=proj_gk)\n",
    "radar_gk = radar_gk[..., 0:2]\n",
    "# reshape\n",
    "radar_gk.shape = (rays * bins, 5, 2)\n",
    "\n",
    "box0 = np.array([[2600000., 5630000.],[2600000., 5640000.],\n",
    "                 [2610000., 5640000.],[2610000., 5630000.],\n",
    "                 [2600000., 5630000.]])\n",
    "\n",
    "box1 = np.array([[2610000., 5630000.],[2610000., 5640000.],\n",
    "                 [2620000., 5640000.],[2620000., 5630000.],\n",
    "                 [2610000., 5630000.]])\n",
    "\n",
    "targets = np.array([box0, box1])\n",
    "\n",
    "zdpoly = wrl.zonalstats.ZonalDataPoly(radar_gk, targets, srs=proj_gk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When calculationg the intersection, also weights are calculated for every source grid feature and attributed to the destination features.\n",
    "\n",
    "With the property ``isecs`` it is possible to retrieve the intersection geometries as numpy array, further get-functions add to the functionality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get intersections as numpy array\n",
    "isecs = zdpoly.isecs\n",
    "# get intersections for target polygon 0\n",
    "isec0 = zdpoly.get_isec(0)\n",
    "# get source indices referring to target polygon 0\n",
    "ind0 = zdpoly.get_source_index(0)\n",
    "\n",
    "print(isecs.shape, isec0.shape, ind0.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are import/export functions using [ESRI-Shapfile Format](https://de.wikipedia.org/wiki/Shapefile) as data format. Next export and import is shown:\n",
    "[wradlib.zonalstats.ZonalDataBase](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalDataBase.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zdpoly.dump_vector('test_zdpoly')\n",
    "zdpoly_new = wrl.zonalstats.ZonalDataPoly('test_zdpoly')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ZonalStats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ZonalStats the [wradlib.zonalstats.ZonalStatsBase](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.ZonalStatsBase.html) class and the two subclasses [wradlib.zonalstats.GridCellsToPoly](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.GridCellsToPoly.html) and [wradlib.zonalstats.GridPointsToPoly](http://wradlib.org/wradlib-docs/latest/generated/wradlib.zonalstats.GridPointsToPoly.html) are available. ZonalStatsBase encapsulates one ZonalData object. Properties for simple access of ZonalData, intersection indices and weights are provided. The following code will add ``mean`` and ``var`` attributes to the target DataSource:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create GridCellsToPoly instance\n",
    "gc = wrl.zonalstats.GridCellsToPoly(zdpoly_new)\n",
    "# create some artificial data for processing using the features indices\n",
    "count = radar_gk.shape[0]\n",
    "data = 1000000. / np.array(range(count))\n",
    "# calculate mean and variance\n",
    "mean = gc.mean(data)\n",
    "var = gc.var(data)\n",
    "\n",
    "print(\"Average:\", mean)\n",
    "print(\"Variance:\", var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can export the resulting zonal statistics to vector and raster files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export to vector GeoJSON\n",
    "gc.zdata.trg.dump_vector('test_zonal_json.geojson', 'GeoJSON')\n",
    "# export 'mean' to raster netCDF\n",
    "gc.zdata.trg.dump_raster('test_zonal_hdr.nc', 'netCDF', 'mean', pixel_size=100.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ZonalStats classes can also be used without any ZonalData by instantiating with precalculated index and weight values. Be sure to use matching ix, w and data arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ix, and weight arrays\n",
    "ix = gc.ix\n",
    "w = gc.w\n",
    "# instantiate new ZonlaStats object\n",
    "gc1 = wrl.zonalstats.GridCellsToPoly(ix=ix, w=w)\n",
    "# caclulate statistics\n",
    "avg = gc1.mean(data)\n",
    "var = gc1.var(data)\n",
    "\n",
    "print(\"Average:\", avg)\n",
    "print(\"Variance:\", var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examples of using Zonal Statistics working with rectangular grids as well as polar grids is shown in the [Zonal Statistics Example](wradlib_zonalstats_example.ipynb) notebook."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}