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)
- Select Dataset: The first argument specifies the dataset. ‘SRTM’ returns the SRTM 90 elevation data.
- Specify Lon: The second argument specifies the lon of the SRTM tile.
- 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}:
If you have any questions, feel free to comment below.
Cheers!
Martin
10 Comments
You can post comments in this post.
Hi Martin,
great blog.
I have a troubles getting the country-wide SRTM download, mosaicking and cropping to work. The download stops after the first itereation because the GeoTiff file can’t be found. I am working on Ubuntu 16.04 and I guess the problem lies in the fact that the getData() function searches for pattern = “.TIF” and my unzipped images end on lower case pattern = “.tif”. I tried to source the getData() function but unfortunately it didn’t work. Do you have any work around suggestions?
Many thanks,
Valentin
Valentin Louis 7 years ago
Hey Valentin! This is interesting,… I never tried the function on Linux.
Try setting your working directory manually before the download. I.e. `setwd(“your/custom/dir/for/srtm/”)`.
Then wrap the `getData()` call inside a try-clause: `try(getData(“…”))`. This should make sure you go through all the iterations.
After everything is downloaded, call `rasters <- list.files(full.names=T, pattern=".tif$")` and stack them using `stack <- stack(rasters)`. Good luck!
Martin 7 years ago
Thanks, Martin.
The try-clause worked, but the stack didn’t because of the different extent. lapply() on the rasters did the job.
Best wishes,
Valentin
Valentin Louis 7 years ago
Happy to hear that you made it work! Cheers Martin
Martin 7 years ago
Awesome tutorial! I’ve been wrestling with various topographic datasets recently, and this method seems super handy. Any tips for how I might incorporate multiple countries? Or better yet, a particular region of the world (e.g. Central America and the Caribbean)? I’m trying to make a topographic map with multiple countries that’s still cropped to the coastline.
-Hollis
Hollis Dahn 7 years ago
I keep getting this error. Please guide me:
Muhammad Waqas 6 years ago
Hey Muhammad,
it seems that the path to the shapefile is wrong. Follow these steps:
1. Make sure you downloaded and unzipped the shapefile: http://www.gis-blog.com/wp-content/uploads/2017/01/srtm.zip
2. Find the location of your shapefile on your computer. For example: “C:/Downloads/srtm/tiles.shp” and copy it into your script instead of the path provided here: country_name shp <- shapefile("srtm/tiles.shp") Hope this works! C
Martin 6 years ago
I have been attempting to carry out this workflow using sf objects, I am struggling to understand how to extract the lon/lat from sf objects to work within the function.
Would you advise the operation of this in sf objects over sp?
My goal is to download SRTM DEM for 121 countries coastlines. I would also like to have this in 1 arc second but trying this at the moment.
Amy 4 years ago
Nice post, I have a question: how to merge 2 countries when the second is within the first bigger country? if so how
Abou 4 years ago
Hey Martin,
thank you for your code and the all the handy information.
I tried your code for Norway and everything seemed to work, but the plottet map doesn’t seem to exceed 60°N so I just got the most southern tip of Norway. Is there any way to get Data above 60°N using the same method in R?
Thank you very much!
Josephine
Josephine 4 years ago
Post A Reply