Tip

For an interactive online version click here: Binder badge

Example: Reading tiled raster data with different zoom levels#

This example will show how one can export raster dataset to individual tiles at differnt zoom levels and read the data via the via the DataCatalog.get_rasterdataset method.

[1]:

# Imports
from hydromt import raster, DataCatalog
from hydromt.log import setuplog
from os.path import join

logger = setuplog('tiling', log_level=20)

# get some elevation data from the data catalog
data_lib = "artifact_data"
data_cat = DataCatalog(data_lib, logger=logger)
source = "merit_hydro"
da0 = data_cat.get_rasterdataset(source, variables=["elevtn"])
da0.raster.shape

2023-02-22 04:43:11,125 - tiling - log - INFO - HydroMT version: 0.7.0
2023-02-22 04:43:11,182 - tiling - data_catalog - INFO - Reading data catalog artifact_data v0.0.8 from archive
2023-02-22 04:43:11,182 - tiling - data_catalog - INFO - Parsing data catalog from /home/runner/.hydromt_data/artifact_data/v0.0.8/data_catalog.yml
2023-02-22 04:43:11,257 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro RasterDataset raster data from /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

[1]:

(1920, 1680)

da0 is gridded data as a xarray.DataArray object. With HydroMT a xarray.DataArray has some extra functionality via .raster This extra functionality does include the ability to write a raster to a tile database (tiling).

tiling raster with XYZ stucture#

First let’s have a look at the XYZ structure. A xarray.DataArray is simple written to a tile database in XYZ structure via .raster.to_xyz_tiles

[2]:

# Write the database in XYZ stucture
name = f"{source}_xyz"
root = join("tmpdir", name)
zoom_levels = [0, 1, 2, 3, 5]  # note that you can skip a zoom level
da0.raster.to_xyz_tiles(
    root=root,
    tile_size=256,
    zoom_levels=zoom_levels,
    driver="GTiff",  # try also 'netcdf4'
    compress="deflate",
)

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20..

Tiles at zoomlevel 3 smaller than tile_size 256
Tiles at zoomlevel 5 smaller than tile_size 256

.30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

The tiles in the ‘merit_hydro_xyz’ folder now contains the zoom_levels as defined above.

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.

A mosaic is created per zoomlevel of these tiles in a .vrt file.

At last, a .yml file is produced which can be read by the DataCatalog of HydroMT.

reading tiled raster data with zoom levels#

With DataCatalog.get_rasterdataset a raster (.vrt) can be retreived. In case of a tile database it can be done for a certain zoomlevel. E.g.

[3]:

# Imports
from hydromt import DataCatalog

# Load the yml into a DataCatalog
data_catalog = DataCatalog(join(root, f"{name}.yml"), logger=logger)

# View the structure of the DataCatalog
data_catalog[name]

2023-02-22 04:43:13,135 - tiling - data_catalog - INFO - Parsing data catalog from tmpdir/merit_hydro_xyz/merit_hydro_xyz.yml

[3]:

crs: 4326
data_type: RasterDataset
driver: raster
filesystem: local
path: /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
zoom_levels:
  0: 0.0008333333333325754
  1: 0.0016666666666651508
  2: 0.0033333333333303017
  3: 0.006666666666660603
  5: 0.026666666666642413

[4]:

# without zoom_level the highest res data is fetched
da0 = data_catalog.get_rasterdataset(name)
da0.raster.shape

2023-02-22 04:43:13,145 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt

[4]:

(1920, 1680)

[5]:

# Request a raster from the Datacatolog based on zoom resolution & unit
da = data_catalog.get_rasterdataset(name, zoom_level=(1/600, 'degree'))

da = data_catalog.get_rasterdataset(name, zoom_level=(1e3, 'meter'))

2023-02-22 04:43:13,170 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
2023-02-22 04:43:13,172 - tiling - data_adapter - INFO - Getting data for zoom_level 1 based on res (0.0016666666666666668, 'degree')
2023-02-22 04:43:13,189 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
2023-02-22 04:43:13,191 - tiling - data_adapter - INFO - Getting data for zoom_level 2 based on res (1000.0, 'meter')

[6]:

# We can also directly request a specific zoom level
da = data_catalog.get_rasterdataset(name, zoom_level=zoom_levels[-1])
print(da.raster.shape)

2023-02-22 04:43:13,212 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
(60, 53)

[7]:

# View the data
import matplotlib.pyplot as plt

fig, (ax, ax1) = plt.subplots(1, 2, figsize=(8, 4))
da0.raster.mask_nodata().plot.imshow(ax=ax, vmin=0, vmax=2500, add_colorbar=False)
ax.set_title(f"zoomlevel {zoom_levels[0]}")
da.raster.mask_nodata().plot.imshow(ax=ax1, vmin=0, vmax=2500, add_colorbar=False)
ax1.set_title(f"zoomlevel {zoom_levels[-1]}")

[7]:

Text(0.5, 1.0, 'zoomlevel 5')
../_images/_examples_working_with_tiled_raster_data_12_1.png

Caching tiled raster datasets#

Tiles of tiled rasterdatasets which are described by a .vrt file can be cached locally (starting from v0.7.0). The requested data tiles will by default be stored to ~/.hydromt_data.

[8]:

# set caching to True
# NOTE this can also be done at initialization of the DataCatalog
data_catalog.cache = True

# request some tiles based on bbox -> note the log messages
da0 = data_catalog.get_rasterdataset(name, bbox=[11.6, 45.3, 12.0, 46.0])

2023-02-22 04:43:14,149 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
2023-02-22 04:43:14,156 - tiling - caching - INFO - Downloading 56 tiles to /home/runner/.hydromt_data/merit_hydro_xyz/merit_hydro_xyz

[9]:

# if we run the same request again we will use the cached files (and download none)
da0 = data_catalog.get_rasterdataset(name, bbox=[11.6, 45.3, 12.0, 46.0])

2023-02-22 04:43:14,188 - tiling - data_catalog - INFO - DataCatalog: Getting merit_hydro_xyz RasterDataset raster data from /home/runner/work/hydromt/hydromt/docs/_examples/tmpdir/merit_hydro_xyz/merit_hydro_xyz_zl{zoom_level}.vrt
2023-02-22 04:43:14,193 - tiling - caching - INFO - Downloading 0 tiles to /home/runner/.hydromt_data/merit_hydro_xyz/merit_hydro_xyz