zonal
is an active package for intersecting vector aggregation units
with large gridded data. While there are many libraries that seek to
tackle this problem (see credits) we needed a library that could handle
large gridded extents storing categorical and continuous data, with
multiple time layers with both many small vector units and few large
units.
We also seek to segment the creation of grid weights from the zonal execution so that the same weight map can be applied across different products with the same structure.
You can install the development version of zonal
from
GitHub with:
# install.packages("remotes")
remotes::install_github("mikejohnson51/zonal")
This is a basic example that takes a NetCDF file containing a 4km grid for the continental USA and daily precipitation for the year 1979 (365 layers). Our goal is to subset this file to the southern USA, and compute daily county level averages. The result is a daily rainfall average for each county.
library(zonal)
library(dplyr)
library(tidyr)
library(ggplot2)
file <- 'to_build/pr_2020.nc'
AOI <- AOI::aoi_get(state = "south", county = "all")
system.time({
pr_zone = execute_zonal(file, geom = AOI, ID = "geoid", join = FALSE)
})
#> user system elapsed
#> 9.135 2.036 11.298
# PET zone: Counties, time slices/ID
dim(pr_zone)
#> [1] 1421 367
x = merge(AOI, pr_zone, by ="geoid")
# Plot Day with the maximum single county max rainfall.
n = colnames(pr_zone)[which(pr_zone[,-1] == max(pr_zone[,-1]), arr.ind = TRUE)[2] + 1]
ggplot(data = x) +
geom_sf(aes(fill = get(n)), color = NA) +
scale_fill_viridis_c() +
theme_void() +
labs(fill = "PR (mm)")
# Plot Day with the maximum county wide rainfall
n2 = names(which.max(colSums(dplyr::select(pr_zone, -geoid))))
terra::plot(terra::rast(file)[[103]])
ggplot() +
geom_sf(data = x, aes(fill = get(n2)), color = NA) +
scale_fill_viridis_c() +
theme_void() +
labs(fill = "PR (mm)")
data = pr_zone %>%
slice_max(rowSums(select(., -geoid))) %>%
pivot_longer(-geoid, names_to = "day", values_to = "prcp") %>%
mutate(day = as.numeric(gsub("V","", day)))
head(data)
#> # A tibble: 6 × 3
#> geoid day prcp
#> <chr> <dbl> <dbl>
#> 1 37175 1 0
#> 2 37175 2 26.7
#> 3 37175 3 13.6
#> 4 37175 4 0.234
#> 5 37175 5 0.144
#> 6 37175 6 2.52
One of the largest limitations of existing utilities is the ability to handle categorical data. Here we show an example for a 1km grid storing land cover data from MODIS. This grid was creating by mos acing 19 MODIS tiles covering CONUS. The summary function for this categorical frequency is “freq”.
file = 'to_build/2019-01-01.tif'
rcl = read.csv("to_build/modis_lc.csv") %>%
dplyr::select(from = Class, to = short)
system.time({
lc = execute_zonal(file, AOI, "geoid", FUN = 'freq', rcl = rcl)
})
#> user system elapsed
#> 6.453 0.459 6.998
Here lets look at a quick intergation of the AOI
/climateR
/zonal
family. The goal is to find monthly mean, normal (1981-2010), rainfall
for all USA counties in the south.
We can do this by (1) defining the southern AOI
, (2) queuing normal
rainfall data from climateR
and (3) passing the returned data through
zonal
:
system.time({
tt = climateR::getTerraClimNormals(AOI, "prcp")
out = zonal::execute_zonal(file = tt[[1]],
geom = AOI,
ID = "geoid",
join = TRUE)
})
#> Spherical geometry (s2) switched off
#> Spherical geometry (s2) switched on
#> user system elapsed
#> 4.755 0.449 6.864
plot(out[paste0("V", 1:12)], border = NA)
- Code style should attempt to follow the tidyverse style guide.
- Please avoid adding significant new dependencies without a documented reason why.
- Please attempt to describe what you want to do prior to contributing by submitting an issue.
- Please follow the typical github fork - pull-request workflow.
- Make sure you use roxygen and run Check before contributing.
Similar R packages:
Logo Artwork: Justin Singh-Mohudpur