{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Reading tiled raster data with different zoom levels"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example will show how one can export raster dataset to individual tiles at differnt zoom levels and read the data via the [DataCatalog.get_rasterdataset](../_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.rst) method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "\n",
    "from hydromt import DataCatalog\n",
    "\n",
    "# get some elevation data from the data catalog\n",
    "data_lib = \"artifact_data=v1.0.0\"\n",
    "data_cat = DataCatalog(data_lib)\n",
    "source = \"merit_hydro\"\n",
    "da0 = data_cat.get_rasterdataset(source, variables=[\"elevtn\"])\n",
    "da0.raster.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`da0` is gridded data as an xarray.DataArray object. \n",
    "With HydroMT an xarray.DataArray has some extra functionality via `.raster`\n",
    "This extra functionality does include the ability to write a raster to a tile database (tiling).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiling raster with XYZ stucture\n",
    "\n",
    "First let's have a look at the XYZ structure.\n",
    "an xarray.DataArray is simple written to a tile database in XYZ structure via .raster.to_xyz_tiles\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the database in XYZ structure\n",
    "name = f\"{source}_xyz\"\n",
    "root = join(\"tmpdir\", name)\n",
    "zoom_levels = [0, 1, 2, 3, 4]\n",
    "da0.raster.to_xyz_tiles(\n",
    "    root=root,\n",
    "    tile_size=256,\n",
    "    zoom_levels=zoom_levels,\n",
    "    gdal_driver=\"GTiff\",  # try also 'netcdf4'\n",
    "    compress=\"deflate\",\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tiles in the 'merit_hydro_xyz' folder now contains the zoom_levels as defined above. \n",
    "\n",
    "zoomlevel 0 is at the scale of the xarray.DataArray (one on one). zoomlevel 1 is downscaled by a factor 2 compared to zoomlevel 0. zoomlevel 3 is downscaled by a factor of 8 compared to zoomlevel 0, etc.\n",
    "\n",
    "A mosaic is created per zoomlevel of these tiles in a .vrt file.\n",
    "\n",
    "At last, a .yml file is produced which can be read by the [DataCatalog](../_generated/hydromt.data_catalog.DataCatalog.rst) of HydroMT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiling raster with OSM stucture\n",
    "\n",
    "Now let's have a look at tiling according to the OSM structure\n",
    "an xarray.DataArray is simple written to a tile database in OSM structure via ``.raster.to_slippy_tiles``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the database in XYZ structure\n",
    "name_osm = f\"{source}_osm\"\n",
    "root_osm = join(\"tmpdir\", name_osm)\n",
    "da0.raster.to_slippy_tiles(\n",
    "    root=root_osm,\n",
    "    driver=\"GTiff\",\n",
    "    reproj_method=\"average\",\n",
    "    write_vrt=True,  # try also 'netcdf4'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tiles in the 'merit_hydro_OSM' folder now contains all the zoom_levels from the minimal starting level (9) until the highest level (0).\n",
    "\n",
    "Every tile, regardless of the zoomlevel, has a resolution of 256 by 256 pixels.\n",
    "\n",
    "Zoomlevel 0 is at the scale of the entire world (one tile) and is the most downscaled. Zoomlevel 9 contains the highest resolution (most tiles) in regards to this tile database.\n",
    "\n",
    "A mosaic is created per zoomlevel of these tiles in a .vrt file, if `write_vrt=True` (only applicable to nc and Gtiff files).\n",
    "At last, if `write_vrt=True`,  a yaml file is produced which can be read by the [DataCatalog](../_generated/hydromt.data_catalog.DataCatalog.rst) of HydroMT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiling raster for a webviewer\n",
    "\n",
    "Finally, let's take a look at tiling of a raster dataset with its use being to view the data in a webviewer.\n",
    "This is easily done with the .raster.to_webviewer_tiles method. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the data in OSM structure, but to images!\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "name_png = f\"{source}_png\"\n",
    "root_png = join(\"tmpdir\", name_png)\n",
    "da0.raster.to_slippy_tiles(\n",
    "    root=root_png,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the images are created, let's take a look at an individual tile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the data\n",
    "from PIL import Image\n",
    "\n",
    "# Create a figure to show the image\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches(2.5, 2.5)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_position([0, 0, 1, 1])\n",
    "ax.axis(\"off\")\n",
    "\n",
    "# Show the image\n",
    "im = Image.open(join(root_png, \"9\", \"273\", \"182.png\"))\n",
    "ax.imshow(im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image itself tends to look like an oil painting, but this is correct.\n",
    "\n",
    "The colors for this image are determined so that they are correctly visually represented according to Terrarium Terrain RGB.\n",
    "\n",
    "If one were to see what this would look like in e.g. QGIS, a local sever is needed. \n",
    "With python this is easily done with the command `python -m http.server 8000` from the command line while within the folder where the tiles are located. In this case that would be 'root_png'.\n",
    "In QGIS, make a new XYZ Tiles connection. For this new connection the URL becomes 'http://localhost:8000/{z}/{x}/{y}.png' and the interpolation is set to Terrarium Terrain RGB.\n",
    "\n",
    "However, if the images are meant to be viewed as is, then a custom colormap can be defined to make them look nice!\n",
    "\n",
    "Let's make another dataset of png's!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "name_png_cmap = f\"{source}_png_cmap\"\n",
    "root_png_cmap = join(\"tmpdir\", name_png_cmap)\n",
    "# let us define a nice color for a terrain image\n",
    "da0.raster.to_slippy_tiles(\n",
    "    root=root_png_cmap,\n",
    "    cmap=\"terrain\",\n",
    "    norm=plt.Normalize(vmin=0, vmax=2000),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets take a look at the improved standard visuals!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the data\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "\n",
    "# Create a figure to show the image\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches(2.5, 2.5)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_position([0, 0, 1, 1])\n",
    "ax.axis(\"off\")\n",
    "\n",
    "# Show the image\n",
    "im = Image.open(join(root_png_cmap, \"9\", \"273\", \"182.png\"))\n",
    "ax.imshow(im)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading tiled raster data with zoom levels\n",
    "\n",
    "With [DataCatalog.get_rasterdataset](../_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.rst) a raster (.vrt) can be retrieved. In case of a tile database it can be done for a certain zoomlevel. E.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hydromt import DataCatalog\n",
    "\n",
    "# Load the yml into a DataCatalog\n",
    "data_catalog = DataCatalog(\n",
    "    [join(root, f\"{name}.yml\"), join(root_osm, f\"{name_osm}.yml\")]\n",
    ")\n",
    "\n",
    "# View the structure of the DataCatalog\n",
    "data_catalog.get_source(name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For the osm dataset\n",
    "data_catalog.get_source(name_osm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# without zoom_level the highest res data is fetched\n",
    "da0 = data_catalog.get_rasterdataset(name)\n",
    "da0.raster.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And for OSM\n",
    "da1 = data_catalog.get_rasterdataset(name_osm, zoom=11, geom=da0.raster.box)\n",
    "da1.raster.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Request a raster from the Datacatolog based on zoom resolution & unit\n",
    "da = data_catalog.get_rasterdataset(name_osm, zoom=(1 / 600, \"degree\"))\n",
    "da = data_catalog.get_rasterdataset(name_osm, zoom=(1e4, \"meter\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also directly request a specific zoom level\n",
    "da = data_catalog.get_rasterdataset(name, zoom=zoom_levels[-1], geom=da0.raster.box)\n",
    "print(da.raster.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the data\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n",
    "da0.raster.mask_nodata().plot.imshow(ax=ax, vmin=0, vmax=2500, add_colorbar=False)\n",
    "ax.set_title(f\"zoomlevel {zoom_levels[0]}\")\n",
    "da.raster.mask_nodata().plot.imshow(ax=ax1, vmin=0, vmax=2500, add_colorbar=False)\n",
    "ax1.set_title(f\"zoomlevel {zoom_levels[-1]}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caching tiled raster datasets\n",
    "\n",
    "Tiles of tiled rasterdatasets which are described by a .vrt file can be cached locally (starting from v0.7.0). \n",
    "The requested data tiles will by default be stored to ~/.hydromt_data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set caching to True\n",
    "# NOTE this can also be done at initialization of the DataCatalog\n",
    "data_catalog.cache = True\n",
    "\n",
    "# request some tiles based on bbox -> note the log messages\n",
    "da0 = data_catalog.get_rasterdataset(name, bbox=[11.6, 45.3, 12.0, 46.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if we run the same request again we will use the cached files (and download none)\n",
    "da0 = data_catalog.get_rasterdataset(name, bbox=[11.6, 45.3, 12.0, 46.0])"
   ]
  }
 ],
 "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.11.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "410b670d1c9629c5111115ab6ddd1f95ea051b0f1536068fc1c23c778943ba68"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}