Test rio-tiler coverage-weighted stats

rio_tiler_zonalstats_reader.py

"""
Adapted from the rio-tiler docs https://cogeotiff.github.io/rio-tiler/advanced/statistics/#area-weighted-statistics
"""

import attr
from typing import Any, Union, Optional, List, Dict

from rio_tiler import io
from rio_tiler.models import BandStatistics
from rio_tiler.constants import WGS84_CRS

from rasterio.crs import CRS

from geojson_pydantic.features import Feature, FeatureCollection
from geojson_pydantic.geometries import Polygon


class Reader(io.Reader):
    """Custom Reader with zonal_statistics method."""

    def zonal_statistics(
        self,
        geojson: Union[FeatureCollection, Feature],
        categorical: bool = False,
        categories: Optional[List[float]] = None,
        percentiles: Optional[List[int]] = None,
        hist_options: Optional[Dict] = None,
        max_size: int = None,
        align_bounds_with_dataset: bool = True,
        shape_crs: Optional[CRS] = WGS84_CRS,
        cover_scale: int = 10,
        **kwargs: Any,
    ) -> Union[FeatureCollection, Feature]:
        """Return statistics from GeoJSON features.

        Args:
            geojson (Feature or FeatureCollection): a GeoJSON Feature or FeatureCollection.
            categorical (bool): treat input data as categorical data. Defaults to False.
            categories (list of numbers, optional): list of categories to return value for.
            percentiles (list of numbers, optional): list of percentile values to calculate. Defaults to `[2, 98]`.
            hist_options (dict, optional): Options to forward to numpy.histogram function.
            max_size (int, optional): Limit the size of the longest dimension of the dataset read, respecting bounds X/Y aspect ratio. Defaults to None.
            kwargs (optional): Options to forward to `self.preview`.

        Returns:
            Feature or FeatureCollection

        """
        kwargs = {**self.options, **kwargs}

        hist_options = hist_options or {}

        fc = geojson
        # We transform the input Feature to a FeatureCollection
        if isinstance(fc, Feature):
            fc = FeatureCollection(type="FeatureCollection", features=[geojson])

        for feature in fc:
            geom = feature.model_dump(exclude_none=True)

            # Get data overlapping with the feature (using Reader.feature method)
            data = self.feature(
                geom,
                shape_crs=shape_crs,
                align_bounds_with_dataset=align_bounds_with_dataset,
                max_size=max_size,
                **kwargs,
            )
            coverage_array = data.get_coverage_array(
                geom,
                shape_crs=shape_crs,
                cover_scale=cover_scale,
            )

            stats = data.statistics(
                categorical=categorical,
                categories=categories,
                percentiles=percentiles,
                hist_options=hist_options,
                coverage=coverage_array,
            )
            
            yield coverage_array, stats

            # Update input feature properties and add the statistics
            feature.properties = feature.properties or {}
            feature.properties.update({"statistics": stats})

        #return fc.features[0] if isinstance(geojson, Feature) else fc

test-rio-tiler-coverage-weighted-stats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c7d0663c-6fd3-4ea5-b805-149b94000361",
   "metadata": {},
   "source": [
    "# Validate zonal statistics calculation with `rio-tiler`\n",
    "\n",
    "In this test, we use rio-tiler's functionality to calculate a `coverage` array and apply that in rio-tiler's statistics calculation, to get **pixel-geometry intersection area weighted** zonal statistics.\n",
    "\n",
    "We choose an equal-area map projection for this case. How the method handles a non-equal-area projection (by on-the-fly reprojecting) is left to other tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7cfdd8f1-87bf-4cfe-b617-08f7d29f2648",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#!pip install --upgrade rio_cogeo rio_tiler geojson_pydantic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "652018dd-92e9-4f54-ae5c-ce93a848eef6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import rasterio\n",
    "import rasterio.warp\n",
    "from rasterio.windows import Window\n",
    "from rasterio.transform import Affine\n",
    "from rasterio.crs import CRS\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from geojson_pydantic.features import Feature\n",
    "\n",
    "from rio_tiler_zonalstats_reader import Reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fb53720c-ec27-436b-98e3-8362086e24b4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_path_cea = \"row_encoded_raster_cea.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cb39e3e3-6578-416e-bba7-973fa3b54346",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "width, height = 10, 10\n",
    "\n",
    "encoded_values = np.arange((width * height)).reshape(width, height)\n",
    "\n",
    "# Cylindrical Equal Area projection\n",
    "crs = CRS.from_epsg(6933)\n",
    "pixel_size = 2e6 / width\n",
    "transform = Affine.translation(0, 0) * Affine.scale(pixel_size, -pixel_size)\n",
    "\n",
    "#crs = CRS.from_epsg(4326)\n",
    "#transform = Affine.translation(0, 0) * Affine.scale(20000, -20000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6dbb9364-e951-48ac-a55a-b2738b36972a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],\n",
       "       [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],\n",
       "       [20, 21, 22, 23, 24, 25, 26, 27, 28, 29],\n",
       "       [30, 31, 32, 33, 34, 35, 36, 37, 38, 39],\n",
       "       [40, 41, 42, 43, 44, 45, 46, 47, 48, 49],\n",
       "       [50, 51, 52, 53, 54, 55, 56, 57, 58, 59],\n",
       "       [60, 61, 62, 63, 64, 65, 66, 67, 68, 69],\n",
       "       [70, 71, 72, 73, 74, 75, 76, 77, 78, 79],\n",
       "       [80, 81, 82, 83, 84, 85, 86, 87, 88, 89],\n",
       "       [90, 91, 92, 93, 94, 95, 96, 97, 98, 99]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "encoded_values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "149272ff-b681-4a2f-9a70-a27a2fc70a10",
   "metadata": {},
   "source": [
    "Sum up row by row progressively increasing the number of pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "488b9d31-6a35-48ec-a35a-645512032ea7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3217.5, 50.0, 64.35)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weighted_sum = 0.0\n",
    "weighted_count = 0.0\n",
    "\n",
    "for i, row in enumerate(encoded_values):\n",
    "    weighted_sum += np.sum(row[:i]) + row[i] * 0.5\n",
    "    weighted_count += len(row[:i]) + 0.5\n",
    "    \n",
    "weighted_sum, weighted_count, (weighted_sum / weighted_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9099ed7a-ca23-4330-a26d-08aec6533afb",
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(\n",
    "    raster_path_cea, \"w\", driver=\"GTiff\",\n",
    "    height=height, width=width,\n",
    "    count=1, dtype=str(encoded_values.dtype),\n",
    "    crs=crs,\n",
    "    transform=transform,\n",
    ") as dst:\n",
    "    dst.write(encoded_values, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "38a9443b-3246-41ed-842b-e95888ce289e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd4AAAGdCAYAAAC8UhIBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAj7ElEQVR4nO3de2xUdRr/8c90gGmr06qYTttQoCRNQKrhUkIEFDZKkxWMhP15Ay/rZQNpQSqJXBYvlV06AbVpfnatKdmQuqTIL1Ejm3ih0VhkkVgLqEED2ZXARG26GtJWLq2dOb8/sLOO5TLTM3Pm29P3K/n+0cM5cx5OC0+f53vO+Xosy7IEAAAckZHuAAAAGElIvAAAOIjECwCAg0i8AAA4iMQLAICDSLwAADiIxAsAgINIvAAAOGiU0yeMRCL67rvv5Pf75fF4nD49AMAGy7LU09OjwsJCZWSkrnY7f/68+vr6bH/OmDFjlJmZmYSIksfxxPvdd9+pqKjI6dMCAJIoFApp3LhxKfns8+fPq3jC1eroDNv+rPz8fJ04ccKo5Ot44vX7/ZKkebpDozTa6dNfWoY33REM4hnt+Lfnijxe82YnPKPMu07ymvfzJK9518kzysTrZN7PuAy6Tv2RPn0U2h79vzwV+vr61NEZ1on2CcrxD/370d0TUfHMk+rr6xvZiXegvTxKozXKY1Di9Zjzgz3A4zHwP0quU3wMvE7KMO86eQz8hdfEX8JNjMmJqcIcf4atxGsq8/4lAgAgKWxFFLaxjE/YiiQvmCQi8QIAjBSRpYiGnnntHJtKJF4AgJEiishOzWrv6NRxX/McAACDUfECAIwUtiyFraG3i+0cm0okXgCAkdw6x0urGQAAB1HxAgCMFJGlsAsrXhIvAMBItJoBAIBtVLwAACO59a7mIVW8r7zyioqLi5WZmamZM2fq448/TnZcAIARLpKEYaKEE+/u3btVVVWlTZs26fDhw7rlllv0+9//XqdOnUpFfAAAuErCibe2tlaPPfaYHn/8cU2ZMkV1dXUqKipSQ0NDKuIDAIxQ4V/uarYzTJTQHG9fX5/a29u1YcOGmO3l5eU6cODARY/p7e1Vb29v9Ovu7u4hhAkAGGnClmyuTpS8WJIpoYr3hx9+UDgcViAQiNkeCATU0dFx0WOCwaByc3Ojo6ioaOjRAgBGDOZ4f+W3CyBblnXJRZE3btyorq6u6AiFQkM5JQAArpBQq/n666+X1+sdVN12dnYOqoIH+Hw++Xy+oUcIABiRIvIorIsXdfEeb6KEKt4xY8Zo5syZamlpidne0tKiOXPmJDUwAMDIFrHsDxMl/AKNtWvX6sEHH1RZWZluvvlmNTY26tSpU1q5cmUq4gMAwFUSTrz33nuvfvzxR23evFnff/+9SktL9c4772jChAmpiA8AMEKFbbaa7RybSkN6ZWRFRYUqKiqSHQsAAFFuTbwskgAAgINYJAEAYKSI5VHEsnFXs41jU4nECwAwEq1mAABgGxUvAMBIYWUobKM+DCcxlmQi8QIAjGTZnOO1mOMFACB+zPECAADbqHgBAEYKWxkKWzbmeN3yrmYAAJwQkUcRG43ZiMzMvLSaAQBwUNoq3lHF4zUqw5x1eq1R3nSHMJiBMVkZ5v2uZo0yMCaveTd1WF4Tr5OJMRn4vTPoZ7y//7x00plzufXmKlrNAAAj2Z/jpdUMAMCIR8ULADDShZurbCySQKsZAID4RWy+MpK7mgEAABUvAMBMbr25isQLADBSRBmufIEGiRcAYKSw5VHYxgpDdo5NJeZ4AQBwEBUvAMBIYZt3NYdpNQMAEL+IlaGIjZurIobeXEWrGQAAB1HxAgCMRKsZAAAHRWTvzuRI8kJJKlrNAAA4iIoXAGAk+y/QMLO2JPECAIxk/5WRZiZeM6MCAMClqHgBAEZiPV4AABzk1lYziRcAYCT7z/GamXjNjAoAAJei4gUAGClieRSx8wINQ5cFJPECAIwUsdlqNvU5XjOjAgDApah4AQBGsr8soJm1JYkXAGCksDwK23gW186xqWTmrwMAALgUFS8AwEhubTWbGRUAYMQL63/t5qGNxPT39+vpp59WcXGxsrKyNGnSJG3evFmRyP9W9rUsS9XV1SosLFRWVpYWLFigo0ePJnQeEi8AAJK2bt2qV199VfX19fr666+1bds2vfDCC3r55Zej+2zbtk21tbWqr69XW1ub8vPztXDhQvX09MR9HlrNAAAjOd1q/uSTT3TXXXdp0aJFkqSJEydq165d+uyzzyRdqHbr6uq0adMmLV26VJLU1NSkQCCg5uZmrVixIq7zUPECAIw0sEiCnZGIefPm6YMPPtDx48clSZ9//rn279+vO+64Q5J04sQJdXR0qLy8PHqMz+fT/PnzdeDAgbjPQ8ULADCSZXNZQOuXY7u7u2O2+3w++Xy+QfuvX79eXV1dmjx5srxer8LhsLZs2aL7779fktTR0SFJCgQCMccFAgGdPHky7rioeAEArlZUVKTc3NzoCAaDF91v9+7d2rlzp5qbm3Xo0CE1NTXpxRdfVFNTU8x+Hk/sLwOWZQ3adjlUvAAAIyVrPd5QKKScnJzo9otVu5L01FNPacOGDbrvvvskSTfeeKNOnjypYDCohx9+WPn5+ZIuVL4FBQXR4zo7OwdVwZeTtsR7fvx1GjUqM12nH8QaZd4bTqwMA2PypjuCwSJ87+Ji4vfOxOsUMfI6pTuC/wn3Ra68U5Ika3WinJycmMR7KWfPnlVGRuzF9nq90ceJiouLlZ+fr5aWFk2fPl2S1NfXp9bWVm3dujXuuKh4AQCQdOedd2rLli0aP368pk6dqsOHD6u2tlaPPvqopAst5qqqKtXU1KikpEQlJSWqqalRdna2li1bFvd5SLwAACOFbS4LmOixL7/8sp555hlVVFSos7NThYWFWrFihZ599tnoPuvWrdO5c+dUUVGh06dPa/bs2dq7d6/8fn/c5/FYlmUlFJlN3d3dys3N1bz5z9FqvgIT23AmtitpNcfHxO+dideJVvPlhfvO69D/e1pdXV1xtW+HYiBPPLH/LvmuHj3kz+n96Wf933lvpzTWoTDo2wkAgPvRagYAGCmiDEVs1Id2jk0lEi8AwEhhy6Owjbua7RybSmb+OgAAgEtR8QIAjJSs53hNQ+IFABjJsrk6kWXS7eC/QuIFABhpYEF7O8ebyMxfBwAAcCkqXgCAkSKWvXnaiKOvh4ofiRcAYKSIzTleO8emkplRAQDgUgkl3mAwqFmzZsnv9ysvL09LlizRsWPHUhUbAGAEi8hje5goocTb2tqqyspKHTx4UC0tLerv71d5ebnOnDmTqvgAACPUwJur7AwTJTTH+95778V8vWPHDuXl5am9vV233nprUgMDAMCNbN1c1dXVJUm67rrrLrlPb2+vent7o193d3fbOSUAYITg5qrfsCxLa9eu1bx581RaWnrJ/YLBoHJzc6OjqKhoqKcEAIwgEXmir40c0nDDHO+vrVq1Sl988YV27dp12f02btyorq6u6AiFQkM9JQAAw96QWs2rV6/Wnj17tG/fPo0bN+6y+/p8Pvl8viEFBwAYuSybdyZbhla8CSVey7K0evVqvfXWW/roo49UXFycqrgAACMcqxNJqqysVHNzs95++235/X51dHRIknJzc5WVlZWSAAEAIxM3V0lqaGhQV1eXFixYoIKCgujYvXt3quIDAMBVEm41AwDgBFrNAAA4yO5rH133OBEAAEgcFS8AwEi0mgEAcJBbEy+tZgAAHETFCwAwklsrXhIvAMBIbk28tJoBAHAQFS8AwEiW7D2La+orn0i8AAAjubXVTOIFABiJxJtkZwrHyDtmTLpOP0jEm+4IBjNxYQ2L6xQXy2veP3i+d/HhOl1euNegYIYpKl4AgJGoeAEAcJBbEy89AwAAHETFCwAwkmV5ZNmoWu0cm0okXgCAkViPFwAA2EbFCwAwkltvriLxAgCM5NY5XlrNAAA4iIoXAGAkWs0AADjIra1mEi8AwEiWzYrX1MTLHC8AAA6i4gUAGMmSZNlYzd7GoSlF4gUAGCkijzy8uQoAANhBxQsAMBJ3NQMA4KCI5ZHHhc/x0moGAMBBVLwAACNZls27mg29rZnECwAwklvneGk1AwDgICpeAICR3FrxkngBAEZy613NJF4AgJHcenMVc7wAADiIihcAYKQLFa+dOd4kBpNEJF4AgJHcenMVrWYAABxExQsAMJIle2vqGtpppuIFAJhpoNVsZyTq22+/1QMPPKCxY8cqOztb06ZNU3t7+69islRdXa3CwkJlZWVpwYIFOnr0aELnIPECACDp9OnTmjt3rkaPHq13331XX331lV566SVdc8010X22bdum2tpa1dfXq62tTfn5+Vq4cKF6enriPg+tZgCAmRzuNW/dulVFRUXasWNHdNvEiRP/93GWpbq6Om3atElLly6VJDU1NSkQCKi5uVkrVqyI6zxUvAAAM9ltM//Sau7u7o4Zvb29Fz3dnj17VFZWprvvvlt5eXmaPn26tm/fHv3zEydOqKOjQ+Xl5dFtPp9P8+fP14EDB+L+a5F4AQBGGnhzlZ0hSUVFRcrNzY2OYDB40fN98803amhoUElJid5//32tXLlSTzzxhF577TVJUkdHhyQpEAjEHBcIBKJ/Fg9azQAAVwuFQsrJyYl+7fP5LrpfJBJRWVmZampqJEnTp0/X0aNH1dDQoIceeii6n8cTe9OWZVmDtl1O2hLvmQKPvD5zHm62DKz9jYzJm+4IBjPzOpn3IIOR18nEmAz8GVeGOT9PkfMRx86VrBdo5OTkxCTeSykoKNANN9wQs23KlCl64403JEn5+fmSLlS+BQUF0X06OzsHVcGXY+CPPQAAujBHa3ckYO7cuTp27FjMtuPHj2vChAmSpOLiYuXn56ulpSX65319fWptbdWcOXPiPg+tZgAAJD355JOaM2eOampqdM899+jTTz9VY2OjGhsbJV1oMVdVVammpkYlJSUqKSlRTU2NsrOztWzZsrjPQ+IFABjJ6WUBZ82apbfeeksbN27U5s2bVVxcrLq6Oi1fvjy6z7p163Tu3DlVVFTo9OnTmj17tvbu3Su/3x/3eUi8AAAzpeGdkYsXL9bixYsv+ecej0fV1dWqrq4ecljM8QIA4CAqXgCAkdy6LCCJFwBgLnOepEoaWs0AADiIihcAYCRazQAAOCkNdzU7gcQLADCU55dh53jzMMcLAICDqHgBAGai1QwAgINcmnhttZqDwWD0pdEAAODKhlzxtrW1qbGxUTfddFMy4wEA4IIhLO036HgDDani/emnn7R8+XJt375d1157bbJjAgAgujqRnWGiISXeyspKLVq0SLfffvsV9+3t7VV3d3fMAABgpEq41fz666/r0KFDamtri2v/YDCo559/PuHAAAAjHDdXSaFQSGvWrNHOnTuVmZkZ1zEbN25UV1dXdIRCoSEFCgAYYQbmeO0MAyVU8ba3t6uzs1MzZ86MbguHw9q3b5/q6+vV29srr9cbc4zP55PP50tOtAAADHMJJd7bbrtNX375Zcy2Rx55RJMnT9b69esHJV0AAIbKY10Ydo43UUKJ1+/3q7S0NGbbVVddpbFjxw7aDgCALS6d4+XNVQAAM7n0OV7bifejjz5KQhgAAIwMVLwAADPRagYAwEEuTbysxwsAgIOoeAEAZnJpxUviBQCYyaV3NdNqBgDAQVS8AAAj8eYqAACc5NI5XlrNAAA4iMQLAICDaDUDAIzkkc053qRFklxpS7znAxFlZEbSdfpBrAwDJwMMXGXR8hp4nUy8g8LE62Rgf8vjNef/gAEeA/8v8Jj083T2vHPn4nEiAABgF61mAICZXHpXM4kXAGAmlyZeWs0AADiIihcAYCTeXAUAgJNoNQMAALuoeAEAZnJpxUviBQAYya1zvLSaAQBwEBUvAMBMLn1lJIkXAGAm5ngBAHAOc7wAAMA2Kl4AgJloNQMA4CCbrWZTEy+tZgAAHETFCwAwE61mAAAc5NLES6sZAAAHUfECAIzEc7wAAMA2Ei8AAA6i1QwAMJNLb64i8QIAjOTWOV4SLwDAXIYmTzuY4wUAwEFUvAAAMzHHCwCAc9w6x0urGQAAB1HxAgDM5NJWMxUvAMBIA61mO8OOYDAoj8ejqqqq6DbLslRdXa3CwkJlZWVpwYIFOnr0aEKfS+IFAOA32tra1NjYqJtuuilm+7Zt21RbW6v6+nq1tbUpPz9fCxcuVE9PT9yfTeIFAJjJSsIYgp9++knLly/X9u3bde211/4vHMtSXV2dNm3apKVLl6q0tFRNTU06e/asmpub4/58Ei8AwExJSrzd3d0xo7e397Knrays1KJFi3T77bfHbD9x4oQ6OjpUXl4e3ebz+TR//nwdOHAg7r8WiRcA4GpFRUXKzc2NjmAweMl9X3/9dR06dOii+3R0dEiSAoFAzPZAIBD9s3ik7a5mb95ZZWRH0nX6QTIyzLv9zes15/oMMPE6jcow7zqNMvB75zXxOhFTXEZ7w+kOIar/TK9OOnSuZD3HGwqFlJOTE93u8/kuun8oFNKaNWu0d+9eZWZmXvpzPZ6Yry3LGrTtcnicCABgpiQ9TpSTkxOTeC+lvb1dnZ2dmjlzZnRbOBzWvn37VF9fr2PHjkm6UPkWFBRE9+ns7BxUBV8OrWYAgJkcvrnqtttu05dffqkjR45ER1lZmZYvX64jR45o0qRJys/PV0tLS/SYvr4+tba2as6cOXGfh4oXAABJfr9fpaWlMduuuuoqjR07Nrq9qqpKNTU1KikpUUlJiWpqapSdna1ly5bFfR4SLwDASCa+q3ndunU6d+6cKioqdPr0ac2ePVt79+6V3++P+zNIvAAAMxnwysiPPvoo5muPx6Pq6mpVV1cP+TOZ4wUAwEFUvAAAI5nYak4GEi8AwEwGtJpTgVYzAAAOouIFAJjJpRUviRcAYCTPL8PO8Sai1QwAgIOoeAEAZqLVDACAc9z6OFHCreZvv/1WDzzwgMaOHavs7GxNmzZN7e3tqYgNADCSObxIglMSqnhPnz6tuXPn6ne/+53effdd5eXl6T//+Y+uueaaFIUHAIC7JJR4t27dqqKiIu3YsSO6beLEicmOCQCACwytWu1IqNW8Z88elZWV6e6771ZeXp6mT5+u7du3X/aY3t5edXd3xwwAAK5kYI7XzjBRQon3m2++UUNDg0pKSvT+++9r5cqVeuKJJ/Taa69d8phgMKjc3NzoKCoqsh00AADDVUKJNxKJaMaMGaqpqdH06dO1YsUK/elPf1JDQ8Mlj9m4caO6urqiIxQK2Q4aADACcHOVVFBQoBtuuCFm25QpU/TGG29c8hifzyefzze06AAAIxaPE0maO3eujh07FrPt+PHjmjBhQlKDAgDArRJKvE8++aQOHjyompoa/fvf/1Zzc7MaGxtVWVmZqvgAACOVS1vNCSXeWbNm6a233tKuXbtUWlqqv/zlL6qrq9Py5ctTFR8AYIRy613NCb8ycvHixVq8eHEqYgEAwPV4VzMAwEwskgAAgINIvAAAOIfHiQAAgG1UvAAAM9FqBgDAOR7Lkscaeva0c2wq0WoGAMBBVLwAADPRagYAwDnc1QwAAGyj4gUAmIlWc3JNyvtRo64yZ53eURmRdIcwyJiM/nSHMAjXKT6jPeZdJ5/XvOtk5vcunO4QBvEZdJ16fT9rn0PnotUMAABso9UMADATrWYAAJzj1lYziRcAYCaXVrzM8QIA4CAqXgCAsUxtF9tB4gUAmMmyLgw7xxuIVjMAAA6i4gUAGIm7mgEAcBJ3NQMAALuoeAEARvJELgw7x5uIxAsAMBOtZgAAYBcVLwDASNzVDACAk1z6Ag0SLwDASG6teJnjBQDAQVS8AAAzufSuZhIvAMBItJoBAIBtVLwAADNxVzMAAM6h1QwAAGyj4gUAmIm7mgEAcA6tZgAAYBsVLwDATBHrwrBzvIFIvAAAMzHHCwCAczyyOcebtEiSizleAAAcROIFAJhp4M1VdkYCgsGgZs2aJb/fr7y8PC1ZskTHjh37TUiWqqurVVhYqKysLC1YsEBHjx5N6DwkXgCAkQYeJ7IzEtHa2qrKykodPHhQLS0t6u/vV3l5uc6cORPdZ9u2baqtrVV9fb3a2tqUn5+vhQsXqqenJ+7zMMcLAICk9957L+brHTt2KC8vT+3t7br11ltlWZbq6uq0adMmLV26VJLU1NSkQCCg5uZmrVixIq7zUPECAMxkJWHY0NXVJUm67rrrJEknTpxQR0eHysvLo/v4fD7Nnz9fBw4ciPtzqXgBAEbyWJY8NlYYGji2u7s7ZrvP55PP57vssZZlae3atZo3b55KS0slSR0dHZKkQCAQs28gENDJkyfjjittibfs2lPyXT06XacfJDPj53SHMIjPwJhGe8LpDmGQTI+J16k/3SEMYuLP+Bh+nuJiUkxnFFZtuoNIUFFRUczXzz33nKqrqy97zKpVq/TFF19o//79g/7M44l9UMmyrEHbLoeKFwBgpsgvw87xkkKhkHJycqKbr1Ttrl69Wnv27NG+ffs0bty46Pb8/HxJFyrfgoKC6PbOzs5BVfDlMMcLADDSQKvZzpCknJycmHGpxGtZllatWqU333xTH374oYqLi2P+vLi4WPn5+WppaYlu6+vrU2trq+bMmRP334uKFwAASZWVlWpubtbbb78tv98fndPNzc1VVlaWPB6PqqqqVFNTo5KSEpWUlKimpkbZ2dlatmxZ3Och8QIAzOTwu5obGhokSQsWLIjZvmPHDv3xj3+UJK1bt07nzp1TRUWFTp8+rdmzZ2vv3r3y+/1xn4fECwAw0xDePjXo+IR2v/L+Ho9H1dXVV7w563JIvAAAIw3l7VO/Pd5E3FwFAICDqHgBAGZyuNXsFBIvAMBInsiFYed4E9FqBgDAQVS8AAAz0WoGAMBBDj/H6xRazQAAOIiKFwBgpGQtC2iahCre/v5+Pf300youLlZWVpYmTZqkzZs3KxIx9NYxAMDwNTDHa2cYKKGKd+vWrXr11VfV1NSkqVOn6rPPPtMjjzyi3NxcrVmzJlUxAgDgGgkl3k8++UR33XWXFi1aJEmaOHGidu3apc8++ywlwQEARjBL9tbjNbPgTazVPG/ePH3wwQc6fvy4JOnzzz/X/v37dccdd1zymN7eXnV3d8cMAACuJFnr8ZomoYp3/fr16urq0uTJk+X1ehUOh7Vlyxbdf//9lzwmGAzq+eeftx0oAGCEsWTzOd6kRZJUCVW8u3fv1s6dO9Xc3KxDhw6pqalJL774opqami55zMaNG9XV1RUdoVDIdtAAAAxXCVW8Tz31lDZs2KD77rtPknTjjTfq5MmTCgaDevjhhy96jM/nk8/nsx8pAGBk4c1V0tmzZ5WREVske71eHicCACRfRJLH5vEGSijx3nnnndqyZYvGjx+vqVOn6vDhw6qtrdWjjz6aqvgAAHCVhBLvyy+/rGeeeUYVFRXq7OxUYWGhVqxYoWeffTZV8QEARii3vrkqocTr9/tVV1enurq6FIUDAMAvXDrHyyIJAAA4iEUSAABmcmnFS+IFAJjJpYmXVjMAAA6i4gUAmInneAEAcA6PEwEA4CTmeAEAgF1UvAAAM0UsyWOjao2YWfGSeAEAZqLVDAAA7Epbxbvg6q91ld+cvJ/p6U93CIOM9ph3L3ymJ5zuEAbJtNOKSpHR6Q7gIjI95vx7GzDawJh8HvO+e6M93nSHENXd7+S/N5sVr8z7v0Gi1QwAMBWtZgAAYBcVLwDATBFLttrF3NUMAEACrMiFYed4A9FqBgDAQVS8AAAzufTmKhIvAMBMzPECAOAgl1a8zPECAOAgKl4AgJks2ax4kxZJUpF4AQBmotUMAADsouIFAJgpEpFk4yUYETNfoEHiBQCYiVYzAACwi4oXAGAml1a8JF4AgJlc+uYqWs0AADiIihcAYCTLisiysbSfnWNTicQLADCTZdlrFzPHCwBAAiybc7yGJl7meAEAcBAVLwDATJGI5LExT8scLwAACaDVDAAA7KLiBQAYyYpEZNloNfM4EQAAiaDVDAAA7KLiBQCYKWJJHvdVvCReAICZLEuSnceJzEy8tJoBAHAQFS8AwEhWxJJlo9VsUfECAJAAK2J/DMErr7yi4uJiZWZmaubMmfr444+T+tci8QIAjGRFLNsjUbt371ZVVZU2bdqkw4cP65ZbbtHvf/97nTp1Kml/LxIvAAC/qK2t1WOPPabHH39cU6ZMUV1dnYqKitTQ0JC0czg+xzvQcz/zk1lvFOm38yLuFBltYEwmXqef7TxukCKj0x3ARfzsSXcEg40yMCafgT/joz3mXKjuX/7vdmL+tN/qtbXQQb9+liR1d3fHbPf5fPL5fIP27+vrU3t7uzZs2BCzvby8XAcOHBhyHL/leOLt6emRJP2fOckr2wEAzurp6VFubm5KPnvMmDHKz8/X/o53bH/W1VdfraKiophtzz33nKqrqwft+8MPPygcDisQCMRsDwQC6ujosB3LAMcTb2FhoUKhkPx+vzw2fovr7u5WUVGRQqGQcnJykhihu3Cd4sN1ig/XKT5uvk6WZamnp0eFhYUpO0dmZqZOnDihvr4+259lWdagXHOxavfXfrv/xT7DDscTb0ZGhsaNG5e0z8vJyXHdD3YqcJ3iw3WKD9cpPm69TqmqdH8tMzNTmZmZKT/Pr11//fXyer2DqtvOzs5BVbAd3FwFAIAutLhnzpyplpaWmO0tLS2aM2dO0s7DCzQAAPjF2rVr9eCDD6qsrEw333yzGhsbderUKa1cuTJp5xi2idfn8+m55567Yq9+pOM6xYfrFB+uU3y4TsPXvffeqx9//FGbN2/W999/r9LSUr3zzjuaMGFC0s7hsUx9pxYAAC7EHC8AAA4i8QIA4CASLwAADiLxAgDgoGGbeFO9bNNwFwwGNWvWLPn9fuXl5WnJkiU6duxYusMyWjAYlMfjUVVVVbpDMc63336rBx54QGPHjlV2dramTZum9vb2dIdllP7+fj399NMqLi5WVlaWJk2apM2bNysSMe/dz0ivYZl4nVi2abhrbW1VZWWlDh48qJaWFvX396u8vFxnzpxJd2hGamtrU2Njo2666aZ0h2Kc06dPa+7cuRo9erTeffddffXVV3rppZd0zTXXpDs0o2zdulWvvvqq6uvr9fXXX2vbtm164YUX9PLLL6c7NBhmWD5ONHv2bM2YMSNmmaYpU6ZoyZIlCgaDaYzMXP/973+Vl5en1tZW3XrrrekOxyg//fSTZsyYoVdeeUV//etfNW3aNNXV1aU7LGNs2LBB//rXv+gqXcHixYsVCAT097//PbrtD3/4g7Kzs/WPf/wjjZHBNMOu4h1Ytqm8vDxme7KXbXKbrq4uSdJ1112X5kjMU1lZqUWLFun2229PdyhG2rNnj8rKynT33XcrLy9P06dP1/bt29MdlnHmzZunDz74QMePH5ckff7559q/f7/uuOOONEcG0wy7N1c5tWyTm1iWpbVr12revHkqLS1NdzhGef3113Xo0CG1tbWlOxRjffPNN2poaNDatWv15z//WZ9++qmeeOIJ+Xw+PfTQQ+kOzxjr169XV1eXJk+eLK/Xq3A4rC1btuj+++9Pd2gwzLBLvANSvWyTm6xatUpffPGF9u/fn+5QjBIKhbRmzRrt3bvX8VVQhpNIJKKysjLV1NRIkqZPn66jR4+qoaGBxPsru3fv1s6dO9Xc3KypU6fqyJEjqqqqUmFhoR5++OF0hweDDLvE69SyTW6xevVq7dmzR/v27Uvqcoxu0N7ers7OTs2cOTO6LRwOa9++faqvr1dvb6+8Xm8aIzRDQUGBbrjhhphtU6ZM0RtvvJGmiMz01FNPacOGDbrvvvskSTfeeKNOnjypYDBI4kWMYTfH69SyTcOdZVlatWqV3nzzTX344YcqLi5Od0jGue222/Tll1/qyJEj0VFWVqbly5fryJEjJN1fzJ07d9CjaMePH0/qS+Pd4OzZs8rIiP0v1ev18jgRBhl2Fa/kzLJNw11lZaWam5v19ttvy+/3RzsEubm5ysrKSnN0ZvD7/YPmvK+66iqNHTuWufBfefLJJzVnzhzV1NTonnvu0aeffqrGxkY1NjamOzSj3HnnndqyZYvGjx+vqVOn6vDhw6qtrdWjjz6a7tBgGmuY+tvf/mZNmDDBGjNmjDVjxgyrtbU13SEZRdJFx44dO9IdmtHmz59vrVmzJt1hGOef//ynVVpaavl8Pmvy5MlWY2NjukMyTnd3t7VmzRpr/PjxVmZmpjVp0iRr06ZNVm9vb7pDg2GG5XO8AAAMV8NujhcAgOGMxAsAgINIvAAAOIjECwCAg0i8AAA4iMQLAICDSLwAADiIxAsAgINIvAAAOIjECwCAg0i8AAA4iMQLAICD/j9knPS7yMOqgAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(encoded_values)\n",
    "plt.colorbar() ;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a0b4dc45-4a5f-415e-9f40-4acc419a26f0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "reader = Reader(raster_path_cea)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d0ec69fe-cf96-4606-9ad3-80ac459ba1d5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "west, north = transform * (0, 0)\n",
    "east, south = transform * (height, height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8d8917ad-e33a-4f4a-9aaf-5466d6c10528",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# shift polygon to the east\n",
    "#east += transform.a * 10.5\n",
    "#west += transform.a * 10.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f0b4dece-3254-4c08-b8ec-395618b9be27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "geojson = {\n",
    "    \"type\": \"Feature\",\n",
    "    \"properties\": {},\n",
    "    \"geometry\": {\n",
    "        \"type\": \"Polygon\",\n",
    "        \"coordinates\": [\n",
    "            [\n",
    "                [west, north],\n",
    "                [east, south],\n",
    "                [west, south],\n",
    "                [west, north]\n",
    "            ]\n",
    "        ]\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b7ab9680-b79b-4b0c-8c93-0d410b959174",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#geom_wgs = rasterio.warp.transform_geom(crs, \"epsg:4326\", geojson[\"geometry\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e64442c2-634e-4ee5-ad40-c90f664775ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "feature = Feature(**geojson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9d9f2810-6cc3-4011-8022-50fdc986e87b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, -2000000.0, 2000000.0, 0.0)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reader.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "49b8848a-0614-4dd7-ad3d-166aa8cfd2b9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'type': 'Feature',\n",
       " 'properties': {},\n",
       " 'geometry': {'type': 'Polygon',\n",
       "  'coordinates': [[[0.0, 0.0],\n",
       "    [2000000.0, -2000000.0],\n",
       "    [0.0, -2000000.0],\n",
       "    [0.0, 0.0]]]}}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geojson"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d5f2605-3573-43f1-a8c7-f429fccf2078",
   "metadata": {},
   "source": [
    "# Calculate coverage array and stats with Reader method\n",
    "\n",
    "Note that the choice of `cover_scale` matters a lot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f2a63e7c-716b-40cb-bc1a-aa4883ebeba6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "coverage, stats = next(reader.zonal_statistics(\n",
    "    feature,\n",
    "    align_bounds_with_dataset=False,\n",
    "    shape_crs=reader.crs,\n",
    "    cover_scale=1000\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "a68a68a4-3736-419c-a246-05f5bac75f50",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5005, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 0.5005, 0.    , 0.    , 0.    , 0.    , 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 0.5005, 0.    , 0.    , 0.    , 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 0.5005, 0.    , 0.    , 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 0.5005, 0.    , 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 1.    , 0.5005, 0.    , 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 0.5005, 0.    ,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 0.5005,\n",
       "        0.    , 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    ,\n",
       "        0.5005, 0.    ],\n",
       "       [1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    ,\n",
       "        1.    , 0.5005]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "06c61894-70dd-49e2-aace-058d7a56fd34",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 10)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coverage.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d3795a87-37f2-4e6d-b0b0-a44c59580b9b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.    , 0.5005, 1.    ])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.unique(coverage)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "09093156-994e-4f85-9203-be2ba781f296",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fa2fef3b250>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAGiCAYAAADUc67xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAllUlEQVR4nO3de3BU9f3/8dcmyC5qEis0kZQQA9ohEi+4sW2AWK9xUKn2ayXewAs4ZsItSWVIxAqiktEqTSsmGAXxBma8UHAmChkdbiIjxKCOMtKiJVtMmoY6CWhNzO75/QHsz+0GyeYke87mPB8z5w8+nE8+72Y6vnm/P59zjsswDEMAAMAycVYHAACA05GMAQCwGMkYAACLkYwBALAYyRgAAIuRjAEAsBjJGAAAi5GMAQCwGMkYAACLkYwBALAYyRgAgKO2bNmiyZMnKzU1VS6XS3/9619POGfz5s3yer3yeDwaNWqUli9fHvG6JGMAAI765ptvdP7552vZsmU9uv/LL7/U1VdfrdzcXDU0NOi+++7TnDlz9Prrr0e0rosPRQAAEM7lcmnt2rW6/vrrj3vP/PnztX79eu3Zsyc4VlBQoI8++kjvv/9+j9caZCbQ3ggEAvrqq6+UkJAgl8sV7eUBACYYhqFDhw4pNTVVcXH911z97rvv1NnZafrnGIYRlmvcbrfcbrfpny1J77//vvLy8kLGrrrqKq1YsULff/+9TjrppB79nKgn46+++kppaWnRXhYA0Id8Pp9GjBjRLz/7u+++U0b6qWpu8Zv+WaeeeqoOHz4cMrZw4UItWrTI9M+WpObmZqWkpISMpaSkqKurS62trRo+fHiPfk7Uk3FCQoIkacSi+xXn8UR7+eMaVbrT6hAAwPa69L22qTb43/L+0NnZqeYWv76sT1diQu+r7/ZDAWV498vn8ykxMTE43ldV8TH/W3kf2/2NpPsb9WR8LLg4j8dWyXiQq2etBABwtKOnjKKxzZiYEGcqGQd/TmJiSDLuS2eccYaam5tDxlpaWjRo0CANHTq0xz8n6skYAICe8BsB+U0cMfYbgb4L5jhycnL05ptvhoxt3LhR2dnZPd4vlni0CQBgUwEZpq9IHT58WLt379bu3bslHXl0affu3WpsbJQklZWVadq0acH7CwoKtH//fpWUlGjPnj1auXKlVqxYoXvvvTeidamMAQC2FFBAZmrb3szetWuXLr300uCfS0pKJEm33367Vq1apaampmBilqSMjAzV1taquLhYTz31lFJTU/WXv/xFN9xwQ0TrkowBADjqkksu0Y+9fmPVqlVhY7/+9a/14YcfmlqXZAwAsCW/Ychv4r1UZuZGG8kYAGBLvd33/eH8WMEBLgAALEZlDACwpYAM+R1SGZOMAQC2RJsaAABEDZUxAMCWnHSauleVcWVlpTIyMuTxeOT1erV169a+jgsA4HCBPrhiRcTJuKamRkVFRVqwYIEaGhqUm5urSZMmhbyRBAAA9FzEyXjp0qWaPn26ZsyYoczMTFVUVCgtLU1VVVX9ER8AwKH8R09Tm7liRUR7xp2dnaqvr1dpaWnIeF5enrZv397tnI6ODnV0dAT/3N7e3oswAQBO4zdk8qtNfRdLf4uoMm5tbZXf71dKSkrIeEpKStj3HI8pLy9XUlJS8EpLS+t9tAAAx2DP+AT+96PShmEc90PTZWVlamtrC14+n683SwIAMGBF1KYeNmyY4uPjw6rglpaWsGr5GLfbLbfb3fsIAQCOFJBLfnVf6PV0fqyIqDIePHiwvF6v6urqQsbr6uo0fvz4Pg0MAOBsAcP8FSsifulHSUmJpk6dquzsbOXk5Ki6ulqNjY0qKCjoj/gAABjwIk7G+fn5OnjwoBYvXqympiZlZWWptrZW6enp/REfAMCh/Cbb1GbmRluvXodZWFiowsLCvo4FAIAgJyVjPhQBAIDF+FAEAMCWAoZLAcPEaWoTc6ONZAwAsCXa1AAAIGqojAEAtuRXnPwmakZ/H8bS30jGAABbMkzuGRvsGQMAYA57xgAAIGqojAEAtuQ34uQ3TOwZD+R3UwMAEA0BuRQw0cANKHayMW1qAAAsZlll/NH/rVRign3+LTBa9vvq1FnFO6wOAQAs46QDXLSpAQC2ZH7PmDY1AADoISpjAIAtHTnAZeJDEbSpAQAwJ2DydZicpgYAAD1GZQwAsCUnHeAiGQMAbCmgOMe89INkDACwJb/hkt/El5fMzI029owBALAYlTEAwJb8Jk9T+2lTAwBgTsCIU8DEAa5ADB3gok0NAIDFqIwBALZEmxoAAIsFZO5EdKDvQul3tKkBALAYlTEAwJbMv/QjdupNkjEAwJbMvw4zdpJx7EQKAMAARWUMALAlvmcMAIDFnNSmJhkDAGzJ/HPGsZOMYydSAAAGKCpjAIAtBQyXAmZe+hFDn1AkGQMAbClgsk0dS88Zx06kAAAMUFTGAABbMv8JxdipN0nGAABb8sslv4lnhc3MjbbY+WcDAAADFJUxAMCWaFMDAGAxv8y1mv19F0q/i51/NgAAMEBRGQMAbIk2NQAAFnPShyJiJ1IAgKMYRz+h2NvL6OV+c2VlpTIyMuTxeOT1erV169Yfvf/ll1/W+eefr5NPPlnDhw/XnXfeqYMHD0a0JskYAICjampqVFRUpAULFqihoUG5ubmaNGmSGhsbu71/27ZtmjZtmqZPn65PP/1Ur776qnbu3KkZM2ZEtC7JGABgS8fa1GauSC1dulTTp0/XjBkzlJmZqYqKCqWlpamqqqrb+3fs2KEzzzxTc+bMUUZGhiZOnKh77rlHu3btimhd9oyP2pe/3OoQwoxWgdUhhDmreIfVIQBwiL76alN7e3vIuNvtltvtDru/s7NT9fX1Ki0tDRnPy8vT9u3bu11j/PjxWrBggWprazVp0iS1tLTotdde0zXXXBNRrFTGAIABLS0tTUlJScGrvLy82/taW1vl9/uVkpISMp6SkqLm5uZu54wfP14vv/yy8vPzNXjwYJ1xxhk67bTT9OSTT0YUI5UxAMCW/CY/oXhsrs/nU2JiYnC8u6r4h1yu0GrcMIywsWM+++wzzZkzRw888ICuuuoqNTU1ad68eSooKNCKFSt6HCvJGABgS33Vpk5MTAxJxsczbNgwxcfHh1XBLS0tYdXyMeXl5ZowYYLmzZsnSTrvvPN0yimnKDc3Vw8//LCGDx/eo1hpUwMAIGnw4MHyer2qq6sLGa+rq9P48eO7nfPtt98qLi40lcbHx0s6UlH3FJUxAMCWAopTwETN2Ju5JSUlmjp1qrKzs5WTk6Pq6mo1NjaqoODIgdqysjIdOHBAL7zwgiRp8uTJuvvuu1VVVRVsUxcVFekXv/iFUlNTe7wuyRgAYEt+wyW/iTZ1b+bm5+fr4MGDWrx4sZqampSVlaXa2lqlp6dLkpqamkKeOb7jjjt06NAhLVu2TL///e912mmn6bLLLtOjjz4a0bokYwAAfqCwsFCFhYXd/t2qVavCxmbPnq3Zs2ebWpNkDACwpb46wBULSMYAAFsyTH61yYihD0WQjAEAtuSXS/5efuzh2PxYETv/bAAAYICiMgYA2FLAMLfvG+j5Y76WIxkDAGwpYHLP2MzcaIudSAEAGKAiSsbl5eW66KKLlJCQoOTkZF1//fX6/PPP+ys2AICDBeQyfcWKiJLx5s2bNXPmTO3YsUN1dXXq6upSXl6evvnmm/6KDwDgUMfewGXmihUR7Rm//fbbIX9+7rnnlJycrPr6el188cV9GhgAAE5h6gBXW1ubJOn0008/7j0dHR3q6OgI/rm9vd3MkgAAh+AAVw8YhqGSkhJNnDhRWVlZx72vvLxcSUlJwSstLa23SwIAHCQgV/CVmL26Buqe8Q/NmjVLH3/8sdasWfOj95WVlamtrS14+Xy+3i4JAMCA1Ks29ezZs7V+/Xpt2bJFI0aM+NF73W633G53r4IDADiXYfJEtBFDlXFEydgwDM2ePVtr167Vpk2blJGR0V9xAQAcjq82HcfMmTO1evVqrVu3TgkJCWpubpYkJSUlaciQIf0SIADAmTjAdRxVVVVqa2vTJZdcouHDhwevmpqa/ooPAIABL+I2NQAA0UCbGgAAi5l9paUjHm0CAAB9g8oYAGBLtKkBALCYk5IxbWoAACxGZQwAsCUnVcYkYwCALTkpGdOmBgDAYlTGAABbMmTuWeFYek0VyRgAYEtOalOTjAEAtkQyhi3sy19udQhhRqvA6hDCnFW8w+oQAMAUkjEAwJaojAEAsJiTkjGPNgEAYDEqYwCALRmGS4aJ6tbM3GgjGQMAbInvGQMAgKihMgYA2JKTDnCRjAEAtuSkPWPa1AAAWIzKGABgS7SpAQCwmJPa1CRjAIAtGSYr41hKxuwZAwBgMSpjAIAtGZIMw9z8WEEyBgDYUkAuuXgDFwAAiAYqYwCALXGaGgAAiwUMl1wOec6YNjUAABajMgYA2JJhmDxNHUPHqUnGAABbctKeMW1qAAAsRmUMALAlJ1XGJGMAgC056TQ1yRgAYEtOOsDFnjEAABajMgYA2NKRytjMnnEfBtPPSMYAAFty0gEu2tQAAFiMyhgAYEuGzH2TOIa61CRjAIA90aYGAABRQ2UMALAnB/WpqYwBAPZ0tE3d20u9bFNXVlYqIyNDHo9HXq9XW7du/dH7Ozo6tGDBAqWnp8vtdmv06NFauXJlRGtSGQMAbMmKN3DV1NSoqKhIlZWVmjBhgp5++mlNmjRJn332mUaOHNntnClTpuhf//qXVqxYobPOOkstLS3q6uqKaF2SMQAARy1dulTTp0/XjBkzJEkVFRXasGGDqqqqVF5eHnb/22+/rc2bN+uLL77Q6aefLkk688wzI16XZIyI7MtfbnUIYUarwOoQwpxVvMPqEICY11enqdvb20PG3W633G532P2dnZ2qr69XaWlpyHheXp62b9/e7Rrr169Xdna2HnvsMb344os65ZRT9Jvf/EYPPfSQhgwZ0uNYScYAAHsyse8bnC8pLS0tZHjhwoVatGhR2O2tra3y+/1KSUkJGU9JSVFzc3O3S3zxxRfatm2bPB6P1q5dq9bWVhUWFuo///lPRPvGJGMAwIDm8/mUmJgY/HN3VfEPuVyh/wAwDCNs7JhAICCXy6WXX35ZSUlJko60un/3u9/pqaee6nF1TDIGANhSXx3gSkxMDEnGxzNs2DDFx8eHVcEtLS1h1fIxw4cP189+9rNgIpakzMxMGYahf/7znzr77LN7FCuPNgEA7MnogysCgwcPltfrVV1dXch4XV2dxo8f3+2cCRMm6KuvvtLhw4eDY3v37lVcXJxGjBjR47VJxgAAHFVSUqJnn31WK1eu1J49e1RcXKzGxkYVFBw5KFpWVqZp06YF77/llls0dOhQ3Xnnnfrss8+0ZcsWzZs3T3fddRcHuAAAsc+Kd1Pn5+fr4MGDWrx4sZqampSVlaXa2lqlp6dLkpqamtTY2Bi8/9RTT1VdXZ1mz56t7OxsDR06VFOmTNHDDz8c0bokYwCAfVnwSsvCwkIVFhZ2+3erVq0KGxszZkxYaztStKkBALAYlTEAwJac9AlFkjEAwJ4c9NUmkjEAwKZcRy8z82MDe8YAAFiMyhgAYE+0qQEAsJiDkrGpNnV5eblcLpeKior6KBwAAJyn15Xxzp07VV1drfPOO68v4wEA4Ig++oRiLOhVZXz48GHdeuuteuaZZ/STn/ykr2MCACD41SYzV6zoVTKeOXOmrrnmGl1xxRUnvLejo0Pt7e0hFwAA+P8iblO/8sor+vDDD7Vz584e3V9eXq4HH3ww4sAAAA7HAa7u+Xw+zZ07Vy+99JI8Hk+P5pSVlamtrS14+Xy+XgUKAHCYY3vGZq4YEVFlXF9fr5aWFnm93uCY3+/Xli1btGzZMnV0dCg+Pj5kjtvtltvt7ptoAQAYgCJKxpdffrk++eSTkLE777xTY8aM0fz588MSMQAAveUyjlxm5seKiJJxQkKCsrKyQsZOOeUUDR06NGwcAABTHLRnzBu4AAD25KDnjE0n402bNvVBGAAAOBeVMQDAnmhTAwBgMQclY75nDACAxaiMAQD25KDKmGQMALAnB52mpk0NAIDFqIwBALbEG7gAALCag/aMaVMDAGAxkjEAABajTQ0AsCWXTO4Z91kk/Y9kjJi3L3+51SGEGa0Cq0MIc1bxDqtDACLDo00AACBaqIwBAPbkoNPUJGMAgD05KBnTpgYAwGJUxgAAW+INXAAAWI02NQAAiBYqYwCAPTmoMiYZAwBsyUl7xrSpAQCwGJUxAMCeHPQ6TJIxAMCe2DMGAMBa7BkDAICooTIGANgTbWoAACxmsk0dS8mYNjUAABajMgYA2BNtagAALOagZEybGgAAi1EZAwBsieeMAQBA1JCMAQCwGG1qAIA9OegAF8kYAGBLTtozJhkDAOwrhhKqGewZAwBgMSpjAIA9sWcMAIC1nLRnTJsaAACLURkDAOyJNjUAANaiTQ0AAKKGZAwAsCejD65eqKysVEZGhjwej7xer7Zu3dqjee+9954GDRqkCy64IOI1ScYAAHuyIBnX1NSoqKhICxYsUENDg3JzczVp0iQ1Njb+6Ly2tjZNmzZNl19+eeSLimQMABjg2tvbQ66Ojo7j3rt06VJNnz5dM2bMUGZmpioqKpSWlqaqqqofXeOee+7RLbfcopycnF7FyAEuoB/sy19udQhhRqvA6hDCnFW8w+oQYGN9dYArLS0tZHzhwoVatGhR2P2dnZ2qr69XaWlpyHheXp62b99+3HWee+457du3Ty+99JIefvjhXsVKMgYA2FMfPdrk8/mUmJgYHHa73d3e3traKr/fr5SUlJDxlJQUNTc3dzvnb3/7m0pLS7V161YNGtT7lEoyBgDYUx8l48TExJBkfCIulyv0xxhG2Jgk+f1+3XLLLXrwwQf185//3ESgJGMAACRJw4YNU3x8fFgV3NLSElYtS9KhQ4e0a9cuNTQ0aNasWZKkQCAgwzA0aNAgbdy4UZdddlmP1iYZAwBsKdov/Rg8eLC8Xq/q6ur029/+NjheV1en6667Luz+xMREffLJJyFjlZWVevfdd/Xaa68pIyOjx2uTjAEA9mTB6zBLSko0depUZWdnKycnR9XV1WpsbFRBwZEDkGVlZTpw4IBeeOEFxcXFKSsrK2R+cnKyPB5P2PiJkIwBADgqPz9fBw8e1OLFi9XU1KSsrCzV1tYqPT1dktTU1HTCZ457w2UYRlTf3tne3q6kpCR9vXeUEhN4zBmIltE1PNoE87qM77VJ69TW1hbRoahIHMsTmbOWKN7t6fXP8Xd8pz3L7uvXWPsKlTEAwJ4c9NUmSlMAACxGZQwAsCcHVcYkYwCALbmOXmbmxwra1AAAWIzKGABgT7SpAQCwVrTfwGWliNvUBw4c0G233aahQ4fq5JNP1gUXXKD6+vr+iA0A4GRGH1wxIqLK+Ouvv9aECRN06aWX6q233lJycrL27dun0047rZ/CAwBg4IsoGT/66KNKS0vTc889Fxw788wz+zomAACOiKHq1oyI2tTr169Xdna2brzxRiUnJ2vcuHF65plnfnROR0eH2tvbQy4AAE7k2J6xmStWRJSMv/jiC1VVVenss8/Whg0bVFBQoDlz5uiFF1447pzy8nIlJSUFr7S0NNNBAwAwkESUjAOBgC688EItWbJE48aN0z333KO7775bVVVVx51TVlamtra24OXz+UwHDQBwAA5wdW/48OE655xzQsYyMzP1+uuvH3eO2+2W2+3uXXQAAMfi0abjmDBhgj7//POQsb179wa/8wgAACIXUTIuLi7Wjh07tGTJEv3973/X6tWrVV1drZkzZ/ZXfAAAp3JQmzqiZHzRRRdp7dq1WrNmjbKysvTQQw+poqJCt956a3/FBwBwKCedpo74dZjXXnutrr322v6IBQAAR+Ld1AAAe+JDEQAAWIxkDACAtXi0CQAARA2VMQDAnmhTAwBgLZdhyGX0PqOamRtttKkBALAYlTEAwJ5oUwMAYC1OUwMAgKihMgYA2BNtagADzb785VaHEGa0CqwOIcxZxTusDgFH0aYGAABRQ2UMALAn2tQAAFjLSW1qkjEAwJ4cVBmzZwwAgMWojAEAthVLrWYzSMYAAHsyjCOXmfkxgjY1AAAWozIGANgSp6kBALAap6kBAEC0UBkDAGzJFThymZkfK0jGAAB7ok0NAACihcoYAGBLnKYGAMBqDnrpB8kYAGBLTqqM2TMGAMBiVMYAAHty0GlqkjEAwJZoUwMAgKihMgYA2BOnqQEAsBZtagAAEDVUxgAAe+I0NQAA1qJNDQAAoobKGABgTwHjyGVmfowgGQMA7Ik9YwAArOWSyT3jPouk/7FnDACAxaiMAQD2xBu4AACwFo82AQDgUJWVlcrIyJDH45HX69XWrVuPe+8bb7yhK6+8Uj/96U+VmJionJwcbdiwIeI1ScYAAHsy+uCKUE1NjYqKirRgwQI1NDQoNzdXkyZNUmNjY7f3b9myRVdeeaVqa2tVX1+vSy+9VJMnT1ZDQ0NE69KmBgDYkssw5DKx73tsbnt7e8i42+2W2+3uds7SpUs1ffp0zZgxQ5JUUVGhDRs2qKqqSuXl5WH3V1RUhPx5yZIlWrdund58802NGzeux7GSjAFYZl/+cqtDCDNaBVaHEOas4h1WhxDT0tLSQv68cOFCLVq0KOy+zs5O1dfXq7S0NGQ8Ly9P27dv79FagUBAhw4d0umnnx5RjCRjAIA9BY5eZuZL8vl8SkxMDA4frypubW2V3+9XSkpKyHhKSoqam5t7tOQTTzyhb775RlOmTIkoVJIxAMCW+qpNnZiYGJKMTzjPFfq6EMMwwsa6s2bNGi1atEjr1q1TcnJyRLGSjAEAkDRs2DDFx8eHVcEtLS1h1fL/qqmp0fTp0/Xqq6/qiiuuiHhtTlMDAOwpyqepBw8eLK/Xq7q6upDxuro6jR8//rjz1qxZozvuuEOrV6/WNddcE9miR1EZAwDsyYI3cJWUlGjq1KnKzs5WTk6Oqqur1djYqIKCIwf7ysrKdODAAb3wwguSjiTiadOm6c9//rN+9atfBavqIUOGKCkpqcfrkowBALZkxRu48vPzdfDgQS1evFhNTU3KyspSbW2t0tPTJUlNTU0hzxw//fTT6urq0syZMzVz5szg+O23365Vq1b1eF2SMQAAP1BYWKjCwsJu/+5/E+ymTZv6ZE2SMQDAnvhQBAAA1nIFjlxm5scKTlMDAGAxKmMAgD3RpgYAwGK9/PJSyPwYQZsaAACLURkDAGypr95NHQsiqoy7urp0//33KyMjQ0OGDNGoUaO0ePFiBQIxdGQNABAbju0Zm7liRESV8aOPPqrly5fr+eef19ixY7Vr1y7deeedSkpK0ty5c/srRgAABrSIkvH777+v6667Lvgi7DPPPFNr1qzRrl27+iU4AICDGTL3PePYKYwja1NPnDhR77zzjvbu3StJ+uijj7Rt2zZdffXVx53T0dGh9vb2kAsAgBM5tmds5ooVEVXG8+fPV1tbm8aMGaP4+Hj5/X498sgjuvnmm487p7y8XA8++KDpQAEADmPI5HPGfRZJv4uoMq6pqdFLL72k1atX68MPP9Tzzz+vxx9/XM8///xx55SVlamtrS14+Xw+00EDADCQRFQZz5s3T6WlpbrpppskSeeee67279+v8vJy3X777d3Ocbvdcrvd5iMFADgLb+Dq3rfffqu4uNBiOj4+nkebAAB9LyDJZXJ+jIgoGU+ePFmPPPKIRo4cqbFjx6qhoUFLly7VXXfd1V/xAQAw4EWUjJ988kn94Q9/UGFhoVpaWpSamqp77rlHDzzwQH/FBwBwKCe9gSuiZJyQkKCKigpVVFT0UzgAABzloD1jPhQBAIDF+FAEAMCeHFQZk4wBAPbkoGRMmxoAAItRGQMA7InnjAEAsBaPNgEAYDX2jAEAQLRQGQMA7ClgSC4T1W0gdipjkjEAwJ5oUwMAgGihMgaAH9iXv9zqEMKMVoHVIQQFvvtOKl0XpdVMVsaKncqYZAwAsCfa1AAAIFqojAEA9hQwZKrVzGlqAABMMgJHLjPzYwRtagAALEZlDACwJwcd4CIZAwDsiT1jAAAs5qDKmD1jAAAsRmUMALAnQyYr4z6LpN+RjAEA9kSbGgAARAuVMQDAngIBSSZe3BGInZd+kIwBAPZEmxoAAEQLlTEAwJ4cVBmTjAEA9uSgN3DRpgYAwGJUxgAAWzKMgAwTn0E0MzfaSMYAAHsyDHOtZvaMAQAwyTC5ZxxDyZg9YwAALEZlDACwp0BAcpnY92XPGAAAk2hTAwCAaKEyBgDYkhEIyDDRpubRJgAAzKJNDQAAooXKGABgTwFDcjmjMiYZAwDsyTAkmXm0KXaSMW1qAAAsRmUMALAlI2DIMNGmNmKoMiYZAwDsyQjIXJs6dh5tok0NALAlI2CYvnqjsrJSGRkZ8ng88nq92rp164/ev3nzZnm9Xnk8Ho0aNUrLly+PeE2SMQAAR9XU1KioqEgLFixQQ0ODcnNzNWnSJDU2NnZ7/5dffqmrr75aubm5amho0H333ac5c+bo9ddfj2jdqLepj/Xw2w/HTvsAAKwU+O47q0MIOhZLNPZju4wOU63mLn0vSWpvbw8Zd7vdcrvd3c5ZunSppk+frhkzZkiSKioqtGHDBlVVVam8vDzs/uXLl2vkyJGqqKiQJGVmZmrXrl16/PHHdcMNN/Q8WCPKfD7fsVeqcHFxcXHF6OXz+fotT/z3v/81zjjjjD6J89RTTw0bW7hwYbfrdnR0GPHx8cYbb7wRMj5nzhzj4osv7nZObm6uMWfOnJCxN954wxg0aJDR2dnZ4//NUa+MU1NT5fP5lJCQIJfL1euf097errS0NPl8PiUmJvZhhAMLv6ee4ffUM/yeemYg/54Mw9ChQ4eUmprab2t4PB59+eWX6uzsNP2zDMMIyzXHq4pbW1vl9/uVkpISMp6SkqLm5uZu5zQ3N3d7f1dXl1pbWzV8+PAexRn1ZBwXF6cRI0b02c9LTEwccP9n7w/8nnqG31PP8HvqmYH6e0pKSur3NTwejzweT7+v053/Td7dJfQT3d/d+I/hABcAAJKGDRum+Pj4sCq4paUlrPo95owzzuj2/kGDBmno0KE9XptkDACApMGDB8vr9aquri5kvK6uTuPHj+92Tk5OTtj9GzduVHZ2tk466aQerx2zydjtdmvhwoXH7f3jCH5PPcPvqWf4PfUMv6fYVVJSomeffVYrV67Unj17VFxcrMbGRhUUFEiSysrKNG3atOD9BQUF2r9/v0pKSrRnzx6tXLlSK1as0L333hvRui7DiKH3hQEA0M8qKyv12GOPqampSVlZWfrTn/6kiy++WJJ0xx136B//+Ic2bdoUvH/z5s0qLi7Wp59+qtTUVM2fPz+YvHuKZAwAgMVitk0NAMBAQTIGAMBiJGMAACxGMgYAwGIxm4wj/cSV05SXl+uiiy5SQkKCkpOTdf311+vzzz+3OixbKy8vl8vlUlFRkdWh2M6BAwd02223aejQoTr55JN1wQUXqL6+3uqwbKWrq0v333+/MjIyNGTIEI0aNUqLFy9WIMBHcXBiMZmMI/3ElRNt3rxZM2fO1I4dO1RXV6euri7l5eXpm2++sTo0W9q5c6eqq6t13nnnWR2K7Xz99deaMGGCTjrpJL311lv67LPP9MQTT+i0006zOjRbefTRR7V8+XItW7ZMe/bs0WOPPaY//vGPevLJJ60ODTEgJh9t+uUvf6kLL7xQVVVVwbHMzExdf/313X7iCtK///1vJScna/PmzcHn5XDE4cOHdeGFF6qyslIPP/ywLrjgguDn0CCVlpbqvffeo/t0Atdee61SUlK0YsWK4NgNN9ygk08+WS+++KKFkSEWxFxl3NnZqfr6euXl5YWM5+Xlafv27RZFZX9tbW2SpNNPP93iSOxn5syZuuaaa3TFFVdYHYotrV+/XtnZ2brxxhuVnJyscePG6ZlnnrE6LNuZOHGi3nnnHe3du1eS9NFHH2nbtm26+uqrLY4MsSDqX20yqzefuHI6wzBUUlKiiRMnKisry+pwbOWVV17Rhx9+qJ07d1odim198cUXqqqqUklJie677z598MEHmjNnjtxud8hrAZ1u/vz5amtr05gxYxQfHy+/369HHnlEN998s9WhIQbEXDI+JtJPXDnZrFmz9PHHH2vbtm1Wh2IrPp9Pc+fO1caNGy37VFssCAQCys7O1pIlSyRJ48aN06effqqqqiqS8Q/U1NTopZde0urVqzV27Fjt3r1bRUVFSk1N1e233251eLC5mEvGvfnElZPNnj1b69ev15YtW/r0O9IDQX19vVpaWuT1eoNjfr9fW7Zs0bJly9TR0aH4+HgLI7SH4cOH65xzzgkZy8zM1Ouvv25RRPY0b948lZaW6qabbpIknXvuudq/f7/Ky8tJxjihmNsz7s0nrpzIMAzNmjVLb7zxht59911lZGRYHZLtXH755frkk0+0e/fu4JWdna1bb71Vu3fvJhEfNWHChLDH4vbu3av09HSLIrKnb7/9VnFxof9JjY+P59Em9EjMVcbSkU9cTZ06VdnZ2crJyVF1dXXIJ65w5FDS6tWrtW7dOiUkJAQ7CUlJSRoyZIjF0dlDQkJC2B76KaecoqFDh7K3/gPFxcUaP368lixZoilTpuiDDz5QdXW1qqurrQ7NViZPnqxHHnlEI0eO1NixY9XQ0KClS5fqrrvusjo0xAIjRj311FNGenq6MXjwYOPCCy80Nm/ebHVItiKp2+u5556zOjRb+/Wvf23MnTvX6jBs58033zSysrIMt9ttjBkzxqiurrY6JNtpb2835s6da4wcOdLweDzGqFGjjAULFhgdHR1Wh4YYEJPPGQMAMJDE3J4xAAADDckYAACLkYwBALAYyRgAAIuRjAEAsBjJGAAAi5GMAQCwGMkYAACLkYwBALAYyRgAAIuRjAEAsNj/A9dR9K9Xmh0sAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(coverage)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b4097e9-e4a2-468f-a56d-a46e5a12c161",
   "metadata": {},
   "source": [
    "# Compare known result with rio-tiler stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "51cb170d-cf34-4616-9ecb-eed7cc5ed735",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'sum': 3217.5, 'count': 50.0, 'mean': 64.35}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{\n",
    "    \"sum\": weighted_sum,\n",
    "    \"count\": weighted_count,\n",
    "    \"mean\": (weighted_sum / weighted_count)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "aa2cdd56-900e-4d58-9b30-a805d5703d77",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'min': 0.0,\n",
       " 'max': 99.0,\n",
       " 'mean': 58.50450000000001,\n",
       " 'count': 50.004999999999995,\n",
       " 'sum': 3217.7475000000004,\n",
       " 'std': 27.244946990221873,\n",
       " 'median': 62.0,\n",
       " 'majority': 0.0,\n",
       " 'minority': 0.0,\n",
       " 'unique': 55.0,\n",
       " 'histogram': [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10],\n",
       "  [0.0,\n",
       "   9.9,\n",
       "   19.8,\n",
       "   29.700000000000003,\n",
       "   39.6,\n",
       "   49.5,\n",
       "   59.400000000000006,\n",
       "   69.3,\n",
       "   79.2,\n",
       "   89.10000000000001,\n",
       "   99.0]],\n",
       " 'valid_percent': 55.0,\n",
       " 'masked_pixels': 45.0,\n",
       " 'valid_pixels': 55.0,\n",
       " 'percentile_2': 10.08,\n",
       " 'percentile_98': 97.92}"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict(stats[\"b1\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "ed1b3ec4-c0b7-4f4c-a127-febbbcb00379",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "64.34851514848516"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3217.7475000000004 / 50.004999999999995"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd08a8d-641e-4e61-af92-f82bf14861ae",
   "metadata": {},
   "source": [
    "Mean should be sum over count and should also be close to the known result.\n",
    "\n",
    "It is not - seems like there are some bugs in the rio-tiler stats: https://github.com/cogeotiff/rio-tiler/issues/680"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d2d7f953-2dfb-4355-bbe4-a5810eba5570",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "ename": "type",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[26], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m s \u001b[38;5;241m=\u001b[39m stats[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb1\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m s\u001b[38;5;241m.\u001b[39mmean \u001b[38;5;241m==\u001b[39m s\u001b[38;5;241m.\u001b[39msum \u001b[38;5;241m/\u001b[39m s\u001b[38;5;241m.\u001b[39mcount\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "assert stats[\"b1\"].mean == stats[\"b1\"].sum / stats[\"b1\"].count"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}