How to download, reproject and export Digital Elevation Model using R ?

Note: Text inside shaded boxes are R commands, while all other texts are explanations. You need only run text inside shaded boxes in your R console to follow this exercise.

This script will show you how to download an SRTM Digital Elevation Model (DEM) of any country. This is accomplished using the raster package in R software.I assume that you have already install R and/or RStudio. We’ll begin the process by installing the raster package including its dependencies, which is necessary for downloading, reprojecting, and exporting the DEM. Check whether the raster package and its dependencies havebeen installed already to save installation time and if not, install them.

packages <- c('raster','rgdal','rgeos')

Warning messages:
1: multiple methods tables found for ‘crop’ 
2: multiple methods tables found for ‘extend’ 

install.packages(setdiff(packages, rownames(installed.packages())))

Now, we load the raster package and use the function getData to download the DEM.

library(raster)

package 㤼㸱raster㤼㸲 was built under R version 4.0.5Loading required package: sp

The getData function requires three major arguments (name, country, and mask). A character representing the name of the data set is provided with the name argument. The country argument uses ISO codes with three letters. You can obtain ISO code by typing getData(‘ISO3’) in the console. Mask is a logical argument for keeping or removing neighboring countries.

getData('ISO3')

getData('alt', country='NPL', mask=TRUE)

class      : RasterLayer 
dimensions : 516, 996, 513936  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 80, 88.3, 26.2, 30.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : NPL_msk_alt.grd 
names      : NPL_msk_alt 
values     : 60, 8682  (min, max)

The SRTM DEM for Nepal has been successfully downloaded using the above code. The data has been saved in your default working directory. Use the following line to determine where it is located.

getwd()

[1] "C:/Users/amthapa/Desktop"

By using the raster function from the raster package, it is possible to import output data into R in its native raster format (.grd).

SRTM_nepal=raster('NPL_msk_alt.grd')

We’ll make a simple default plot to see how the data looks. However, you can make more sophisticated plots using packages such as ggplot2, rasterVis, sp, etc., which is beyond the scope of this exercise.

plot(SRTM_nepal)

Now, we will be reprojecting our raster from wgs84 to utm and exporting it as GeoTiff. In fact, you can transform the raster into any projection system by using the epsg code of that projection system. Let’s do this. To do this, we use the projectRaster function from the raster package. R offers help for every function by using the ? followed by the function name. In order to get help for the projectRaster* function, type ?projectRaster in the console. There are several mandatory arguments with the projectRaster function, including from (raster to be reprojected), res (output resolution), crs (coordinate reference system), and format (output raster format). You can also specify a filename if you want to save the projected raster.

?projectRaster
nepal_srtm_utm45n=projectRaster(from=SRTM_nepal, res=90, crs='+init=epsg:32645', method="bilinear", 
              filename="nepal_srtm_dem_90m_utm45n.tif",format='GTiff',overwrite=TRUE)
plot(nepal_srtm_utm45n)

The reprojected DEM can be found in your working directory as GeoTiff files. This data can be explored in any GIS software you are familiar with.