NWM Channel Routing Application • zonal

The purpose here is to demonstrate a real world application of zonal to summarize the NWM channel routing parameters to a in-development hydrofabric.

library(terra)
library(sf)
library(zonal)
library(paint)

Define Aggregations Geometries and Geogrid

Here we are using a version of the ngen hydrofabric:

geom <- read_sf('hydrofabric.gpkg', "catchments")
paint(geom)

## sf [18041, 4] 
## active geometry column: geom (POLYGON)
## crs: 5070 (NAD83 / Conus Albers)
## crs unit: metre 
## ID        chr cat-1 cat-2 cat-4 cat-5 cat-6 cat-7
## area_sqkm dbl 12.457576 267.083595 8.319214 9.278138 60.577~
## toID      chr nex-2 nex-3 nex-5 nex-6 nex-7 nex-8
## geom      sfc POLY 2,024B POLY 9,064B POLY 1,656B POLY 1,81~

And the soilproperties_CONUS_FullRouting file from NWM v2.16. Can download from here

(f = rast('soilproperties_CONUS_FullRouting.nc'))

## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named west_east BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
## [1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named south_north BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"

## class       : SpatRaster 
## dimensions  : 3840, 4608, 45  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 0.5, 4608, 0.5, 3840  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## sources     : soilproperties_CONUS_FullRouting.nc:bexp  (4 layers) 
##               soilproperties_CONUS_FullRouting.nc:cwpvt  
##               soilproperties_CONUS_FullRouting.nc:dksat  (4 layers) 
##               ... and 15 more source(s)
## varnames    : bexp 
##               cwpvt 
##               dksat 
##               ...
## names       : bexp_~ime=1, bexp_~ime=1, bexp_~ime=1, bexp_~ime=1, cwpvt_Time=1, dksat~ime=1, ...

Note that when the layered attributes are separated, there are 45 raster layers. Also note that the spatial metadata for the NWM file is out of wack (see the resolution, extent, and CRS). We will fix this using some functionality from wrfhydroSubsetter to identify the layer structure.