This is an example of using the spatial.R code to aggregate gridded weather data. Specifically, we aggregate gridded maximum temperature data to Brazilian municipalities.
First, load the required packages. Note that this library, research-common, is in my local folder ~/projects/research-common, but it could be anywhere.
In [16]:
library(PBSmapping)
library(ncdf4)
source("~/projects/research-common/R/distance.R")
source("~/projects/research-common/R/spatial.R", chdir=T)
Load the shapefile. I use PBSmapping, which represents shapefiles as lists of points.
In [28]:
shape <- importShapefile("~/data/political/brazil/BRA_adm2/BRA_adm2")
Load the gridded data.
In [8]:
ncdf <- nc_open("~/data/agmerra/monthly/tasmax_agmerra_1980-2010.mm.nc4")
longitude <- ncvar_get(ncdf, "lon")
latitude <- ncvar_get(ncdf, "lat")
The prepareEvents function identifies the center of each grid cell, which will be used to make the averages. In a loop over multiple files, it's best to only do this once.
In [10]:
events <- prepareEvents(longitude, latitude, attr(shape, "projection"))
Get out the weather data itself. It's important that this has dimensions of the order longitude, latitude, time.
In [25]:
lonlattimeraster <- ncvar_get(ncdf, "tasmax")
Also get out the time variable, for producing output.
In [14]:
time <- ncvar_get(ncdf, "time")
Iterate through all of the polygons, constructing averages and adding them to the result dataframe.
In [26]:
values <- data.frame(PID=c(), time=c(), average=c())
for (pid in 1:10) {
print(pid)
average <-
spaceTimeRasterAverage(events, lonlattimeraster,
subset(shape, PID == pid))
values <- rbind(values, data.frame(PID=pid, time, average))
}
Here's what the final result looks like.
In [27]:
values