Basic example WRT to TWBD

[mikejohnson51]

key libraries (the last two are lynker/owp developed)

pacman::p_load(sf, terra, dplyr, opendap.catalog, zonal)

read in the area of interest (AOI) from the GDB Justin sent

AOI = read_sf('/Users/mjohnson/Downloads/glfc_n/glfc_n_grid.gdb',
            'glfc_n_grid_poly010620') %>% 
  filter(IBOUND1 == 1)

Get PET data from Gridmet for 2020-01-01

Here we can use any of the 15,000 datasets indexed in our catalog in the same way. Once we align our catalog to the STAC specification, we will also be able to get most imagery datasets like Landsat, sentinel, etc.

That work is being documented here, and will guide some of the New Zealand work.

x = dap(catalog = search('gridMET PET'), 
        AOI = AOI, 
        startDate = "2020-01-01")
#> source:   http://thredds.northwestknowledge.net:8080/thredds/dodsC/agg... 
#> varname(s):
#>    > daily_mean_reference_evapotranspiration_grass [mm] (pet)
#> ==================================================
#> diminsions:  81, 64, 1 (names: lon,lat,day)
#> resolution:  0.042, 0.042, 1 days
#> extent:      -96.91, -93.54, 28.21, 30.88 (xmin, xmax, ymin, ymax)
#> crs:         +proj=longlat +a=6378137 +f=0.00335281066474748 +p...
#> time:        2020-01-01 to 2020-01-01
#> ==================================================
#> values: 5,184 (vars*X*Y*T)

terra::plot(x$daily_mean_reference_evapotranspiration_grass_total)

Putting arbitrary gridded data onto the model grid

The AOI sent along isn't "rectilinear" in that it conforms to the raster/GDAL view of the world. However, we have developed ways to quickly summarize big datasets (both in domain and feature count).

Here we are determining the area weighted mean PET for each model cell. the fun used can be any standard summary function (min, mean, mode, freq, ect) or a user defined function:

o = zonal::execute_zonal(x, AOI, fun = "mean", ID = 'CELL_ID')

plot(o['mean'], border = NA)

Relation to the Hydrofabric

Hydrofabrics are delivered VPU by VPU (VPU = NHDPlusV2 vector processing unit). This Texas AOI falls in VPU-12. We can point the function subset_hy to the VPU-12 geopackage and extract the flowpaths, catchments, poi (points of interest) for the AOI:

nl = subset_hy('/Volumes/Transcend/ngen/CONUS-hydrofabric/refactor/refactor_12.gpkg',  AOI)

sapply(nl, nrow)
#>              fp            cats            pois nexus_locations 
#>            5150            5295             588             396

We can further explore what type of POIs we have indexed (gages, thermoelectric plants, HUC12 outlets, waterbodies, and so on)

table(nl$nexus_locations$type)
#> 
#> Gages HUC12   NID    TE  WBIn WBOut 
#>    84   257    17    16    15     7

Flexible zonal stats and data

Just to emphasize the ability to reuse flexible data, we can generate the mean, catchment wide PET on 2020-01-01 simply by switching out the areal units:

o2 = zonal::execute_zonal(x, nl$cats, fun = "mean", ID = 'ID')
plot(o2['mean'], border = NA)

^{Created on 2022-08-15 by the reprex package (v2.0.1)}