Today I was looking for an high resolution Digital Elevation Model for a Region of Belgium. Until now, I was using the raster::getData("alt", country = "BE", mask = TRUE) which only provides a low resolution DEM (5 km).

The belgian National Geographic institute provides DEM’s at high resolution (80m and 20m). These are available from their download portal. Unfortunately I have not found yet a trick to download their data automatically (RSelenium?) as they ask you to click on some buttons to accept their terms of services.

After some researches on the web, I’ve found this helpful tutorial from gis-blog.com where the author proposed a little script to automate the download process of SRTM tiles based on the intersections of the extent of your region of interest and the extent of the SRTM tiles. It also relies on the raster::getData() function but rather than passing "alt" as dataset name, it uses "SRTM".

While trying the script, the getData always threw me the exact same error as reported on this stackoverflow thread. After entering in debug mode (debug(getData)), I’ve found that the first URL from which the function tries to download the SRTM was broken while the third one pointing to http://srtm.csi.cgiar.org/ was valid. So I decided to modify the function according to my needs.

I first forked the CRAN/raster mirror repo and edited the .SRTM function so that it first tries to download SRTM data from http://srtm.csi.cgiar.org/. I also had to change the file extension from “.TIF” to “.tif” as the files were created with a lowercase extension. You can find the full log of my edit on commit log of my raster fork repo. The last step was to make this adapted raster package available to my R installation instead of the CRAN version. For this, I simply used

devtools::install_github("pokyah/raster")
library(raster)

As I was finally able to automatically download the SRTM tiles from http://srtm.csi.cgiar.org/, the last remaining part of the task was to put all my code in a function. I wanted it to work not only for my current region of interest, so I decided to generalize it. my custom function allow you to construct a terrain raster stack with the desired resolution and CRS for your region of interest. To make it work, you first need to download the srtm tils.shp as mentioned in the gis-blog plost. Here is its source code with an example :

# @param country_code.chr a character specifying the ISO contrycode. Ex : BE for belgium
# @param NAME_1.chr a character specifying the NAME_1 value for lower than country level information
# @param aggregation_factor.num a numeric specifying the aggregation factor to get the desired spatial resolution
# @param EPSG.chr a character specifying the EPSG code of the desired Coordiante Reference System (CRS)
# @param path.chr a character specifying the path where to dowload the SRTM data
build_highRes_terrain_rasters.fun <- function(country_code.chr, NAME_1.chr=NULL, aggregation_factor.num=NULL, EPSG.chr=NULL, path.chr) {
  # Path to downloaded SRTM Tiles refs
  srtm.tiles.ref <- raster::shapefile("<PATH_TO_DOWNLOADED_TILES_REFS")

  # Get country geometry first
  extent.sp <- raster::getData('GADM', country=country_code.chr, level=1)

  if(!is.null(NAME_1.chr)){
    extent.sp <- subset(extent.sp, NAME_1 == NAME_1.chr)
  }

  # Intersect extent geometry and tile grid
  intersects <- rgeos::gIntersects(extent.sp, srtm.tiles.ref, byid=T)
  tiles      <- srtm.tiles.ref[intersects[,1],]

  # Download tiles using getData
  # inspired from https://www.gis-blog.com/download-srtm-for-an-entire-country/
  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', #data are downloaded from http://www.cgiar-csi.org/. See getData do of pokyah/raster repo on github
                    lon=lon,
                    lat=lat,
                    download = FALSE,
                    path = path.chr)

    srtm_list[[i]] <- tile
  }

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

  # Crop tiles to extent borders
  extent.elevation.ras <- raster::crop(srtm_mosaic.ras, extent.sp)
  extent.elevation.ras <- raster::mask(extent.elevation.ras, extent.sp)

  # transform to desired CRS
  if(!is.null(EPSG.chr)){
    raster::projectRaster(extent.elevation.ras, crs = toString((dplyr::filter(rgdal::make_EPSG(), code==EPSG.chr))$prj4))
  }

  # aggregate to lower resolution
  # inspired from https://stackoverflow.com/questions/32278825/how-to-change-the-resolution-of-a-raster-layer-in-r
  if(!is.null(aggregation_factor.num)){
    extent.elevation.ras <- raster::aggregate(extent.elevation.ras, fact=aggregation_factor.num)
  }

  # compute the slope from the elevation
  # inspired from https://rpubs.com/etiennebr/visualraster
  extent.slope.ras <- raster::terrain(extent.elevation.ras, opt="slope", unit="degrees")
  extent.aspect.ras <- raster::terrain(extent.elevation.ras, opt="aspect", unit="degrees")
  extent.roughness.ras <- raster::terrain(extent.elevation.ras, opt="roughness")

  # compute the aspect from the elevation
  extent.terrain.ras = raster::stack(
    extent.elevation.ras,
    extent.slope.ras,
    extent.aspect.ras,
    extent.roughness.ras)
}

Back to posts