Example: Working with raster data#

# get the GHS population raster dataset from the catalog
da = data_catalog.get_rasterdataset("ghs_pop_2015").rename("population")
da.raster.mask_nodata().reset_coords(drop=True).plot(vmax=200)

# coordinate reference system
da.raster.crs

# geospatial transform, see https://www.perrygeo.com/python-affine-transforms.html
da.raster.transform

# names of x- and y dimensions
(da.raster.x_dim, da.raster.y_dim)

# nodata value (or fillvalue) of a specific variable
da.raster.nodata

from hydromt.gis import utm_crs

utm = utm_crs(da.raster.bounds)
print(f"Destination CRS: {utm}")
da_pop = da.raster.mask_nodata()
da_pop_dens = da_pop.raster.density_grid().rename("population_dens")  # pop.m-2
da_pop_dens_utm = da_pop_dens.raster.reproject(
    method="bilinear", dst_crs=utm, dst_res=250
)
da_pop_utm = da_pop_dens_utm * da_pop_dens_utm.raster.area_grid()  # pop
bias = (da_pop.sum() / da_pop_utm.sum()).compute().item()
print(f"Error: {(1-bias)*100:.3f}%")
da_pop_utm_adj = da_pop_utm * bias  # bias correct
da_pop_utm_adj.name = "population"
da_pop_utm_adj.reset_coords(drop=True).plot(vmax=200)

import numpy as np
import xarray as xr

gdf = data_catalog.get_geodataframe("gadm_level3", variables=["NAME_3"])
ds = xr.merge(
    [
        da_pop.raster.zonal_stats(gdf, stats=["sum"]) / 1e3,  # count [pop x1000]
        da_pop_dens.raster.zonal_stats(gdf, stats=["max", np.nanmax])
        * 1e6,  # density [pop km-2]
    ]
)
for var in ds.data_vars:
    gdf[var] = ds[var]

gdf.sort_values("population_sum", ascending=False).head()

gdf.plot(
    "population_sum",
    scheme="NaturalBreaks",
    legend=True,
    legend_kwds=dict(fmt="{:.0f}", title="population [x1000]"),
    figsize=(6, 6),
)

da_soil = data_catalog.get_rasterdataset(
    "soilgrids", bbox=[12.7, 45.6, 13, 45.8], variables=["ph_sl1"]
)
da_soil.raster.mask_nodata().plot()

# Note that data is only interpolated leaving a small nodata gap in the lower right corner. This can be extrapolated with the 'nearest' method.
da_soil.raster.interpolate_na(method="linear").raster.mask_nodata().plot()

# step 1: read the data and mask nodata values
bbox = [12.2, 45.3, 13, 45.8]
da_dem = data_catalog.get_rasterdataset("merit_hydro", variables=["elevtn"], bbox=bbox)
da_dem = da_dem.raster.mask_nodata()
da_dem.attrs.update(units="m+EGM96")
print(f"resolution MERIT Hydro: {da_dem.raster.res[0]:.05f}")
da_bath = data_catalog.get_rasterdataset("gebco").raster.mask_nodata()
print(f"resolution GEBCO: {da_bath.raster.res[0]:.05f}")
da_mdt = data_catalog.get_rasterdataset("mdt_cnes_cls18").raster.mask_nodata()
print(f"resolution MDT: {da_mdt.raster.res[0]:.05f}")

plot_kwargs = dict(
    vmin=-50,
    vmax=50,
    cmap="coolwarm",
)
da_dem.plot(**plot_kwargs)

# step 2: convert GEBCO to EGM96 ref and reproject to MERIT Hydro grid
da_bath_egm = da_bath + da_mdt.raster.reproject_like(da_bath)
da_bath_reproj = da_bath_egm.raster.reproject_like(da_dem, method="cubic")
print(f"resolution reprojected GEBCO: {da_bath_reproj.raster.res[0]:.05f}")
da_dem_merged = da_dem.where(da_dem.notnull(), da_bath_reproj)
da_dem_merged.plot(**plot_kwargs)