{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Example: Reading raster data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "This example illustrates the how to read raster data using the HydroMT [DataCatalog](../_generated/hydromt.data_catalog.DataCatalog.rst) with the `raster`, `netcdf` and `raster_tindex` drivers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hydromt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download artifacts for the Piave basin to `~/.hydromt_data/`.\n",
    "data_catalog = hydromt.DataCatalog(data_libs=[\"artifact_data=v1.0.0\"])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Rasterio driver"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "To read raster data and parse it into a [xarray Dataset or DataArray](https://xarray.pydata.org/en/stable/user-guide/data-structures.html) we use the [get_rasterdataset()](../_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.rst) method. All `driver_kwargs` in the data catalog yaml file will be passed to this method. The `raster` driver supports all [GDAL data formats ](http://www.gdal.org/formats_list.html), including the often used GeoTiff of Cloud Optimized GeoTiff (COG) formats. Tiled datasets can also be passed as a [virtual raster tileset (vrt) file](https://gdal.org/drivers/raster/vrt.html). \n",
    "\n",
    "As an example we will use the [MERIT Hydro](http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro) dataset which is a set of GeoTiff files with identical grids, one for each variable of the datasets. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# inspect data source entry in data catalog yaml file\n",
    "data_catalog.get_source(\"merit_hydro\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "We can load any RasterDataset using [DataCatalog.get_rasterdataset()](../_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.rst). Note that if we don't provide any arguments it returns the full dataset with nine data variables and for the full spatial domain. Only the data coordinates are actually read, the data variables are still lazy [Dask arrays](https://docs.dask.org/en/stable/array.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = data_catalog.get_rasterdataset(\"merit_hydro\")\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "The data can be visualized with the [.plot()](https://docs.xarray.dev/en/latest/generated/xarray.DataArray.plot.html) xarray method. We replace all nodata values with NaNs with [.raster.mask_nodata()](../_generated/hydromt.gis.DataArray.raster.mask_nodata.rst)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds[\"elevtn\"].raster.mask_nodata().plot(cmap=\"terrain\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "We can request a (spatial) subset data by providing additional `variables` and `bbox` / `geom` arguments. Note that these return a smaller spatial extent and just two data variables. The variables argument is especially useful if each variable of the dataset is saved in a separate file and the `{variable}` key is used in the path argument of the data source (see above) to limit which files are actually read. If a single variable is requested a DataArray instead of a Dataset is returned unless the `single_var_as_array` argument is set to False (True by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = [11.70, 45.35, 12.95, 46.70]\n",
    "ds = data_catalog.get_rasterdataset(\n",
    "    \"merit_hydro\", bbox=bbox, variables=[\"elevtn\"], single_var_as_array=True\n",
    ")\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "TIP: To write a dataset back to a stack of raster in a single folder use the [.raster.to_mapstack()](../_generated/hydromt.gis.Dataset.raster.to_mapstack.rst) method."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Xarray driver"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Many gridded datasets with a third dimension (e.g. time) are saved in netcdf or zarr files, which can be read with Xarray. This data is read using the [xarray.open_mfdataset()](https://docs.xarray.dev/en/latest/generated/xarray.open_mfdataset.html) method. These formats are flexible and therefore \n",
    "HydroMT is not always able to read the geospatial attributes such as the CRS from the data and it has to be set through the data catalog [yaml file](../guides/advanced_user/data_prepare_cat.rst).\n",
    "\n",
    "If the data is stored per year or month, the `{year}` and `{month}` keys can be used in the path argument of a data source in the data catalog yaml file to speed up the reading of a temporal subset of the data using the `date_tuple` argument of [DataCatalog.get_rasterdataset()](../_generated/hydromt.data_catalog.DataCatalog.get_rasterdataset.rst) (not in this example).\n",
    "\n",
    "As example we use the [ERA5](https://doi.org/10.24381/cds.bd0915c6) dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note the crs argument as this is missing in the original data\n",
    "data_catalog.get_source(\"era5_hourly\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the some units are converted\n",
    "ds = data_catalog.get_rasterdataset(\"era5_hourly\")\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Raster_tindex Resolver"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "If the raster data is tiled but for each tile a different CRS is used (for instance a different UTM projection for each UTM zone), this dataset cannot be described using a VRT file. In this case a vector file can be built to use a raster tile index using [gdaltindex](https://gdal.org/programs/gdaltindex.html). To read the data into a single `xarray.Dataset` the data needs to be reprojected and mosaicked to a single CRS while reading. As this type of data cannot be loaded lazily the method is typically used with an area of interest for which the data is loaded and combined. \n",
    "\n",
    "As example we use the [GRWL mask](https://doi.org/10.5281/zenodo.1297434) raster tiles for which we have created a tileindex using the aforementioned *gdaltindex* command line tool. Note that the path points to the GeoPackage output of the *gdaltindex* tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_catalog.get_source(\"grwl_mask\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the tileindex is a GeoPackage vector file\n",
    "# with an attribute column 'location' (see also the tileindex argument under driver_kwargs) containing the (relative) paths to the raster file data\n",
    "import geopandas as gpd\n",
    "\n",
    "tindex_path = data_catalog.get_source(\"grwl_mask\").full_uri\n",
    "print(tindex_path)\n",
    "gpd.read_file(tindex_path, rows=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this returns a DataArray (single variable) wit a mosaic of several files (see source_file attribute)\n",
    "ds = data_catalog.get_rasterdataset(\"grwl_mask\", bbox=bbox)\n",
    "ds"
   ]
  }
 ],
 "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": "3808d5b5b54949c7a0a707a38b0a689040fa9c90ab139a050e41373880719ab1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}