{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "25db29cc",
   "metadata": {},
   "source": [
    "## Accessing High Resolution Electricity Access (HREA) data with the Planetary Computer STAC API\n",
    "\n",
    "The HREA project aims to provide open access to new indicators of electricity access and reliability across the world. Leveraging VIIRS satellite imagery with computational methods, these high-resolution data provide new tools to track progress towards reliable and sustainable energy access across the world.\n",
    "\n",
    "This notebook provides an example of accessing HREA data using the Planetary Computer STAC API."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "202b9ace",
   "metadata": {},
   "source": [
    "### Environment setup\n",
    "\n",
    "This notebook works with or without an API key, but you will be given more permissive access to the data with an API key. The Planetary Computer Hub is pre-configured to use your API key."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7c8c9a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.colors as colors\n",
    "import matplotlib.pyplot as plt\n",
    "import planetary_computer\n",
    "import rasterio\n",
    "import rioxarray\n",
    "\n",
    "import pystac_client\n",
    "from rasterio.plot import show"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94749e92",
   "metadata": {},
   "source": [
    "### Data access\n",
    "\n",
    "The datasets hosted by the Planetary Computer are available from [Azure Blob Storage](https://docs.microsoft.com/en-us/azure/storage/blobs/). We'll use [pystac-client](https://pystac-client.readthedocs.io/) to search the Planetary Computer's [STAC API](https://planetarycomputer.microsoft.com/api/stac/v1/docs) for the subset of the data that we care about, and then we'll load the data directly from Azure Blob Storage. We'll specify a `modifier` so that we can access the data stored in the Planetary Computer's private Blob Storage Containers. See [Reading from the STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) and [Using tokens for data access](https://planetarycomputer.microsoft.com/docs/concepts/sas/) for more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e5198dec",
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "    modifier=planetary_computer.sign_inplace,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eb0ea8f-6c10-48d9-b103-f20aa14ad3af",
   "metadata": {},
   "source": [
    "The HREA dataset covers all of Africa as well as Ecuador. Let's pick up an area of interest that covers Djibouti and query the Planetary Computer API for data coverage for the year 2019."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0d76a979",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Returned 2 Items\n"
     ]
    }
   ],
   "source": [
    "area_of_interest = {\"type\": \"Point\", \"coordinates\": (42.4841, 11.7101)}\n",
    "\n",
    "search = catalog.search(\n",
    "    collections=[\"hrea\"], intersects=area_of_interest, datetime=\"2019-12-31\"\n",
    ")\n",
    "\n",
    "# Check how many items were returned, there could be more pages of results as well\n",
    "items = search.item_collection()\n",
    "print(f\"Returned {len(items)} Items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6796a3a8",
   "metadata": {},
   "source": [
    "We found 3 items for our search. We'll grab jsut the one for Djibouti and see what data assets are available on it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8daca8f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lightscore: Probability of electrification\n",
      "light-composite: Nighttime light annual composite\n",
      "night-proportion: Proportion of nights a settlement is brighter than uninhabited areas\n",
      "estimated-brightness: Estimated brightness levels\n"
     ]
    }
   ],
   "source": [
    "(item,) = [x for x in items if \"Djibouti\" in x.id]\n",
    "data_assets = [\n",
    "    f\"{key}: {asset.title}\"\n",
    "    for key, asset in item.assets.items()\n",
    "    if \"data\" in asset.roles\n",
    "]\n",
    "\n",
    "print(*data_assets, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8dd36b8",
   "metadata": {},
   "source": [
    "### Plotting the data\n",
    "\n",
    "Let's pick the variable `light-composite`, and read in the entire GeoTIFF to plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "56f0f1f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-hrea-hrea-example.ipynb/5.png\"/>"
      ],
      "text/plain": [
       "<Figure size 1400x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "light_comp_asset = item.assets[\"light-composite\"]\n",
    "data_array = rioxarray.open_rasterio(light_comp_asset.href)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(14, 7), dpi=100)\n",
    "show(\n",
    "    data_array.data,\n",
    "    ax=ax,\n",
    "    norm=colors.PowerNorm(1, vmin=0.01, vmax=1.4),\n",
    "    cmap=\"magma\",\n",
    "    title=\"Djibouti (2019)\",\n",
    ")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db330c2f",
   "metadata": {},
   "source": [
    "### Read a window\n",
    "\n",
    "Cloud Optimized GeoTIFFs (COGs) allows us to efficiently download and read sections of a file, rather than the entire file, when only part of the region is required. The COGs are stored on disk with an internal set of windows. You can read sections of any shape and size, but reading them in the file-defined window size is most efficient. Let's read the same asset, but this time only request the second window. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "59989159",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available windows:\n",
      "((0, 0), Window(col_off=0, row_off=0, width=256, height=256))\n",
      "((0, 1), Window(col_off=256, row_off=0, width=142, height=256))\n",
      "((1, 0), Window(col_off=0, row_off=256, width=256, height=165))\n",
      "((1, 1), Window(col_off=256, row_off=256, width=142, height=165))\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-hrea-hrea-example.ipynb/6.png\"/>"
      ],
      "text/plain": [
       "<Figure size 1400x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Reading only the second window of the file, as an example\n",
    "i_window = 2\n",
    "with rasterio.open(light_comp_asset.href) as src:\n",
    "    windows = list(src.block_windows())\n",
    "    print(\"Available windows:\", *windows, sep=\"\\n\")\n",
    "    _, window = windows[i_window]\n",
    "    section = data_array.rio.isel_window(window)\n",
    "\n",
    "fig, xsection = plt.subplots(1, 1, figsize=(14, 7))\n",
    "show(\n",
    "    section.data,\n",
    "    ax=xsection,\n",
    "    norm=colors.PowerNorm(1, vmin=0.01, vmax=1.4),\n",
    "    cmap=\"magma\",\n",
    "    title=\"Reading a single window\",\n",
    ")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f51acaa4",
   "metadata": {},
   "source": [
    "### Zoom in on a region within the retrieved window\n",
    "\n",
    "Let's plot the region around the city of Dikhil, situated within that second data window, around this bounding box (in x/y coordinates, which is latitude / longitude):\n",
    "\n",
    "```\n",
    "(42.345868941491204, 11.079694223371735, 42.40420227530527, 11.138027557181712)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "75090aea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-hrea-hrea-example.ipynb/7.png\"/>"
      ],
      "text/plain": [
       "<Figure size 1400x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, xsection = plt.subplots(1, 1, figsize=(14, 7))\n",
    "show(\n",
    "    section.sel(\n",
    "        x=slice(42.345868941491204, 42.40420227530527),\n",
    "        y=slice(11.138027557181712, 11.079694223371735),\n",
    "    ).data,\n",
    "    ax=xsection,\n",
    "    norm=colors.PowerNorm(1, vmin=0.01, vmax=1.4),\n",
    "    cmap=\"magma\",\n",
    "    title=\"Dikhil (2019)\",\n",
    ")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d809171",
   "metadata": {},
   "source": [
    "### Plot change over time\n",
    "\n",
    "The HREA dataset goes back several years. Let's search again for the same area, but this time over a longer temporal span."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "41319df7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Returned 8 Items:\n"
     ]
    }
   ],
   "source": [
    "search = catalog.search(\n",
    "    collections=[\"hrea\"], intersects=area_of_interest, datetime=\"2012-12-31/2019-12-31\"\n",
    ")\n",
    "items = [item.to_dict() for item in search.get_items() if \"Djibouti\" in item.id]\n",
    "print(f\"Returned {len(items)} Items:\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9421ec36",
   "metadata": {},
   "source": [
    "We got 8 items this time, each corresponding to single year. To plot the change of light intensity over time, we'll open the same asset on each of these year-items and read in the window with Dikhil. Since we're using multiple items, we'll use `stackstac` to stack them together into a single DataArray for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d033e014",
   "metadata": {},
   "outputs": [],
   "source": [
    "import stackstac\n",
    "\n",
    "bounds_latlon = (\n",
    "    42.345868941491204,\n",
    "    11.079694223371735,\n",
    "    42.40420227530527,\n",
    "    11.138027557181712,\n",
    ")\n",
    "dikhil = (\n",
    "    stackstac.stack(items, assets=[\"light-composite\"], bounds_latlon=bounds_latlon)\n",
    "    .squeeze()\n",
    "    .compute()\n",
    "    .quantile(0.9, dim=[\"y\", \"x\"])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2e89d342",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-hrea-hrea-example.ipynb/10.png\"/>"
      ],
      "text/plain": [
       "<Figure size 1200x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "dikhil.plot(ax=ax)\n",
    "\n",
    "ax.set(title=\"Dikhil composite light output\", ylabel=\"Annual light output, normalize\");"
   ]
  }
 ],
 "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.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {
     "d653d72f680241f5aca06f6787af40ba": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "layout": "IPY_MODEL_f740e45becec4081a84ea3efec0fb70a"
      }
     },
     "f740e45becec4081a84ea3efec0fb70a": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {}
     }
    },
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}