Happy new year 2017!

It’s January 1st and I made the new years resolution to post more this year. 🙂 So here is a short tutorial on how to download SRTM (Shuttle Radar Topography Mission) 90m resolution data (3 arc seconds) for an entire country using the {raster} and {rgeos} R packages.

Basics (Download single SRTM tile)

As described in a previous post it’s possible to download the 90m SRTM using the getData() function as follows:

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 a sginle SRTM Tile somewhere around Vienna:

Advanced Download (entire country)

Specifying the longitude and latitude for the data download is often not very handy. Especially when we want to download SRTM data for an entire country. Usually a country consists of more than one tile and looking up the coordinates of each tile is a pain.

Here is a short script where you only have to specify the ISO country code and the tiles are automatically being mosaiced and plotted. You can find the code and all the data needed on github.

Preparation: Before we can start, download the shapefile which contains the the SRTM tile grid. It can be found here or on github. Without this file, the code below won’t work.

Once you have downloaded the shapefile, we can start the script. First load the necessary libraries which are:

library(raster)
library(rgeos)
library(rasterVis)

Then specify the country for your SRTM data download:

#Specify target ISO country code and path to downloaded shapefile
country_name <- "AUT"                        #Austria
shp          <- shapefile("srtm/tiles.shp")  #Path to downloaded SRTM Tiles

After these initial steps, you can run the code from below and lean back. The script automatically determines all the tiles that are needed for the download, mosaics and crops them into a single file:

#Get country geometry first
country <- getData("GADM", 
                   country = country_name, 
                   level=0)

#Intersect country geometry and 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('SRTM', 
                  lon=lon, 
                  lat=lat)
  
  srtm_list[[i]] <- tile
}

#Mosaic tiles
srtm_list$fun <- mean 
srtm_mosaic   <- do.call(mosaic, srtm_list)

#Crop tiles to country borders
srtm_crop     <- mask(srtm_mosaic, country)

#Plot
p <- levelplot(srtm_mosaic)
p + layer(sp.lines(country, 
                   lwd=0.8, 
                   col='darkgray'))

Here the result for Austria plotted with {rasterVis}:

Austria SRTM

If you have any questions, feel free to comment below.

Cheers!

Martin