Downloading environmental data

Now we can download the bioclim data directly from Worldclim. The environmental data are made up of 19 bioclimatic variables, the details of which you can read about here: http://www.worldclim.org/bioclim

The getData() function can download the bioclim data at a variety of resolutions. Valid resolutions are 0.5, 2.5, 5, and 10 (minutes of a degree). 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.

If you are a species of fairly restricted range you may be best of downloading the 0.5 resolution data but if the species is more broadly distributed it is probably best to use the 2.5 resolution.

env <- getData("worldclim", var="bio", res=2.5)

Once downloaded you can look at the resulting data by just running the command ‘env’

This will show you the RasterStack of downloaded data - this is a stack of 19 layers of raster data, one for each of the bioclim variables. You can also see some other basic information such as the spatial resolution of each cell and the extent of the rasters. As a side note rasters have to have the same extent, resolution and projection in order to be in a stack together.

env

## class      : RasterStack 
## dimensions : 3600, 8640, 31104000, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names      :  bio1,  bio2,  bio3,  bio4,  bio5,  bio6,  bio7,  bio8,  bio9, bio10, bio11, bio12, bio13, bio14, bio15, ... 
## min values :  -272,    35,     8,    64,   -86,  -559,    53,  -278,  -501,  -127,  -506,     0,     0,     0,     0, ... 
## max values :   312,   201,    96, 22704,   489,   258,   725,   376,   365,   382,   289, 10577,  2437,   697,   265, ...

Attaching names to the data and plotting it

We can change the names of the layers in the RasterStack to make them a bit more meaningful. Here I have written out the names of the bioclim variables and then attached them to the RasterStack using the names() function.

We can also plot the environmental data, If we just plot them all we get a small plot of global data of the first 16 layers.

bioclim_names <-
  c(
    "Annual_Mean_Temp",
    "Mean_Diurnal_Range",
    "Isothermality",
    "Temp_Seasonality",
    "Max_Temp_Warmest Month",
    "Min_Temp_Coldest_Month",
    "Temp_Annual_Range",
    "Mean_Temp_Wettest_Quarter",
    "Mean_Temp_Driest_Quarter",
    "Mean_Temp_Warmest_Quarter",
    "Mean_Temp_Coldest_Quarter",
    "Annual_Precip",
    "Precip_Wettest_Month",
    "Precip_Driest_Month",
    "Precip_Seasonality",
    "Precip_Wettest_Quarter",
    "Precip_Driest_Quarter",
    "Precip_Warmest_Quarter",
    "Precip_Coldest_Quarter"
  )

names(env) <- bioclim_names

plot(env)

We can look in more detail at the rasters within the stack - here are the first and twelfth layers, annual mean temperature and annual precipitation, respectively.

plot(env[[1]], main = names(env[[1]]))

plot(env[[12]], main = names(env[12]))

You might notice that some of the values seem a bit odd - in particular the temperature values which are a lot larger than they should be. This is to do with how the numbers are stored - it takes up less space to save an integer than a decimal so in order to maintain the accuracy of the data they are stored as values which are ten times larger than the ‘real’ values. We can fix this later by dividing these temperature layers by ten.

Cropping the environmental data

The data is currently at a global level, it needs to be cropped to the area of interest so that Maxent does not try to run the model for the whole world as that will take a really long time.

We can read back in the location data from Lab 1 and then crop the raster data to the extent of those locations.

capg <- read.csv("~/SDM_Course/Species_locs/capra_locs.csv")

To do this we need to create an extent object, which is a set of coordinates in this order - minimum longitude, maximum longitude, minimum latitude and finally maximum latitude.

Here we have taken the maximum and minimum values from the location points then extenced them by 1 degree so that the rasters will be cropped to an area slightly larger than the extent of the location points.

buff <- 1   #a buffer of one degree around the raster

xmin <- min(capg$lon) - buff
xmax <- max(capg$lon) + buff
ymin <- min(capg$lat) - buff
ymax <- max(capg$lat) + buff

e <- extent(xmin, xmax, ymin, ymax)

envcrop <- crop(env, e)

If we look at the cropped RasterStack by running the command ‘envcrop’ we should see that all of the information is still there but the extent values have changed

envcrop

## class      : RasterBrick 
## dimensions : 162, 236, 38232, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : 5, 14.83333, 42.91667, 49.66667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : Annual_Mean_Temp, Mean_Diurnal_Range, Isothermality, Temp_Seasonality, Max_Temp_Warmest.Month, Min_Temp_Coldest_Month, Temp_Annual_Range, Mean_Temp_Wettest_Quarter, Mean_Temp_Driest_Quarter, Mean_Temp_Warmest_Quarter, Mean_Temp_Coldest_Quarter, Annual_Precip, Precip_Wettest_Month, Precip_Driest_Month, Precip_Seasonality, ... 
## min values :              -87,                 42,            22,             4733,                      4,                   -174,               178,                      -139,                     -107,                       -25,                      -143,           579,                   67,                   9,                  7, ... 
## max values :              159,                109,            38,             7855,                    308,                     62,               326,                       213,                      233,                       236,                        98,          2543,                  247,                 172,                 47, ...

Now we can plot this cropped area and add the occurrence points on top.

plot(envcrop[[1]], main = "Annual Mean Temperature")
map(
  'worldHires',
  xlim = c(xmin, xmax),
  ylim = c(ymin, ymax),
  fill = F,
  add = T
)

points(capg$lon, capg$lat, pch = "+")

Dividing the temperature data by 10

Now we can divide the temperature layers by 10 to get the ‘real’ values. The temperature layers are the ones which are numbered, 1,2,5,6,7,8,9,10 and 11. To do this we create a vector called ‘ten_div’, selecting the temperature layers, then write a ‘for loop’ which is essentially saying for each of ‘layer’ in ‘ten_div’ take that layer and divide it by ten and then put it back in the stack in the same location or layer it came from.

ten_div <-
  c(1, 2, 5, 6, 7, 8, 9, 10, 11)  #the layers we want to divide by ten

for (layer in ten_div) {
  envcrop[[layer]] <- envcrop[[layer]] / 10
}

Saving selected bioclim layers in ASCII format

Now you need to pick which of the bioclimatic variables you want to include in the Maxent model, try and pick ones which are most ecologically relevant to your species. Here I have picked the bioclimatic variables 1, 4, 11, 12, 14 and 19. They need to be exported as ASCII files as this is the file type required by Maxent.

layers <- c(1, 4, 11, 12, 14, 19)

writeRaster(
  stack(envcrop[[layers]]),
  paste0("bio", layers), #this names the files with a prefix of bio and then the number of the bioclim layer
  bylayer = TRUE,
  format = 'ascii',
  overwrite = T
)

Now you have a folder with your environmental ASCII layers in, which you can now use with your species locations to build a model in Maxent.