diff --git a/DESCRIPTION b/DESCRIPTION index b7c101c8..dfc5a769 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,7 @@ Authors@R: email = "kyle.messier@nih.gov")) Description: Near Real Time air pollution model results and code produced by the SET group. It is fully tested, versioned, and open source and open access. Depends: R (>= 4.1.0) -Imports: dplyr, sf, stats, terra, methods, BART, data.table, spData +Imports: dplyr, sf, stats, terra, methods, BART, data.table, spData, exactextractr, ggplot2, scatterpie Suggests: covr, withr, knitr, rmarkdown, testthat (>= 3.0.0), stars, sftime Encoding: UTF-8 VignetteBuilder: knitr, rmarkdown diff --git a/R/calc_nlcd_ratio.R b/R/calc_nlcd_ratio.R new file mode 100644 index 00000000..08fce25a --- /dev/null +++ b/R/calc_nlcd_ratio.R @@ -0,0 +1,86 @@ +# needs terra, exactextractr, spData, sf + +globalVariables(c("us_states", "%>%")) + + +#' Compute land cover classes ratio in circle buffers around points +#' +#' @param data_vect terra::SpatVector of points geometry +#' @param buf_radius numeric (non-negative) giving the +#' radius of buffer around points +#' @param year numeric giving the year of NLCD data used +#' @param nlcd_path character giving nlcd data path +#' @export +calc_nlcd_ratio <- function(data_vect, + buf_radius = 1000, + year = 2021, + nlcd_path) { + # check inputs + if (!is.numeric(buf_radius)) { + stop("buf_radius is not a numeric.") + } + if (buf_radius <= 0) { + stop("buf_radius has not a likely value.") + } + if (!is.numeric(year)) { + stop("year is not a numeric.") + } + if (class(data_vect)[1] != "SpatVector") { + stop("data_vect is not a terra::SpatVector.") + } + if (!is.character(nlcd_path)) { + stop("nlcd_path is not a character.") + } + if (!file.exists(nlcd_path)) { + stop("nlcd_path does not exist.") + } + # open nlcd file corresponding to the year + nlcd_file <- list.files(nlcd_path, + pattern = paste0("nlcd_", year, "_.*.tif$"), + full.names = TRUE) + if (length(nlcd_file) == 0) { + stop("NLCD data not available for this year.") + } + nlcd <- rast(nlcd_file) + # select points within mainland US and reproject on nlcd crs if necessary + # need spData library + data(us_states, package = "spData") + us_main <- sf::st_union(us_states) %>% + terra::vect() %>% + terra::project(y = crs(data_vect)) + data_vect_b <- data_vect %>% + terra::intersect(x = us_main) + if (!terra::same.crs(data_vect_b, nlcd)) { + data_vect_b <- terra::project(data_vect_b, crs(nlcd)) + } + # create circle buffers with buf_radius + bufs_pol <- terra::buffer(data_vect_b, width = buf_radius) %>% + sf::st_as_sf() + # ratio of each nlcd class per buffer + nlcd_at_bufs <- exactextractr::exact_extract(nlcd, + st_geometry(bufs_pol), + fun = "frac", + stack_apply = TRUE, + progress = FALSE) + # select only the columns of interest + nlcd_at_bufs <- nlcd_at_bufs[names(nlcd_at_bufs)[grepl("frac_", + names(nlcd_at_bufs))]] + # change column names + fpath <- system.file("extdata", "nlcd_classes.csv", package = "NRTAPmodel") + nlcd_classes <- read.csv(fpath) + nlcd_names <- names(nlcd_at_bufs) %>% + sub(pattern = "frac_", replacement = "") %>% + as.numeric() + nlcd_names <- nlcd_classes[nlcd_classes$value %in% nlcd_names, c("class")] + new_names <- sapply( + nlcd_names, + function(x) { + paste0("frac_", x, "_", year, "_", buf_radius, "m") + } + ) + names(nlcd_at_bufs) <- new_names + # merge data_vect with nlcd class fractions (and reproject) + new_data_vect <- cbind(data_vect_b, nlcd_at_bufs) %>% + terra::project(crs(data_vect)) + return(new_data_vect) +} diff --git a/inst/extdata/nlcd_classes.csv b/inst/extdata/nlcd_classes.csv new file mode 100644 index 00000000..90dd7b9f --- /dev/null +++ b/inst/extdata/nlcd_classes.csv @@ -0,0 +1,17 @@ +"","value","class","names","col" +"1",0,"Unc","Unclassified","white" +"2",11,"WTR","Open Water","#476ba1" +"3",21,"OSD","Developed, Open Space","#decaca" +"4",22,"LID","Developed, Low Intensity","#d99482" +"5",23,"MID","Developed, Medium Intensity","#ee0000" +"6",24,"HID","Developed, High Intensity","#ab0000" +"7",31,"BRN","Barren Land","#b3aea3" +"8",41,"DFO","Deciduous Forest","#68ab63" +"9",42,"EFO","Evergreen Forest","#1c6330" +"10",43,"MFO","Mixed Forest","#b5ca8f" +"11",52,"SHB","Shrub/Scrub","#ccba7d" +"12",71,"GRS","Herbaceous","#e3e3c2" +"13",81,"PAS","Hay/Pasture","#dcd93d" +"14",82,"CRP","Cultivated Crops","#ab7028" +"15",90,"WDW","Woody Wetlands","#bad9eb" +"16",95,"EHW","Emergent Herbaceous Wetlands","#70a3ba" diff --git a/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif b/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif new file mode 100644 index 00000000..ec459f54 Binary files /dev/null and b/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif differ diff --git a/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif.aux.xml b/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif.aux.xml new file mode 100644 index 00000000..85b33d64 --- /dev/null +++ b/tests/testdata/nlcd_2019_land_cover_l48_20210604.tif.aux.xml @@ -0,0 +1,2346 @@ + + + NLCD Land Cover Class + + + value + 0 + 0 + + + NLCD Land Cover Class + 2 + 0 + + + Histogram + 1 + 0 + + + Red + 0 + 0 + + + Green + 0 + 0 + + + Blue + 0 + 0 + + + Opacity + 0 + 0 + + + 0 + Unclassified + 7853863229 + 0 + 0 + 0 + 0 + + + 1 + + 0 + 0 + 0 + 0 + 255 + + + 2 + + 0 + 0 + 0 + 0 + 255 + + + 3 + + 0 + 0 + 0 + 0 + 255 + + + 4 + + 0 + 0 + 0 + 0 + 255 + + + 5 + + 0 + 0 + 0 + 0 + 255 + + + 6 + + 0 + 0 + 0 + 0 + 255 + + + 7 + + 0 + 0 + 0 + 0 + 255 + + + 8 + + 0 + 0 + 0 + 0 + 255 + + + 9 + + 0 + 0 + 0 + 0 + 255 + + + 10 + + 0 + 0 + 0 + 0 + 255 + + + 11 + Open Water + 472399232 + 70 + 107 + 159 + 255 + + + 12 + Perennial Snow/Ice + 962418 + 209 + 222 + 248 + 255 + + + 13 + + 0 + 0 + 0 + 0 + 255 + + + 14 + + 0 + 0 + 0 + 0 + 255 + + + 15 + + 0 + 0 + 0 + 0 + 255 + + + 16 + + 0 + 0 + 0 + 0 + 255 + + + 17 + + 0 + 0 + 0 + 0 + 255 + + + 18 + + 0 + 0 + 0 + 0 + 255 + + + 19 + + 0 + 0 + 0 + 0 + 255 + + + 20 + + 0 + 0 + 0 + 0 + 255 + + + 21 + Developed, Open Space + 240566180 + 222 + 197 + 197 + 255 + + + 22 + Developed, Low Intensity + 153288747 + 217 + 146 + 130 + 255 + + + 23 + Developed, Medium Intensity + 92578072 + 235 + 0 + 0 + 255 + + + 24 + Developed, High Intensity + 33121466 + 171 + 0 + 0 + 255 + + + 25 + + 0 + 0 + 0 + 0 + 255 + + + 26 + + 0 + 0 + 0 + 0 + 255 + + + 27 + + 0 + 0 + 0 + 0 + 255 + + + 28 + + 0 + 0 + 0 + 0 + 255 + + + 29 + + 0 + 0 + 0 + 0 + 255 + + + 30 + + 0 + 0 + 0 + 0 + 255 + + + 31 + Barren Land + 87406005 + 179 + 172 + 159 + 255 + + + 32 + + 0 + 0 + 0 + 0 + 255 + + + 33 + + 0 + 0 + 0 + 0 + 255 + + + 34 + + 0 + 0 + 0 + 0 + 255 + + + 35 + + 0 + 0 + 0 + 0 + 255 + + + 36 + + 0 + 0 + 0 + 0 + 255 + + + 37 + + 0 + 0 + 0 + 0 + 255 + + + 38 + + 0 + 0 + 0 + 0 + 255 + + + 39 + + 0 + 0 + 0 + 0 + 255 + + + 40 + + 0 + 0 + 0 + 0 + 255 + + + 41 + Deciduous Forest + 833976610 + 104 + 171 + 95 + 255 + + + 42 + Evergreen Forest + 1033039764 + 28 + 95 + 44 + 255 + + + 43 + Mixed Forest + 305029988 + 181 + 197 + 143 + 255 + + + 44 + + 0 + 0 + 0 + 0 + 255 + + + 45 + + 0 + 0 + 0 + 0 + 255 + + + 46 + + 0 + 0 + 0 + 0 + 255 + + + 47 + + 0 + 0 + 0 + 0 + 255 + + + 48 + + 0 + 0 + 0 + 0 + 255 + + + 49 + + 0 + 0 + 0 + 0 + 255 + + + 50 + + 0 + 0 + 0 + 0 + 255 + + + 51 + + 0 + 0 + 0 + 0 + 255 + + + 52 + Shrub/Scrub + 1961779404 + 204 + 184 + 121 + 255 + + + 53 + + 0 + 0 + 0 + 0 + 255 + + + 54 + + 0 + 0 + 0 + 0 + 255 + + + 55 + + 0 + 0 + 0 + 0 + 255 + + + 56 + + 0 + 0 + 0 + 0 + 255 + + + 57 + + 0 + 0 + 0 + 0 + 255 + + + 58 + + 0 + 0 + 0 + 0 + 255 + + + 59 + + 0 + 0 + 0 + 0 + 255 + + + 60 + + 0 + 0 + 0 + 0 + 255 + + + 61 + + 0 + 0 + 0 + 0 + 255 + + + 62 + + 0 + 0 + 0 + 0 + 255 + + + 63 + + 0 + 0 + 0 + 0 + 255 + + + 64 + + 0 + 0 + 0 + 0 + 255 + + + 65 + + 0 + 0 + 0 + 0 + 255 + + + 66 + + 0 + 0 + 0 + 0 + 255 + + + 67 + + 0 + 0 + 0 + 0 + 255 + + + 68 + + 0 + 0 + 0 + 0 + 255 + + + 69 + + 0 + 0 + 0 + 0 + 255 + + + 70 + + 0 + 0 + 0 + 0 + 255 + + + 71 + Herbaceous + 1198000354 + 223 + 223 + 194 + 255 + + + 72 + + 0 + 0 + 0 + 0 + 255 + + + 73 + + 0 + 0 + 0 + 0 + 255 + + + 74 + + 0 + 0 + 0 + 0 + 255 + + + 75 + + 0 + 0 + 0 + 0 + 255 + + + 76 + + 0 + 0 + 0 + 0 + 255 + + + 77 + + 0 + 0 + 0 + 0 + 255 + + + 78 + + 0 + 0 + 0 + 0 + 255 + + + 79 + + 0 + 0 + 0 + 0 + 255 + + + 80 + + 0 + 0 + 0 + 0 + 255 + + + 81 + Hay/Pasture + 560647664 + 220 + 217 + 57 + 255 + + + 82 + Cultivated Crops + 1464715609 + 171 + 108 + 40 + 255 + + + 83 + + 0 + 0 + 0 + 0 + 255 + + + 84 + + 0 + 0 + 0 + 0 + 255 + + + 85 + + 0 + 0 + 0 + 0 + 255 + + + 86 + + 0 + 0 + 0 + 0 + 255 + + + 87 + + 0 + 0 + 0 + 0 + 255 + + + 88 + + 0 + 0 + 0 + 0 + 255 + + + 89 + + 0 + 0 + 0 + 0 + 255 + + + 90 + Woody Wetlands + 403631293 + 184 + 217 + 235 + 255 + + + 91 + + 0 + 0 + 0 + 0 + 255 + + + 92 + + 0 + 0 + 0 + 0 + 255 + + + 93 + + 0 + 0 + 0 + 0 + 255 + + + 94 + + 0 + 0 + 0 + 0 + 255 + + + 95 + Emergent Herbaceous Wetlands + 137098525 + 108 + 159 + 184 + 255 + + + 96 + + 0 + 96 + 96 + 96 + 255 + + + 97 + + 0 + 97 + 97 + 97 + 255 + + + 98 + + 0 + 98 + 98 + 98 + 255 + + + 99 + + 0 + 99 + 99 + 99 + 255 + + + 100 + + 0 + 100 + 100 + 100 + 255 + + + 101 + + 0 + 101 + 101 + 101 + 255 + + + 102 + + 0 + 102 + 102 + 102 + 255 + + + 103 + + 0 + 103 + 103 + 103 + 255 + + + 104 + + 0 + 104 + 104 + 104 + 255 + + + 105 + + 0 + 105 + 105 + 105 + 255 + + + 106 + + 0 + 106 + 106 + 106 + 255 + + + 107 + + 0 + 107 + 107 + 107 + 255 + + + 108 + + 0 + 108 + 108 + 108 + 255 + + + 109 + + 0 + 109 + 109 + 109 + 255 + + + 110 + + 0 + 110 + 110 + 110 + 255 + + + 111 + + 0 + 111 + 111 + 111 + 255 + + + 112 + + 0 + 112 + 112 + 112 + 255 + + + 113 + + 0 + 113 + 113 + 113 + 255 + + + 114 + + 0 + 114 + 114 + 114 + 255 + + + 115 + + 0 + 115 + 115 + 115 + 255 + + + 116 + + 0 + 116 + 116 + 116 + 255 + + + 117 + + 0 + 117 + 117 + 117 + 255 + + + 118 + + 0 + 118 + 118 + 118 + 255 + + + 119 + + 0 + 119 + 119 + 119 + 255 + + + 120 + + 0 + 120 + 120 + 120 + 255 + + + 121 + + 0 + 121 + 121 + 121 + 255 + + + 122 + + 0 + 122 + 122 + 122 + 255 + + + 123 + + 0 + 123 + 123 + 123 + 255 + + + 124 + + 0 + 124 + 124 + 124 + 255 + + + 125 + + 0 + 125 + 125 + 125 + 255 + + + 126 + + 0 + 126 + 126 + 126 + 255 + + + 127 + + 0 + 127 + 127 + 127 + 255 + + + 128 + + 0 + 128 + 128 + 128 + 255 + + + 129 + + 0 + 129 + 129 + 129 + 255 + + + 130 + + 0 + 130 + 130 + 130 + 255 + + + 131 + + 0 + 131 + 131 + 131 + 255 + + + 132 + + 0 + 132 + 132 + 132 + 255 + + + 133 + + 0 + 133 + 133 + 133 + 255 + + + 134 + + 0 + 134 + 134 + 134 + 255 + + + 135 + + 0 + 135 + 135 + 135 + 255 + + + 136 + + 0 + 136 + 136 + 136 + 255 + + + 137 + + 0 + 137 + 137 + 137 + 255 + + + 138 + + 0 + 138 + 138 + 138 + 255 + + + 139 + + 0 + 139 + 139 + 139 + 255 + + + 140 + + 0 + 140 + 140 + 140 + 255 + + + 141 + + 0 + 141 + 141 + 141 + 255 + + + 142 + + 0 + 142 + 142 + 142 + 255 + + + 143 + + 0 + 143 + 143 + 143 + 255 + + + 144 + + 0 + 144 + 144 + 144 + 255 + + + 145 + + 0 + 145 + 145 + 145 + 255 + + + 146 + + 0 + 146 + 146 + 146 + 255 + + + 147 + + 0 + 147 + 147 + 147 + 255 + + + 148 + + 0 + 148 + 148 + 148 + 255 + + + 149 + + 0 + 149 + 149 + 149 + 255 + + + 150 + + 0 + 150 + 150 + 150 + 255 + + + 151 + + 0 + 151 + 151 + 151 + 255 + + + 152 + + 0 + 152 + 152 + 152 + 255 + + + 153 + + 0 + 153 + 153 + 153 + 255 + + + 154 + + 0 + 154 + 154 + 154 + 255 + + + 155 + + 0 + 155 + 155 + 155 + 255 + + + 156 + + 0 + 156 + 156 + 156 + 255 + + + 157 + + 0 + 157 + 157 + 157 + 255 + + + 158 + + 0 + 158 + 158 + 158 + 255 + + + 159 + + 0 + 159 + 159 + 159 + 255 + + + 160 + + 0 + 160 + 160 + 160 + 255 + + + 161 + + 0 + 161 + 161 + 161 + 255 + + + 162 + + 0 + 162 + 162 + 162 + 255 + + + 163 + + 0 + 163 + 163 + 163 + 255 + + + 164 + + 0 + 164 + 164 + 164 + 255 + + + 165 + + 0 + 165 + 165 + 165 + 255 + + + 166 + + 0 + 166 + 166 + 166 + 255 + + + 167 + + 0 + 167 + 167 + 167 + 255 + + + 168 + + 0 + 168 + 168 + 168 + 255 + + + 169 + + 0 + 169 + 169 + 169 + 255 + + + 170 + + 0 + 170 + 170 + 170 + 255 + + + 171 + + 0 + 171 + 171 + 171 + 255 + + + 172 + + 0 + 172 + 172 + 172 + 255 + + + 173 + + 0 + 173 + 173 + 173 + 255 + + + 174 + + 0 + 174 + 174 + 174 + 255 + + + 175 + + 0 + 175 + 175 + 175 + 255 + + + 176 + + 0 + 176 + 176 + 176 + 255 + + + 177 + + 0 + 177 + 177 + 177 + 255 + + + 178 + + 0 + 178 + 178 + 178 + 255 + + + 179 + + 0 + 179 + 179 + 179 + 255 + + + 180 + + 0 + 180 + 180 + 180 + 255 + + + 181 + + 0 + 181 + 181 + 181 + 255 + + + 182 + + 0 + 182 + 182 + 182 + 255 + + + 183 + + 0 + 183 + 183 + 183 + 255 + + + 184 + + 0 + 184 + 184 + 184 + 255 + + + 185 + + 0 + 185 + 185 + 185 + 255 + + + 186 + + 0 + 186 + 186 + 186 + 255 + + + 187 + + 0 + 187 + 187 + 187 + 255 + + + 188 + + 0 + 188 + 188 + 188 + 255 + + + 189 + + 0 + 189 + 189 + 189 + 255 + + + 190 + + 0 + 190 + 190 + 190 + 255 + + + 191 + + 0 + 191 + 191 + 191 + 255 + + + 192 + + 0 + 192 + 192 + 192 + 255 + + + 193 + + 0 + 193 + 193 + 193 + 255 + + + 194 + + 0 + 194 + 194 + 194 + 255 + + + 195 + + 0 + 195 + 195 + 195 + 255 + + + 196 + + 0 + 196 + 196 + 196 + 255 + + + 197 + + 0 + 197 + 197 + 197 + 255 + + + 198 + + 0 + 198 + 198 + 198 + 255 + + + 199 + + 0 + 199 + 199 + 199 + 255 + + + 200 + + 0 + 200 + 200 + 200 + 255 + + + 201 + + 0 + 201 + 201 + 201 + 255 + + + 202 + + 0 + 202 + 202 + 202 + 255 + + + 203 + + 0 + 203 + 203 + 203 + 255 + + + 204 + + 0 + 204 + 204 + 204 + 255 + + + 205 + + 0 + 205 + 205 + 205 + 255 + + + 206 + + 0 + 206 + 206 + 206 + 255 + + + 207 + + 0 + 207 + 207 + 207 + 255 + + + 208 + + 0 + 208 + 208 + 208 + 255 + + + 209 + + 0 + 209 + 209 + 209 + 255 + + + 210 + + 0 + 210 + 210 + 210 + 255 + + + 211 + + 0 + 211 + 211 + 211 + 255 + + + 212 + + 0 + 212 + 212 + 212 + 255 + + + 213 + + 0 + 213 + 213 + 213 + 255 + + + 214 + + 0 + 214 + 214 + 214 + 255 + + + 215 + + 0 + 215 + 215 + 215 + 255 + + + 216 + + 0 + 216 + 216 + 216 + 255 + + + 217 + + 0 + 217 + 217 + 217 + 255 + + + 218 + + 0 + 218 + 218 + 218 + 255 + + + 219 + + 0 + 219 + 219 + 219 + 255 + + + 220 + + 0 + 220 + 220 + 220 + 255 + + + 221 + + 0 + 221 + 221 + 221 + 255 + + + 222 + + 0 + 222 + 222 + 222 + 255 + + + 223 + + 0 + 223 + 223 + 223 + 255 + + + 224 + + 0 + 224 + 224 + 224 + 255 + + + 225 + + 0 + 225 + 225 + 225 + 255 + + + 226 + + 0 + 226 + 226 + 226 + 255 + + + 227 + + 0 + 227 + 227 + 227 + 255 + + + 228 + + 0 + 228 + 228 + 228 + 255 + + + 229 + + 0 + 229 + 229 + 229 + 255 + + + 230 + + 0 + 230 + 230 + 230 + 255 + + + 231 + + 0 + 231 + 231 + 231 + 255 + + + 232 + + 0 + 232 + 232 + 232 + 255 + + + 233 + + 0 + 233 + 233 + 233 + 255 + + + 234 + + 0 + 234 + 234 + 234 + 255 + + + 235 + + 0 + 235 + 235 + 235 + 255 + + + 236 + + 0 + 236 + 236 + 236 + 255 + + + 237 + + 0 + 237 + 237 + 237 + 255 + + + 238 + + 0 + 238 + 238 + 238 + 255 + + + 239 + + 0 + 239 + 239 + 239 + 255 + + + 240 + + 0 + 240 + 240 + 240 + 255 + + + 241 + + 0 + 241 + 241 + 241 + 255 + + + 242 + + 0 + 242 + 242 + 242 + 255 + + + 243 + + 0 + 243 + 243 + 243 + 255 + + + 244 + + 0 + 244 + 244 + 244 + 255 + + + 245 + + 0 + 245 + 245 + 245 + 255 + + + 246 + + 0 + 246 + 246 + 246 + 255 + + + 247 + + 0 + 247 + 247 + 247 + 255 + + + 248 + + 0 + 248 + 248 + 248 + 255 + + + 249 + + 0 + 249 + 249 + 249 + 255 + + + 250 + + 0 + 250 + 250 + 250 + 255 + + + 251 + + 0 + 251 + 251 + 251 + 255 + + + 252 + + 0 + 252 + 252 + 252 + 255 + + + 253 + + 0 + 253 + 253 + 253 + 255 + + + 254 + + 0 + 254 + 254 + 254 + 255 + + + 255 + + 0 + 255 + 255 + 255 + 255 + + + + diff --git a/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif b/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif new file mode 100644 index 00000000..0fcb38b5 Binary files /dev/null and b/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif differ diff --git a/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif.aux.xml b/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif.aux.xml new file mode 100644 index 00000000..b25aa006 --- /dev/null +++ b/tests/testdata/nlcd_2021_land_cover_l48_20230630.tif.aux.xml @@ -0,0 +1,2346 @@ + + + NLCD Land Cover Class + + + value + 0 + 0 + + + NLCD Land Cover Class + 2 + 0 + + + Histogram + 1 + 0 + + + Red + 0 + 0 + + + Green + 0 + 0 + + + Blue + 0 + 0 + + + Opacity + 0 + 0 + + + 0 + Unclassified + 7853863229 + 0 + 0 + 0 + 0 + + + 1 + + 0 + 0 + 0 + 0 + 255 + + + 2 + + 0 + 0 + 0 + 0 + 255 + + + 3 + + 0 + 0 + 0 + 0 + 255 + + + 4 + + 0 + 0 + 0 + 0 + 255 + + + 5 + + 0 + 0 + 0 + 0 + 255 + + + 6 + + 0 + 0 + 0 + 0 + 255 + + + 7 + + 0 + 0 + 0 + 0 + 255 + + + 8 + + 0 + 0 + 0 + 0 + 255 + + + 9 + + 0 + 0 + 0 + 0 + 255 + + + 10 + + 0 + 0 + 0 + 0 + 255 + + + 11 + Open Water + 466021535 + 70 + 107 + 159 + 255 + + + 12 + Perennial Snow/Ice + 959819 + 209 + 222 + 248 + 255 + + + 13 + + 0 + 0 + 0 + 0 + 255 + + + 14 + + 0 + 0 + 0 + 0 + 255 + + + 15 + + 0 + 0 + 0 + 0 + 255 + + + 16 + + 0 + 0 + 0 + 0 + 255 + + + 17 + + 0 + 0 + 0 + 0 + 255 + + + 18 + + 0 + 0 + 0 + 0 + 255 + + + 19 + + 0 + 0 + 0 + 0 + 255 + + + 20 + + 0 + 0 + 0 + 0 + 255 + + + 21 + Developed, Open Space + 240377730 + 222 + 197 + 197 + 255 + + + 22 + Developed, Low Intensity + 154252969 + 217 + 146 + 130 + 255 + + + 23 + Developed, Medium Intensity + 94385151 + 235 + 0 + 0 + 255 + + + 24 + Developed, High Intensity + 33594478 + 171 + 0 + 0 + 255 + + + 25 + + 0 + 0 + 0 + 0 + 255 + + + 26 + + 0 + 0 + 0 + 0 + 255 + + + 27 + + 0 + 0 + 0 + 0 + 255 + + + 28 + + 0 + 0 + 0 + 0 + 255 + + + 29 + + 0 + 0 + 0 + 0 + 255 + + + 30 + + 0 + 0 + 0 + 0 + 255 + + + 31 + Barren Land + 88976651 + 179 + 172 + 159 + 255 + + + 32 + + 0 + 0 + 0 + 0 + 255 + + + 33 + + 0 + 0 + 0 + 0 + 255 + + + 34 + + 0 + 0 + 0 + 0 + 255 + + + 35 + + 0 + 0 + 0 + 0 + 255 + + + 36 + + 0 + 0 + 0 + 0 + 255 + + + 37 + + 0 + 0 + 0 + 0 + 255 + + + 38 + + 0 + 0 + 0 + 0 + 255 + + + 39 + + 0 + 0 + 0 + 0 + 255 + + + 40 + + 0 + 0 + 0 + 0 + 255 + + + 41 + Deciduous Forest + 832916302 + 104 + 171 + 95 + 255 + + + 42 + Evergreen Forest + 1027738936 + 28 + 95 + 44 + 255 + + + 43 + Mixed Forest + 304954315 + 181 + 197 + 143 + 255 + + + 44 + + 0 + 0 + 0 + 0 + 255 + + + 45 + + 0 + 0 + 0 + 0 + 255 + + + 46 + + 0 + 0 + 0 + 0 + 255 + + + 47 + + 0 + 0 + 0 + 0 + 255 + + + 48 + + 0 + 0 + 0 + 0 + 255 + + + 49 + + 0 + 0 + 0 + 0 + 255 + + + 50 + + 0 + 0 + 0 + 0 + 255 + + + 51 + + 0 + 0 + 0 + 0 + 255 + + + 52 + Shrub/Scrub + 1953262413 + 204 + 184 + 121 + 255 + + + 53 + + 0 + 0 + 0 + 0 + 255 + + + 54 + + 0 + 0 + 0 + 0 + 255 + + + 55 + + 0 + 0 + 0 + 0 + 255 + + + 56 + + 0 + 0 + 0 + 0 + 255 + + + 57 + + 0 + 0 + 0 + 0 + 255 + + + 58 + + 0 + 0 + 0 + 0 + 255 + + + 59 + + 0 + 0 + 0 + 0 + 255 + + + 60 + + 0 + 0 + 0 + 0 + 255 + + + 61 + + 0 + 0 + 0 + 0 + 255 + + + 62 + + 0 + 0 + 0 + 0 + 255 + + + 63 + + 0 + 0 + 0 + 0 + 255 + + + 64 + + 0 + 0 + 0 + 0 + 255 + + + 65 + + 0 + 0 + 0 + 0 + 255 + + + 66 + + 0 + 0 + 0 + 0 + 255 + + + 67 + + 0 + 0 + 0 + 0 + 255 + + + 68 + + 0 + 0 + 0 + 0 + 255 + + + 69 + + 0 + 0 + 0 + 0 + 255 + + + 70 + + 0 + 0 + 0 + 0 + 255 + + + 71 + Herbaceous + 1223488270 + 223 + 223 + 194 + 255 + + + 72 + + 0 + 0 + 0 + 0 + 255 + + + 73 + + 0 + 0 + 0 + 0 + 255 + + + 74 + + 0 + 0 + 0 + 0 + 255 + + + 75 + + 0 + 0 + 0 + 0 + 255 + + + 76 + + 0 + 0 + 0 + 0 + 255 + + + 77 + + 0 + 0 + 0 + 0 + 255 + + + 78 + + 0 + 0 + 0 + 0 + 255 + + + 79 + + 0 + 0 + 0 + 0 + 255 + + + 80 + + 0 + 0 + 0 + 0 + 255 + + + 81 + Hay/Pasture + 558088764 + 220 + 217 + 57 + 255 + + + 82 + Cultivated Crops + 1453539695 + 171 + 108 + 40 + 255 + + + 83 + + 0 + 0 + 0 + 0 + 255 + + + 84 + + 0 + 0 + 0 + 0 + 255 + + + 85 + + 0 + 0 + 0 + 0 + 255 + + + 86 + + 0 + 0 + 0 + 0 + 255 + + + 87 + + 0 + 0 + 0 + 0 + 255 + + + 88 + + 0 + 0 + 0 + 0 + 255 + + + 89 + + 0 + 0 + 0 + 0 + 255 + + + 90 + Woody Wetlands + 403110859 + 184 + 217 + 235 + 255 + + + 91 + + 0 + 0 + 0 + 0 + 255 + + + 92 + + 0 + 0 + 0 + 0 + 255 + + + 93 + + 0 + 0 + 0 + 0 + 255 + + + 94 + + 0 + 0 + 0 + 0 + 255 + + + 95 + Emergent Herbaceous Wetlands + 142573444 + 108 + 159 + 184 + 255 + + + 96 + + 0 + 96 + 96 + 96 + 255 + + + 97 + + 0 + 97 + 97 + 97 + 255 + + + 98 + + 0 + 98 + 98 + 98 + 255 + + + 99 + + 0 + 99 + 99 + 99 + 255 + + + 100 + + 0 + 100 + 100 + 100 + 255 + + + 101 + + 0 + 101 + 101 + 101 + 255 + + + 102 + + 0 + 102 + 102 + 102 + 255 + + + 103 + + 0 + 103 + 103 + 103 + 255 + + + 104 + + 0 + 104 + 104 + 104 + 255 + + + 105 + + 0 + 105 + 105 + 105 + 255 + + + 106 + + 0 + 106 + 106 + 106 + 255 + + + 107 + + 0 + 107 + 107 + 107 + 255 + + + 108 + + 0 + 108 + 108 + 108 + 255 + + + 109 + + 0 + 109 + 109 + 109 + 255 + + + 110 + + 0 + 110 + 110 + 110 + 255 + + + 111 + + 0 + 111 + 111 + 111 + 255 + + + 112 + + 0 + 112 + 112 + 112 + 255 + + + 113 + + 0 + 113 + 113 + 113 + 255 + + + 114 + + 0 + 114 + 114 + 114 + 255 + + + 115 + + 0 + 115 + 115 + 115 + 255 + + + 116 + + 0 + 116 + 116 + 116 + 255 + + + 117 + + 0 + 117 + 117 + 117 + 255 + + + 118 + + 0 + 118 + 118 + 118 + 255 + + + 119 + + 0 + 119 + 119 + 119 + 255 + + + 120 + + 0 + 120 + 120 + 120 + 255 + + + 121 + + 0 + 121 + 121 + 121 + 255 + + + 122 + + 0 + 122 + 122 + 122 + 255 + + + 123 + + 0 + 123 + 123 + 123 + 255 + + + 124 + + 0 + 124 + 124 + 124 + 255 + + + 125 + + 0 + 125 + 125 + 125 + 255 + + + 126 + + 0 + 126 + 126 + 126 + 255 + + + 127 + + 0 + 127 + 127 + 127 + 255 + + + 128 + + 0 + 128 + 128 + 128 + 255 + + + 129 + + 0 + 129 + 129 + 129 + 255 + + + 130 + + 0 + 130 + 130 + 130 + 255 + + + 131 + + 0 + 131 + 131 + 131 + 255 + + + 132 + + 0 + 132 + 132 + 132 + 255 + + + 133 + + 0 + 133 + 133 + 133 + 255 + + + 134 + + 0 + 134 + 134 + 134 + 255 + + + 135 + + 0 + 135 + 135 + 135 + 255 + + + 136 + + 0 + 136 + 136 + 136 + 255 + + + 137 + + 0 + 137 + 137 + 137 + 255 + + + 138 + + 0 + 138 + 138 + 138 + 255 + + + 139 + + 0 + 139 + 139 + 139 + 255 + + + 140 + + 0 + 140 + 140 + 140 + 255 + + + 141 + + 0 + 141 + 141 + 141 + 255 + + + 142 + + 0 + 142 + 142 + 142 + 255 + + + 143 + + 0 + 143 + 143 + 143 + 255 + + + 144 + + 0 + 144 + 144 + 144 + 255 + + + 145 + + 0 + 145 + 145 + 145 + 255 + + + 146 + + 0 + 146 + 146 + 146 + 255 + + + 147 + + 0 + 147 + 147 + 147 + 255 + + + 148 + + 0 + 148 + 148 + 148 + 255 + + + 149 + + 0 + 149 + 149 + 149 + 255 + + + 150 + + 0 + 150 + 150 + 150 + 255 + + + 151 + + 0 + 151 + 151 + 151 + 255 + + + 152 + + 0 + 152 + 152 + 152 + 255 + + + 153 + + 0 + 153 + 153 + 153 + 255 + + + 154 + + 0 + 154 + 154 + 154 + 255 + + + 155 + + 0 + 155 + 155 + 155 + 255 + + + 156 + + 0 + 156 + 156 + 156 + 255 + + + 157 + + 0 + 157 + 157 + 157 + 255 + + + 158 + + 0 + 158 + 158 + 158 + 255 + + + 159 + + 0 + 159 + 159 + 159 + 255 + + + 160 + + 0 + 160 + 160 + 160 + 255 + + + 161 + + 0 + 161 + 161 + 161 + 255 + + + 162 + + 0 + 162 + 162 + 162 + 255 + + + 163 + + 0 + 163 + 163 + 163 + 255 + + + 164 + + 0 + 164 + 164 + 164 + 255 + + + 165 + + 0 + 165 + 165 + 165 + 255 + + + 166 + + 0 + 166 + 166 + 166 + 255 + + + 167 + + 0 + 167 + 167 + 167 + 255 + + + 168 + + 0 + 168 + 168 + 168 + 255 + + + 169 + + 0 + 169 + 169 + 169 + 255 + + + 170 + + 0 + 170 + 170 + 170 + 255 + + + 171 + + 0 + 171 + 171 + 171 + 255 + + + 172 + + 0 + 172 + 172 + 172 + 255 + + + 173 + + 0 + 173 + 173 + 173 + 255 + + + 174 + + 0 + 174 + 174 + 174 + 255 + + + 175 + + 0 + 175 + 175 + 175 + 255 + + + 176 + + 0 + 176 + 176 + 176 + 255 + + + 177 + + 0 + 177 + 177 + 177 + 255 + + + 178 + + 0 + 178 + 178 + 178 + 255 + + + 179 + + 0 + 179 + 179 + 179 + 255 + + + 180 + + 0 + 180 + 180 + 180 + 255 + + + 181 + + 0 + 181 + 181 + 181 + 255 + + + 182 + + 0 + 182 + 182 + 182 + 255 + + + 183 + + 0 + 183 + 183 + 183 + 255 + + + 184 + + 0 + 184 + 184 + 184 + 255 + + + 185 + + 0 + 185 + 185 + 185 + 255 + + + 186 + + 0 + 186 + 186 + 186 + 255 + + + 187 + + 0 + 187 + 187 + 187 + 255 + + + 188 + + 0 + 188 + 188 + 188 + 255 + + + 189 + + 0 + 189 + 189 + 189 + 255 + + + 190 + + 0 + 190 + 190 + 190 + 255 + + + 191 + + 0 + 191 + 191 + 191 + 255 + + + 192 + + 0 + 192 + 192 + 192 + 255 + + + 193 + + 0 + 193 + 193 + 193 + 255 + + + 194 + + 0 + 194 + 194 + 194 + 255 + + + 195 + + 0 + 195 + 195 + 195 + 255 + + + 196 + + 0 + 196 + 196 + 196 + 255 + + + 197 + + 0 + 197 + 197 + 197 + 255 + + + 198 + + 0 + 198 + 198 + 198 + 255 + + + 199 + + 0 + 199 + 199 + 199 + 255 + + + 200 + + 0 + 200 + 200 + 200 + 255 + + + 201 + + 0 + 201 + 201 + 201 + 255 + + + 202 + + 0 + 202 + 202 + 202 + 255 + + + 203 + + 0 + 203 + 203 + 203 + 255 + + + 204 + + 0 + 204 + 204 + 204 + 255 + + + 205 + + 0 + 205 + 205 + 205 + 255 + + + 206 + + 0 + 206 + 206 + 206 + 255 + + + 207 + + 0 + 207 + 207 + 207 + 255 + + + 208 + + 0 + 208 + 208 + 208 + 255 + + + 209 + + 0 + 209 + 209 + 209 + 255 + + + 210 + + 0 + 210 + 210 + 210 + 255 + + + 211 + + 0 + 211 + 211 + 211 + 255 + + + 212 + + 0 + 212 + 212 + 212 + 255 + + + 213 + + 0 + 213 + 213 + 213 + 255 + + + 214 + + 0 + 214 + 214 + 214 + 255 + + + 215 + + 0 + 215 + 215 + 215 + 255 + + + 216 + + 0 + 216 + 216 + 216 + 255 + + + 217 + + 0 + 217 + 217 + 217 + 255 + + + 218 + + 0 + 218 + 218 + 218 + 255 + + + 219 + + 0 + 219 + 219 + 219 + 255 + + + 220 + + 0 + 220 + 220 + 220 + 255 + + + 221 + + 0 + 221 + 221 + 221 + 255 + + + 222 + + 0 + 222 + 222 + 222 + 255 + + + 223 + + 0 + 223 + 223 + 223 + 255 + + + 224 + + 0 + 224 + 224 + 224 + 255 + + + 225 + + 0 + 225 + 225 + 225 + 255 + + + 226 + + 0 + 226 + 226 + 226 + 255 + + + 227 + + 0 + 227 + 227 + 227 + 255 + + + 228 + + 0 + 228 + 228 + 228 + 255 + + + 229 + + 0 + 229 + 229 + 229 + 255 + + + 230 + + 0 + 230 + 230 + 230 + 255 + + + 231 + + 0 + 231 + 231 + 231 + 255 + + + 232 + + 0 + 232 + 232 + 232 + 255 + + + 233 + + 0 + 233 + 233 + 233 + 255 + + + 234 + + 0 + 234 + 234 + 234 + 255 + + + 235 + + 0 + 235 + 235 + 235 + 255 + + + 236 + + 0 + 236 + 236 + 236 + 255 + + + 237 + + 0 + 237 + 237 + 237 + 255 + + + 238 + + 0 + 238 + 238 + 238 + 255 + + + 239 + + 0 + 239 + 239 + 239 + 255 + + + 240 + + 0 + 240 + 240 + 240 + 255 + + + 241 + + 0 + 241 + 241 + 241 + 255 + + + 242 + + 0 + 242 + 242 + 242 + 255 + + + 243 + + 0 + 243 + 243 + 243 + 255 + + + 244 + + 0 + 244 + 244 + 244 + 255 + + + 245 + + 0 + 245 + 245 + 245 + 255 + + + 246 + + 0 + 246 + 246 + 246 + 255 + + + 247 + + 0 + 247 + 247 + 247 + 255 + + + 248 + + 0 + 248 + 248 + 248 + 255 + + + 249 + + 0 + 249 + 249 + 249 + 255 + + + 250 + + 0 + 250 + 250 + 250 + 255 + + + 251 + + 0 + 251 + 251 + 251 + 255 + + + 252 + + 0 + 252 + 252 + 252 + 255 + + + 253 + + 0 + 253 + 253 + 253 + 255 + + + 254 + + 0 + 254 + 254 + 254 + 255 + + + 255 + + 0 + 255 + 255 + 255 + 255 + + + + diff --git a/tests/testthat/test-calc-nlcd-ratio.R b/tests/testthat/test-calc-nlcd-ratio.R new file mode 100644 index 00000000..a4e2e58e --- /dev/null +++ b/tests/testthat/test-calc-nlcd-ratio.R @@ -0,0 +1,99 @@ +withr::local_package("spData") + +test_that("Check extract_nlcd_ratio works", { + point_us1 <- cbind(lon = -114.7, lat = 38.9, dem = 40) + point_us2 <- cbind(lon = -114, lat = 39, dem = 15) + point_ak <- cbind(lon = -155.997, lat = 69.3884, dem = 100) # alaska + point_fr <- cbind(lon = 2.957, lat = 43.976, dem = 15) # france + eg_data <- rbind(point_us1, point_us2, point_ak, point_fr) %>% + as.data.frame() %>% + vect(., crs = "EPSG:4326") + getwd() + path_testdata <- "../testdata/" + # CHECK INPUT (error message) + # -- buf_radius is numeric + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + buf_radius = "1000", + nlcd_path = path_testdata), + "buf_radius is not a numeric." + ) + # -- buf_radius has likely value + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + buf_radius = -3, + nlcd_path = path_testdata), + "buf_radius has not a likely value." + ) + # -- year is numeric + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + year = "2019", + nlcd_path = path_testdata), + "year is not a numeric." + ) + # -- year has likely value + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + year = 20192, + nlcd_path = path_testdata), + "NLCD data not available for this year." + ) + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + year = 2020, + nlcd_path = path_testdata), + "NLCD data not available for this year." + ) + # -- data_vect is a SpatVector + expect_error( + calc_nlcd_ratio(data_vect = 12, + nlcd_path = path_testdata), + "data_vect is not a terra::SpatVector." + ) + # -- nlcd_path is not a character + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + nlcd_path = 2), + "nlcd_path is not a character." + ) + # -- nlcd_path does not exist + nice_sentence <- "That's one small step for a man, a giant leap for mankind." + expect_error( + calc_nlcd_ratio(data_vect = eg_data, + nlcd_path = nice_sentence), + "nlcd_path does not exist." + ) + + # CHECK OUTPUT + year <- 2021 + buf_radius <- 3000 + expect_no_error(calc_nlcd_ratio( + data_vect = eg_data, + year = year, + buf_radius = buf_radius, + nlcd_path = path_testdata + )) + output <- calc_nlcd_ratio( + data_vect = eg_data, + year = year, + buf_radius = buf_radius, + nlcd_path = path_testdata + ) + # -- returns a SpatVector + expect_equal(class(output)[1], "SpatVector") + # -- crs is the same than input + expect_true(terra::same.crs(eg_data, output)) + # -- out-of-mainland-US points removed (France and Alaska) + expect_equal(nrow(output), 2) + # -- initial names are still in the output SpatVector + expect_true(all(names(eg_data) %in% names(output))) + # -- check the value of some of the points in the US + expect_equal(output$frac_EFO_2021_3000m[1], 0.7940682, tolerance = 1e-7) + expect_equal(output$frac_SHB_2021_3000m[2], 0.9987249, tolerance = 1e-7) + # -- class fraction rows should sum to 1 + expect_equal(rowSums(as.data.frame(output[, 2:ncol(output)])), + rep(1, 2), + tolerance = 1e-7 + ) +}) diff --git a/vignettes/calc_nlcd_ratio_monitors.Rmd b/vignettes/calc_nlcd_ratio_monitors.Rmd new file mode 100644 index 00000000..f71202fd --- /dev/null +++ b/vignettes/calc_nlcd_ratio_monitors.Rmd @@ -0,0 +1,157 @@ +--- +title: "NLCD ratio: calculation at monitor locations" +author: "Eva Marques" +output: rmarkdown::html_vignette +date: "2023-11-27" +vignette: > + %\VignetteIndexEntry{NLCD ratio: calculation at monitor locations} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, message=F} +pkgs <- c( + "terra", + "exactextractr", + "tidyverse", + "data.table", + "sf", + "spData", + "ggplot2", + "scatterpie" +) +sapply(pkgs, library, character.only = TRUE) +``` + + +## Open NLCD data + +```{r, eval = F} +nlcd_2019 <- rast("/Volumes/set/NLCD/nlcd_2019_land_cover_l48_20210604.tif") +nlcd_2021 <- rast("/Volumes/set/NLCD/nlcd_2021_land_cover_l48_20230630.tif") + +# check that unit is in meters +linearUnits(nlcd_2021) +same.crs(nlcd_2019, nlcd_2021) +``` + +Test to see if data is well loaded + +```{r, eval = F} +lat <- c(1900000, 2000000, 2000000, 1900000) +lon <- c(-1600000, -1600000, -1500000, -1500000) +ext <- vect(cbind(lon, lat), type = "polygons", crs = crs(nlcd_2021)) +nlcd_small <- terra::crop(nlcd_2021, ext, mask = TRUE) +plot(nlcd_2021) +plot(ext, add = TRUE, col = "black") +plot(nlcd_small) + +point1 <- cbind(-2e6, 3e6) +terra::extract(nlcd_2021, point1) + +point2 <- cbind(0, 1500000) +terra::extract(nlcd_2021, point2) +``` + + +## Open monitors data + +```{r} +path_aqs <- "/Volumes/set/Projects/NRT-AP-Model/input/aqs/" +aqs_mon <- fread(paste0(path_aqs, "aqs_monitors.csv")) +aqs_mon <- aqs_mon[`Parameter Code` == 88101, ] # PM2.5 monitors only +``` + +Process monitors + +```{r} +aqs_mon <- aqs_mon %>% + rename(lon = Longitude) %>% + rename(lat = Latitude) + +source("/Volumes/set/Projects/NRT-AP-Model/R/manipulate_spacetime_data.R") +a <- project_dt(aqs_mon[Datum == "NAD83"], "EPSG:4269", "EPSG:4326") +a$lon_ori <- NULL +a$lat_ori <- NULL +b <- aqs_mon[Datum != "NAD83"] +aqs_mon <- rbind(a, b) %>% + terra::vect(., + geom = c("lon", "lat"), + crs = "EPSG:4326", + keepgeom = FALSE + ) +``` + + +## NLCD classes + +```{r} +nlcd_classes <- list(value = c(0, 11, 21, 22, 23, 24, 31, 41, 42, 43, 52, + 71, 81, 82, 90, 95), + class = c("Unc", "WTR", "OSD", "LID", "MID", "HID", + "BRN", "DFO", "EFO", "MFO", "SHB", + "GRS", "PAS", "CRP", "WDW", "EHW"), + names = c("Unclassified", + "Open Water", + "Developed, Open Space", + "Developed, Low Intensity", + "Developed, Medium Intensity", + "Developed, High Intensity", + "Barren Land", + "Deciduous Forest", + "Evergreen Forest", + "Mixed Forest", + "Shrub/Scrub", + "Herbaceous", + "Hay/Pasture", + "Cultivated Crops", + "Woody Wetlands", + "Emergent Herbaceous Wetlands"), + col = c("white", "#476ba1", "#decaca", "#d99482", + "#ee0000", "#ab0000", "#b3aea3", "#68ab63", + "#1c6330", "#b5ca8f", "#ccba7d", "#e3e3c2", + "#dcd93d", "#ab7028", "#bad9eb", "#70a3ba")) +nlcd_classes <- as.data.frame(nlcd_classes) + +nlcd_classes + +# write in ./inst/extdata +write.csv(nlcd_classes, "../inst/extdata/nlcd_classes.csv") +``` + + +## Compute NLCD ratio in 100m buffers + +```{r} +source("../R/calc_nlcd_ratio.R") +start <- Sys.time() +aqs_mon_nlcd <- calc_nlcd_ratio( + data_vect = aqs_mon, + buf_radius = 100, + year = 2021, + nlcd_path = "/Volumes/set/NLCD/" +) +end <- Sys.time() +end - start +``` + +Plot pie chart + +```{r, fig.width = 12} +map <- map_data("county") +df <- as.data.frame(aqs_mon_nlcd, geom = "XY") +names(df)[names(df) == "x"] <- "long" +names(df)[names(df) == "y"] <- "lat" +frac <- df[names(df)[grepl("frac_", names(df))]] +df <- cbind(df[, c("long", "lat")], frac) + +g <- ggplot() + + geom_polygon(data = map, aes(x = long, y = lat, group = group)) + + geom_scatterpie(aes(x = long, y = lat, r = .5), + data = df, + cols = names(df)[3:ncol(df)], color = NA + ) + + scale_fill_manual(values = nlcd_classes$col, labels = nlcd_classes$names) + + coord_equal() +g +```