Get the rainfall data
Now that we have objects for the states we can create a raster grid to represent the 0.5 x 0.5 degree grid that is the NASA-POWER data and select only cells that fall within the two states of interest.
POWER source
Create a grid of 0.5 x 0.5 arc degrees and extract the x, y values from it for each state to use the coordinates to query the POWER data.
# create a global 0.5 x 0.5 degree raster object
r <- raster(xmn = -180,
xmx = 180,
ymn = -90,
ymx = 90,
resolution = 0.5)
values(r) <- 1:ncell(r)
plot(r, main = "Full global raster at 0.5 x 0.5 degrees")
# Extract the two states, first crop by bounding box, then mask the raster
PR_coords <- crop(r, PR)
RS_coords <- crop(r, RS)
PR_coords <- mask(PR_coords, PR)
plot(PR_coords, main = "Paraná")
# add the sf object, note the use of indexing to plot only the first column
# to avoid warning messages
plot(PR[, 1], col = NA, add = TRUE)
RS_coords <- mask(RS_coords, RS)
plot(RS_coords, main = "Rio Grande do Sul")
plot(RS[, 1], col = NA, add = TRUE)
# extract the centroid values of the cells to use querying the POWER data
PR_coords <- as.data.frame(xyFromCell(PR_coords, 1:ncell(PR_coords)))
RS_coords <- as.data.frame(xyFromCell(RS_coords, 1:ncell(RS_coords)))
names(PR_coords) <- names(RS_coords) <- c("LON", "LAT")
coords <- rbind(PR_coords, RS_coords)WARNING This step is time intensive. WARNING
Using nested for() loops, query the NASA-POWER database to gather precipitation data for the states where rust was reported and save a CSV file of the rainfall.
power <- vector(mode = "list", 2) # hold two growing seasons
precip <- vector(mode = "list", nrow(coords)) # hold the cells
seasons <- list(
c("2014-11-01", "2015-03-31"),
c("2015-11-01", "2016-03-31")
)
for (i in seq_along(seasons)) {
# two "seasons" (outer loop 2x)
season <- seasons[[i]]
# inner loop for each pair coords
for (j in seq_along(1:nrow(coords))) {
NA_df <-
data.frame(
LON = coords[1, 1],
LAT = coords[1, 2],
YEAR = NA,
MM = NA,
DD = NA,
DOY = NA,
YYYYMMDD = NA,
PRECTOT = NA
)
p_get_power <- possibly(get_power, otherwise = NA_df)
# 312 coordinate pairs (inner loop 312x)
site <- as.numeric(coords[j,])
power_precip <- p_get_power(
community = "AG",
lonlat = site,
pars = "PRECTOT",
dates = season,
temporal_average = "DAILY"
)
precip[[j]] <- power_precip
Sys.sleep(5) # wait 5 seconds between requests so we don't hammer the server
}
precip_df <- bind_rows(precip)
power[[i]] <- precip_df
}
power_df <- bind_rows(power)
# adds states to rows
power_df <- data.frame(STATE = c(rep("PR", nrow(PR_coords)),
rep("RS", nrow(RS_coords))),
power_df)