Instantly share code, notes, and snippets.

Embed
What would you like to do?
Uncertainty in global climate databases: comparing Worldclim with Iberian climate atlas

Precipitation: Worldclim vs Iberian Climatic Atlas

Francisco Rodriguez-Sanchez 2017-05-16

library(knitr)
opts_chunk$set(cache = TRUE, message = FALSE)
library(raster)
library(rgis)
sp.extent <- c(-10, 5, 35, 44)

Worldclim v1.4

tile1 <- getData('worldclim', download = TRUE, var = 'prec', res = 0.5, 
                     lon = -5, lat = 36)
tile2 <- getData('worldclim', download = TRUE, var = 'prec', res = 0.5, 
               lon = 0, lat = 36)
wc1 <- merge(tile1, tile2)

worldclim1 <- sum(wc1)
worldclim1 <- crop(worldclim1, sp.extent)
#plot(worldclim1)

Worldclim v2

pmonth <- stack(unzip("C:/Users/FRS/Dropbox/GIS.layers/WORLDCLIM/v2/Present/wc2.0_30s_prec.zip", exdir = tempdir())[-13])
pmonth <- crop(pmonth, sp.extent)
worldclim2 <- sum(pmonth)
#plot(worldclim2)

Iberian Climate Atlas

atlas <- raster(ungzip("C:/Users/FRS/Dropbox/GIS.layers/AtlasClimaticoDigitalPenIb/pluvio_anual.nc.gz", ext = "gz"))
atlas.geo <- projectRaster(atlas, worldclim1)
atlas.geo <- atlas.geo/10
#plot(atlas.geo)

Compare precipitation in selected localities

locs <- cbind(x = c(-5.3541, -2.9125, -5.0875, -5.8625, -8.1458, -1.7875, -8.6541), 
              y = c(36.6958, 37.8375, 40.1958, 40.07916, 40.04583, 43.2208, 42.6791))
P.worldclimv1.4 <- worldclim1[cellFromXY(worldclim1, locs)]
P.worldclimv2 <- worldclim2[cellFromXY(worldclim2, locs)]
P.atlas <- atlas.geo[cellFromXY(atlas.geo, locs)]
knitr::kable(data.frame(locs, P.worldclimv1.4, P.worldclimv2, P.atlas), digits = c(4,4,0,0,0))
x y P.worldclimv1.4 P.worldclimv2 P.atlas
-5.3541 36.6958 824 796 1237
-2.9125 37.8375 539 509 1046
-5.0875 40.1958 366 377 1259
-5.8625 40.0792 599 814 1216
-8.1458 40.0458 1203 1228 1877
-1.7875 43.2208 1346 1330 2512
-8.6541 42.6791 1240 1480 1966

Compare maps

compar <- stack(worldclim1, worldclim2, atlas.geo)
names(compar) <- c("Worldclimv1.4", "Worldclimv2", "IberianAtlas")

locs.sp <- as.data.frame(locs)
coordinates(locs.sp) <- c("x", "y")
proj4string(locs.sp) <- projection(worldclim2)

library(rasterVis)
levelplot(compar, par.settings = viridisTheme, main = "Annual Precipitation") +
  layer(sp.points(locs.sp, pch = 1, col = "red", cex = 2))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment