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

About This Author

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.

30 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 3 years ago Reply


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

      Martin 3 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 3 years ago Reply


  • What’s the best way to clip the resulting raster to a shapefile? i.e. if I only want Austria to show?

    Jenny 2 years ago Reply


    • Use the crop function:

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

      Martin 2 years 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 2 years 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 2 years 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 2 years ago Reply


    • I am using QGIS by the way.

      Bela Arora 2 years 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 1 year ago Reply


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

    Dev 1 year 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 1 year ago Reply


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

    rommel rojas 1 year 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 1 year 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 1 year 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 1 year 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 12 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 12 months ago Reply


  • Hi Martin, this basically downloads raster data for the current year. Is there any way to download data for a particular or a set of years?
    Thanks

    Shaurya 5 months ago Reply


    • The data represent the average of the years 1970-2000. It’s not possible to download single years BUT there is a possibility to download future projections of 2050 and 2070:

      “To get (projected) future climate data (CMIP5), you must provide arguments var and res as above. Only resolutions 2.5, 5, and 10 are currently available. In addition, you need to provide model, rcp and year. For example:”

      getData('CMIP5', var='tmin', res=10, rcp=85, model='AC', year=70)

      Martin 5 months ago Reply


  • Thank you very much. This was extremely helpful with regards to obtaining elevational data. Do you have any code for an overlay function? For example I would like to overlay a species richness map onto a tile of srtm in order to determine which species occur at which altitudinal band (200m) altitudinal band

    Nikhail Arumoogum 3 months ago Reply


  • Its a great site. I want to learn 2 things.
    1- supervised and unsupervised classification of Satellite data
    2- Regression Kriging

    Muhammad Mohsin Waqas 2 months ago Reply


    • Dear Muhammad!
      Thanks for your feedback. I will try to write up a short post on unsupervised classifcation with R in the upcoming days.
      Cheers Martin

      Martin 2 months ago Reply


      • Any idea about the regression kriging

        Muhammad Waqas 2 months ago Reply


  • Hello Martin,

    I want to use the function ‘GADM’ for Madagascar but it tells me
    “trying URL ‘http://biogeo.ucdavis.edu/data/gadm2.8/rds/MDG_adm0.rds’
    Content type ‘ êú ‘ length 734446 bytes (717 KB)
    downloaded 717 KB

    Could not download file — perhaps it does not exist
    NULL”

    Thanks

    Lucas Chapuis 1 month ago Reply


    • Hey Lucas! I checked the link and the download works for me. Can you try again? Maybe it was only temporarily disabled!

      Martin 1 month ago Reply


Post A Reply

*