[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())




#Specify target ISO country code and path to downloaded shapefile

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 <- getData("GADM", 
                   country = country_name, 

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

#Download tiles
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, 

  srtm_list[[i]] <- tile


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

#Load raster mosaic into R 
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) 

write.csv(AFGmeansdf, file= "AUTmeansdf.csv")
  • 1
    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):


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

## download worldclim 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:

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))
  download.file(url, destfile)
  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
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
  if (!dir.exists("data/gadm")) dir.create("data/gadm")
  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
  • 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.