R {raster}: Data acquisition – SRTM, Worldclim, Global adm. boundaries

Hey!

Today I will show how powerful the R {raster} package is on another example. The raster package is not only a great tool for raster processing and calculation but also very useful for data acquisition. With the function getData() you can download the following data directly into R and process it:

  • SRTM 90 (elevation data with 90m resolution between latitude  -60 and 60)
  • World Climate Data (Tmin, Tmax, Precip, BioClim)
  • Global adm. boundaries (different levels)

Global adm. boundaries

Lets start with the global adm. boundaries:

Install the raster package and load it first. To be able to use the getData() function to acquire data about global amd. boundaries we have to specify three arguments:

install.packages("raster")
library(raster)
austria0 <- getData('GADM', country='AUT', level=0)
  1. Select Dataset: The first argument specifies the dataset. ‘GADM’ returns the global administrative boundaries.
  2. Select country: The second argument provides the country name of the boundaries by using its ISO A3 country code (more info here)
  3. Specify level: The third argument specifies the level of of administrative subdivision (0=country, 1=first level subdivision).

The code above returns the boundaries for Austria for the level 0. Lets compare them to the Level 1 subdivision by plotting both of them:

#Get Data
austria0 <- getData('GADM' , country="AUT", level=0)
austria1 <- getData('GADM' , country="AUT", level=1)

#Plot
par(mfrow(2,1))
plot(austria0, main="Adm. Boundaries Austria Level 0")
plot(austria1, main="Adm. Boundaries Austria Level 1")

austria boundaries

World Climate

Lets do the same with the World Climate data, here you also have to specify three arguments:

climate <- getData('worldclim', var='bio', res=2.5)
  1. Select Dataset: The first argument specifies the dataset. ‘worldclim’ returns the World Climate Data.
  2. Select variable: The second argument specifies the variable: ‘tmin’, ‘tmax’, ‘prec’ and ‘bio’ (more info here).
  3. Specify resolution:  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.

The code above returns a raster with the 18 bioclimate variables covering the whole world with a resoltion of 2.5 minutes of degrees:

BIO1 = Annual Mean Temperature
BIO2 = Mean Diurnal Range (Mean of monthly (max temp – min temp))
BIO3 = Isothermality (BIO2/BIO7) (* 100)
BIO4 = Temperature Seasonality (standard deviation *100)
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO10 = Mean Temperature of Warmest Quarter
BIO11 = Mean Temperature of Coldest Quarter
BIO12 = Annual Precipitation
BIO13 = Precipitation of Wettest Month
BIO14 = Precipitation of Driest Month
BIO15 = Precipitation Seasonality (Coefficient of Variation)
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter

Lets plot the first indicator “Annual Mean Temperature”:

#Plot
plot(climate$bio1, main="Annual Mean Temperature")

Mean annual temperature

SRTM 90 Elevation

Last but not least, lets have a look at the SRTM 90 Data. We will use the getData() function one last time:

srtm <- getData('SRTM', lon=16, lat=48)
  1. Select Dataset: The first argument specifies the dataset. ‘SRTM’ returns the SRTM 90 elevation data.
  2. Specify Lon: The second argument specifies the lon of the SRTM tile.
  3. Specify Lat:  The second argument specifies the lat of the SRTM tile.

The code above will return one SRTM Tile somewhere around Vienna. Let’s plot the adm. boundaries of Austria together with the SRTM Tile in one plot:

plot(srtm)
plot(austria0, add=TRUE)

SRTM Austria
As you can see, not all of Austria is covered by this tile. Lets download two more eastern tiles and mosaic them to get the full extent of Austria. You can get more info on the tile extent here.

#Download two more tiles
srtm2 <- getData('SRTM', lon=13, lat=48)
srtm3 <- getData('SRTM', lon=9, lat=48)

#Mosaic/merge srtm tiles
srtmmosaic <- mosaic(srtm, srtm2, srtm3, fun=mean)

Lets plot the result:

plot(srtmmosaic, main="Elevation (SRTM)")
plot(austria0, add=TRUE)

SRTM Austria Total

I hope this gave you a good first impression on how to acquire spatial data with R.  More useful tutorials will follow, until then you can get more information here:

http://www.inside-r.org/packages/cran/raster/docs/getData

Martin

Martin was born in Czech Republic and studied at the University of Natural Resources and Life Sciences, Vienna. He is currently working at GeoVille - an Earth Observation Company based in Austria, specialised in Land Monitoring. His main interests are: Open-source applications like R, (geospatial) statistics and data-management, web-mapping and visualization. He loves travelling, geocaching, photography and sports.

22 Comments

You can post comments in this post.


  • Thank you immensely for your effort. I was stuck in getting data in R for my analysis and stumbled on your work and the ‘mosaic’ function really did save my day.

    Nwankwo Emmanuel 2 years ago Reply


    • Hey! I am always happy to hear that I could help.
      Cheers Martin

      Martin 2 years ago Reply


  • Thanks a lot for sharing. Never will I get the geo-data in R without you. But there’re something I feel confused. When combining the global-adm. boundaries- data with SRTM data, are we certain abou they are using the same spatial reference? If possible, is there any solution in R to cliping the srtm data and make it match with the boundary? Thanks

    Xue Yin 2 years ago Reply


  • Best way to clip resulting raster to shapefile? i.e. if I only want Austria to show?

    jenny 1 year ago Reply


    • Use the crop function:

      image < - crop(srtmmosaic, austria0)
      plot(image)
      

      Martin 1 year ago Reply


  • Hi,

    I am mapping the alps using ‘getData(“SRTM”, lon=12, lat=48), but when i am plotting it i have huge blank part at each side of the map and I do not know how to remove it. do you know a way to do so?

    Thank you!
    Sylvie

    Sylvie Pighini 1 year ago Reply


    • If you are working in RStudio, simply try to adjust the width of the plot window (by dragging it with the mouse). The white stripes should then dissapear.

      Martin 1 year ago Reply


  • Hi Martin,

    I am trying to extract statistic values (min, max, median, mean, range) from BioClim rasters. When I divide the entire world raster by 10 [for temperature variable (Bio 1)], the values calculated, eg for mean T for my polygon, changes for some of the polygons.

    Another thing: When I try to clip raster to polygons, it gives error message “Cutline feature without a geometry”. The CRS is same for my shapefiles and raster.

    What am I doing wrong?

    Thanks!
    B.

    Bela Arora 1 year ago Reply


    • I am using QGIS by the way.

      Bela Arora 1 year ago Reply


  • Is there a way to get WorldClim data for a specific country or do you always need to download the layer for the entire world?

    Sarah 1 year ago Reply


    • Unfortuately you always have to download the layer for the entire world. But you can also download the Adm Boundaries of your target country and then use this shapefile to crop the WorldClim layer to the desired extent.

      Martin 11 months ago Reply


  • Hi Martin
    How to use get data for downloading SRTM or alt data for whole world?
    Liebe Gruße
    Dev

    Dev 11 months ago Reply


    • You can use the following link: http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp and then click on “Enable Mouse Drag” and draw a rectangle over the entore globe. Doing this manually is probably faster than programmatically. Cheers and happy new year!

      Martin 11 months ago Reply


  • hello martin
    how i can get paleo climate (past variables) with R…?

    rommel rojas 10 months ago Reply


  • Hi Martin,

    I am trying to download elevation data of Belgium using code:
    bel <- getData('alt', country='BEL', mask=TRUE)

    However it always shows http status 403 forbidden. Any suggestions?

    Xin 10 months ago Reply


  • Hi Martin? Do you know how could I make a rasterstack right after using the codes above and getting the worldclim data?

    MG 10 months ago Reply


    • Hey!
      The output of the worldclim getData call is a rasterStack already. If you want to stack the SRTM and wordclim in one stack, you need to harmonize the resolution and extent of both datasets.
      Have a look into these functions `resample`, `crop` and `stack`. Note: Also the projection of datasets have to be the same.

      Martin 10 months ago Reply


  • Hi Martin,
    I have 32 years precipitation data in tiff format. I have 32 tiff raster files. Every tiff/raster is a single year. Every year/tiff is comprised of 365/6 days/bands. I would like to calculate mean annual precipitation. How can I first sum precipitation of every day in to annual sum of precipitation? Then I want to calculate the mean annual of precipitation. I want to use R.
    Best regards,
    Kibru

    Kibru 7 months ago Reply


    • Dear Kibru!
      What exactely do you mean by 365/6 days/bands?
      You need to use loops here, in combination with the calc function from the raster R package:

      Here a short overview on how to use the calc function here:

      library(raster)
      #yearly sum

      yearly_stack <-stack(list.files("pathToYearOne", full.names=T, pattern="tif$")) yearly_sum <- calc(yearly_stack, sum, na.rm=T) yearly_mean <- calc(yearly_stack, mean na.rm=T) If you can tell me more about your data/folder structure, I might be able to help you with the set up of the loop. If processing takes too long, you should look into the clusterR function. Cheers Martin

      Martin 6 months ago Reply


Post A Reply

*