read NASADEM .slope file into R and change format to .tif file

Question

I downloaded some tiles using NASA EARTHDATA for the NASADEM Slope and Curvation Global 1 arc second V001 data set. When I extracted the files, I saw the filenames followed the pattern: "TileLocation.Variable". For example, a tile with slope data located on the eastern Mediteranian is called: "n36e035.slope". I was surprised the file extension was ".slope" and not ".tif".

When I tried reading the file into R with r <- raster::raster("Filepath/n36e035.slope") I get an error, because the file structure is not a tif or geotiff. I want to read multiple tiles for each variable (slope, curvature, etc.), merge, crop to my study area, then write out the combined raster to my local device as a .tif file. That way the DEM file format and data structure will match the other rasters I have.

My preferred language is R, but if there's another open-source way to change the file format from this NASA-specific extension to a .tif then I can use that. I tried looking online, but all the tutorials used Google Earth Engine or Arc and didn't allow me to save the combined .tif files locally.

You can download the n36e35 zip here and the n35e35 zip here. You may need to log-in to NASA EARTHDATA to view and download the DEM tiles. The overview is here and the user guide is available here, but the user guide is more about how the data set was made, not how to read it in or change the data format. One strange note is that the DEM this data set is based off of has an .hgt file extension, which I can easily read into R with the raster::raster function.

Answer 1

Regrettably, NASA does not provide header files so you need to create them yourself.

To help with that, I added a function makeVRT to the terra version 1.5-9. That is currently the development version, you can install it with install.packages('terra', repos='https://rspatial.r-universe.dev'). I use that function in demhdr below (after looking at this github repo by David Shean) that sets the specific parameters for these files. Parameters are also give here but note that on that page the datatype for SWB is incorrectly given as 8-bit signed integer, whereas it should be 8-bit unsigned integer.

demhdr <- function(filename, ...) {
    f <- basename(filename)
    stopifnot(tools::file_ext(f) != "vrt")
    sign <- ifelse(substr(f, 1, 1) == "n", 1, -1)
    lat <- sign * as.numeric(substr(f, 2, 3))
    sign <- ifelse(substr(f, 4, 4) == "e", 1, -1)
    lon <- sign * as.numeric(substr(f, 5, 7))
    if (grepl("aspect", f) || grepl("slope", f)) {
        datatype <- "INT2U"
    } else if (grepl("swb", f)) {
        datatype <- "INT1U"
    } else {
        datatype <- "FLT4S"
    }
    name <- unlist(strsplit(f, "\\."))[2]
    terra::makeVRT(filename, 3601, 3601, 1, xmin=lon, ymin=lat, xres=1/3600,
           lyrnms=name, datatype=datatype, byteorder="MSB", ...) 
}

For a folder with these files that do not end on ".vrt":

ff <- grep(list.files("."), pattern="\\.vrt$", invert=TRUE, value=TRUE)
ff
#[1] "n37e037.aspect" "n37e037.planc"  "n37e037.profc"  "n37e037.slope" 
#[5] "n37e037.swb"   

You can use the demhdr function like this:

fvrt <- sapply(ff, demhdr, USE.NAMES=FALSE)
fvrt
#"n37e037.aspect.vrt"  "n37e037.planc.vrt"  "n37e037.profc.vrt" 
# "n37e037.slope.vrt"    "n37e037.swb.vrt" 

And then, with files for a single tile, you can do

library(terra)
r <- rast(fvrt)
r
#class       : SpatRaster 
#dimensions  : 3601, 3601, 5  (nrow, ncol, nlyr)
#resolution  : 0.0002777778, 0.0002777778  (x, y)
#extent      : 36.99986, 38.00014, 36.99986, 38.00014  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : n37e037.aspect.vrt  
#              n37e037.planc.vrt  
#              n37e037.profc.vrt  
#              ... and 2 more source(s)
#names       : aspect, planc, profc, slope, swb 

Note the very unfortunate georeferencing of NASA SRTM data. The data would have lined up with other lon/lat raster data, and would have been much more usable if the extent would have been 37.0, 38.0, 37.0, 38.0 and the number of rows and columns would have been 3600. But that is not the case.

plot(r)

This tile did not seem to have no data values; in other tiles you may need to set it with the NAflag argument in makeVRT or by using NAflag(x) <- on single layer SpatRaster

For aspect, it looks like you could use scale=0.01 to get values between 0 and 360 degrees)

To merge many tiles for say aspect, you should be able to do something like

fasp <- grep("aspect", fvrt, value=TRUE)
x <- src(lapply(fasp, rast))
m <- merge(x)

or make a new VRT file that combines the tiles, like this

vrt(fasp, "aspect.vrt")
m <- rast("aspect.vrt")

To read the files with raster

library(raster)
s <- stack(fvrt)