Extract mean values of climate variables by level 2 administrative area country by country using R?

[Working with R version 3.2.2 in a Mac computer.]

I am trying to calculate mean values for climatic variables (http://www.worldclim.org/) for spatial objects in R. The challenge is that I'm trying to calculate these means for every level 2 administrative area in the world (http://www.gadm.org/), and I need an efficient way of calculating the statistics.

At the moment, I am trying to get this information country by country. Notice that in order to save time, I am just trying to do it for Austria (later on I will create a loop to do it for every country). Before starting I had to download the SRTM data and part of the script from here: https://github.com/sikli/srtm_country.

I think I miss some command to tell R to join tiles, boundaries, and climate data when I try to extract mean values for the climate variables.

In what follows, I present the procedure I am following.

``````#Clear all
rm(list = ls())

library(raster)
library(rgeos)
library(rasterVis)

setwd('/.../srtm_country-master')

#------------SETTINGS--------------

country_name <- "AUT"                             #Austria
shp          <- shapefile("/.../srtm_country-master/srtm/tiles.shp")       #Path to SRTM Tiles (can be found in subfolder srtm)

#------------EXECUTE FROM HERE--------------

#Get country geometry first
country = country_name,
level=2)

#Intersect country geometry with tile grid
intersects <- gIntersects(country, shp, byid=T)
tiles      <- shp[intersects[,1],]

srtm_list  <- list()
for(i in 1:length(tiles)) {
lon <- extent(tiles[i,])[1]  + (extent(tiles[i,])[2] - extent(tiles[i,])[1]) / 2
lat <- extent(tiles[i,])[3]  + (extent(tiles[i,])[4] - extent(tiles[i,])[3]) / 2

tile <- getData('worldclim', var='bio', res=10,
lon=lon,
lat=lat)

srtm_list[[i]] <- tile
}
``````

THIS IS WHAT IT IS NOT WORKING -->

``````#Create a list with all exported .bil files from the climate data folder
ras_lst <- list.files("/.../wc10/",full.names=TRUE, pattern=".bil\$")

climatemosaic <- stack(ras_lst)

#Extract mean values
AFGmeans <- extract(climatemosaic, country, fun=mean, na.rm=TRUE, small=TRUE, layer=1, nl=19, sp=TRUE)
AFGmeansdf <- as.data.frame(AFGmeans)

setwd("/.../srtm_country-master/Exports")
write.csv(AFGmeansdf, file= "AUTmeansdf.csv")
``````
• How is it not working? Are you getting error messages? Zeroes? – Spacedman Feb 6 '17 at 12:47

Unfortunately, you haven't included the error message associated with the above code, so it's not quite straightforward to guess what exactly went wrong. Anyway, when working with `res = 10` (i.e., spatial resolution equals 10 minutes of a degree), such an iterative approach to dowload WorldClim data for each tile separately is not required. As described in `?getData`,

In the case of res=0.5, you must also provide a lon and lat argument for a tile; for the lower resolutions global data will be downloaded.

Consequently, your code can be reduced as follows (which works on my machine):

``````library(raster)

## get country geometry
if (!dir.exists("data")) dir.create("data")
aut <- getData("GADM", country = "AUT", level = 2, path = "data")

path <- "data/worldclim"
if (!dir.exists(path)) dir.create(path)
bio <- getData("worldclim", var = "bio", res = 10, path = path)

## extract values
val <- extract(bio, aut, fun = mean, na.rm = TRUE, sp = TRUE)

## write to file
dat <- data.frame(val)
write.csv(dat, file = "AUTmeansdf.csv", quote = FALSE, row.names = FALSE)
``````

Unfortunately, the velox package, which provides a quite fast extraction method (see `?VeloxRaster_extract`) and could hence come in handy for your purposes, is not compatible with `na.rm = TRUE` (yet). I already issued a feature request on GitHub, but this will likely take some time until available. If you should be interested anyway, have a look at the official README which describes the process in detail. Using `velox`, the code (section "write to file" excluded) would look as follows:

``````library(velox)
vlx <- velox(bio)
xtr <- vlx\$extract(aut, fun = mean)
colnames(xtr) <- names(bio)
val <- cbind(aut, xtr)
``````

Update for 30 arc-seconds data:

In order to achieve a spatial resolution of 30 seconds (i.e., `res = 0.5`) and stick with this country-based approach, I would suggest the following:

First, download data for the entire WorldClim tile grid used for 30 arc-seconds resolution manually. If your goal is to run this analysis globally, you'll need it anyway.

``````## download worldclim 'bio' variables globally and extract
lst <- lapply(c("bio1-9", "bio10-19"), function(i) {
url <- paste0("http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/",
i, "_30s_bil.zip")
destfile <- paste0("data/wc0.5/", basename(url))
unzip(destfile, exdir = dirname(destfile))
})
``````

Next, import the files into R. Since `unzip` returns an invisible vector of extracted file names that can easily be saved as an object, this is straightforward.

``````## import files into r
fls <- unlist(lst)
fls <- fls[grep(".bil\$", fls)]
rst <- stack(fls)
``````

Finally, loop over all ISO3 codes of the countries you want to process (see `getData('ISO3')`) and extract the corresponding mean values per level-2 administrative area. Since this might take some time, I recommend to enable parallel processing e.g. through parallel and `parLapply` (or, alternatively, foreach and `%dopar%`).

``````## enable parallel processing and load/export relevant packages/objects
library(parallel)
cl <- makePSOCKcluster(detectCores() - 1) # use all but one core available

clusterEvalQ(cl, library(raster))
clusterExport(cl, "rst")

## loop over countries of interest (in this example, austria and spain)
xtr <- parLapply(cl, c("AUT", "ESP"), function(country) {

# get country data
cnt <- getData("GADM", country = country, level = 2, path = "data/gadm")

# extract mean value per level-2 administrative area
val <- extract(rst, cnt, fun = mean, na.rm = TRUE, sp = TRUE)
dat <- data.frame(val)

# write to file and return
write.csv(dat, file = paste0(country, "meansdf.csv"),
quote = FALSE, row.names = FALSE)

return(dat) # or 'val'
})

## deregister parallel backend
stopCluster(cl)
``````
• Thank you very much, I will check all of the points you made! I will probably need to increase the resolution, so it is likely I am going to use the iteration to download each tile. – Marina Feb 7 '17 at 7:31
• @Marina please refer to the above update for 30 arc-seconds WorldClim data. Hope that clarifies! – fdetsch Feb 7 '17 at 10:30