将 NASADEM .slope 文件读入 R 并将格式更改为 .tif 文件

Question

我使用 NASA EARTHDATA for the NASADEM Slope and Curvation Global 1 arc second V001 数据集下载了一些图块。当我提取文件时，我看到文件名遵循以下模式：“TileLocation.Variable”。例如，具有位于东地中海的坡度数据的瓦片称为：“n36e035.slope”。我很惊讶文件扩展名是“.slope”而不是“.tif”。

当我尝试使用 r <- raster::raster("Filepath/n36e035.slope") 将文件读入 R 时出现错误，因为文件结构不是 tif 或 geotiff。我想为每个变量（斜率、曲率等）读取多个图块，合并、裁剪到我的研究区域，然后将组合栅格作为 .tif 文件写入我的本地设备。这样 DEM 文件格式和数据结构将匹配我拥有的其他栅格。

我的首选语言是 R，但如果有另一种开源方法可以将文件格式从这个 NASA 特定的扩展名更改为 .tif，那么我可以使用它。我尝试在线查找，但所有教程都使用 Google Earth Engine 或 Arc，并且不允许我在本地保存组合的 .tif 文件。

您可以下载 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，但用户指南更多的是关于如何制作数据集，而不是如何读取或更改数据格式。一个奇怪的注意事项是，此数据集所基于的 DEM 具有 .hgt 文件扩展名，我可以使用 raster::raster 函数轻松将其读入 R。

Answer 1

遗憾的是，NASA 不提供头文件，因此您需要自己创建它们。

为了解决这个问题，我在 terra version 1.5-9 中添加了一个函数 makeVRT。那是目前的开发版本，你可以用install.packages('terra', repos='https://rspatial.r-universe.dev')安装它。我在下面的 demhdr 中使用了该函数（在查看了这个 github repo by David Shean) that sets the specific parameters for these files. Parameters are also give here 但请注意，在该页面上，SWB 的数据类型被错误地指定为 8 位有符号整数，而它应该是 8 位无符号整数.

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", ...) 
}

对于包含这些不以“.vrt”结尾的文件的文件夹：

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

您可以像这样使用 demhdr 函数：

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

然后，对于单个图块的文件，您可以执行

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

请注意 NASA SRTM 数据的非常不幸的地理配准。数据将与其他 lon/lat 栅格数据对齐，如果范围是 37.0, 38.0, 37.0, 38.0 并且行数和列 应该是 3600。但事实并非如此。

plot(r)

这个图块似乎没有数据值；在其他图块中，您可能需要使用 makeVRT 中的 NAflag 参数或在单层 SpatRaster

中使用 NAflag(x) <- 来设置它

对于方面，看起来你可以使用 scale=0.01 来获得 0 到 360 度之间的值）

要为某个方面合并多个图块，您应该能够执行类似

的操作

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

或者制作一个新的 VRT 文件，像这样组合这些图块

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

阅读文件 raster

library(raster)
s <- stack(fvrt)

将 NASADEM .slope 文件读入 R 并将格式更改为 .tif 文件

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

tiff

r

raster

gdal

rgdal