generated from NOAA-OWP/owp-open-source-project-template
-
Notifications
You must be signed in to change notification settings - Fork 4
/
README.Rmd
165 lines (120 loc) · 5.54 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%",
warning = FALSE,
message = FALSE
)
library(dplyr)
library(tidyr)
library(ggplot2)
library(terra)
```
# zonal <img src="man/figures/logo.png" align="right" alt="" width="120" />
<!-- badges: start -->
[](https://github.com/mikejohnson51/zonal/actions/workflows/R-CMD-check.yaml)
[](https://www.repostatus.org/#active)
[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
[](#)
[](https://choosealicense.com/licenses/mit/)
[](https://github.com/mikejohnson51/zonal/actions/workflows/pkgdown.yaml)
<!-- badges: end -->
`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.
The package offers 3 main options through a common syntax:
1. The ability to pregenerate weighting grids and applying those over large datasets distrbuted across files (e.g. 30 years of daily data stored in annula files). Rapid data summarization is supported by collapse and data.table.
2. Thin wrappers over pure exact_extract() when appropriate (e.g. calling a core `exactrextractr` function)
3. Flexible custom functions that are easily applied over multi layer files (e.g. geometric means and circular means)
## Installation
You can install the development version of `zonal` with:
``` r
# install.packages("remotes")
remotes::install_github("mikejohnson51/zonal")
```
## Example
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.
```{r example}
library(zonal)
AOI <- AOI::aoi_get(state = "south", county = "all")
d = rast("to_build/pr_2022.nc")
system.time({
pr_zone <- execute_zonal(data = d,
geom = AOI,
ID = "fip_code",
join = TRUE)
})
```
### Daily maximum mean rainfall in the South?
```{r}
n = names(which.max(colSums(as.data.frame(pr_zone)[,grepl('precipitation', names(pr_zone))])))
ggplot(data = pr_zone) +
geom_sf(aes(fill = get(n)), color = NA) +
scale_fill_viridis_c() +
theme_void() +
labs(fill = "PR (mm)")
```
### Timeseries of county with maximum annual rainfall
```{r}
data <- pr_zone %>%
as.data.frame() %>%
slice_max(rowSums(select(., c(starts_with('mean'))))) %>%
select(fip_code, starts_with('mean')) %>%
pivot_longer(-fip_code, names_to = "day", values_to = "prcp") %>%
mutate(day = as.numeric(gsub("mean.precipitation_amount_day.", "", day)))
head(data)
```
```{r, echo = FALSE}
ggplot(data, aes(x = day, y = prcp)) +
geom_line() +
labs(
x = "Date", y = "Mean Rainfall",
title = paste("GEOID: ", data$fip_code[1])
) +
theme_bw()
```
# 1km Landcover Grid (Categorical)
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 mosacing 19 MODIS tiles covering CONUS. The summary function for this categorical frequency is "freq".
```{r}
system.time({
lc <- execute_zonal(data = rast("to_build/2019-01-01.tif"),
geom = AOI,
ID = "fip_code",
fun = "frac")
})
```
## Zonal and opendap.catalog
Here lets look at a quick integration of the `AOI`/`opendap.catalog`/`zonal` family. The goal is to find monthly mean, normal (1981-2010), rainfall for all USA counties in the south.
```{r, message = FALSE, warning = FALSE}
library(climateR)
AOI <- AOI::aoi_get(state = "FL", county = "all")
system.time({
data <- climateR::dap(
URL = "https://cida.usgs.gov/thredds/dodsC/bcsd_obs",
AOI = AOI,
startDate = "1995-01-01",
verbose = FALSE,
varname = "pr"
) |>
execute_zonal(geom = AOI, ID = "fip_code", join = TRUE)
})
plot(data[grepl("mean", names(data))], border = NA)
```
----
## Getting involved
1. Code style should attempt to follow the tidyverse style guide.
2. Please avoid adding significant new dependencies without a documented reason why.
3. Please attempt to describe what you want to do prior to contributing by submitting an issue.
4. Please follow the typical github fork - pull-request workflow.
5. Make sure you use roxygen and run Check before contributing.
----
## Credits and references
Similar R packages / Core Dependencies:
1. [exactexactr](https://github.com/isciences/exactextractr)
4. [sf](https://github.com/r-spatial/sf)
5. [terra](https://github.com/rspatial/raster)
**Logo Artwork:** [Justin Singh-Mohudpur](https://www.justinsingh.me/about/)