如何分区域计算DEM高程均值与方差

“ DEM是提取地貌特征的重要手段，高程的均值和方差能够反映一个区域的高程和地貌起伏情况，是重要的地貌特征参数，本文以批量计算四川各县DEM高程均值与方差为例介绍计算方法。”
首先，我们要下载DEM，目前全球30m的DEM已经更新到了ASTER DEMv3，下载方法可以参考ASTER Global Digital Elevation Model全球DEM下载方法然后，我们对下载下来的DEM数据进行镶嵌和裁剪，根据需要进行分割。
最后，在R语言中计算DEM的高程均值和方差等参数，并生成表格。R语言相关基础知识请参考：R语言安装部署基础
R语言常用代码集

具体操作过程如下：

01

ArcGIS分割DEM

—

DEM的下载和拼接方法在ASTER Global Digital Elevation Model全球DEM下载方法文中已有详细说明，在这里不再赘述。

我们先来准备用于分割DEM的矢量，在ArcGIS中加载四川县域矢量，为了把DEM切成一个县城一块，我们需要把整个的县域矢量分为各县的单独矢量。

使用ArcToolbox-Analysis Tools-Extract-Split将整个的矢量分成独立的小矢量。

分拆后的矢量如下：

这样，就可以用这些分拆后的县分块矢量来分割DEM了。为了方便分割DEM，我们使用Model Builder建立一个处理模型。对矢量进行循环，分割整个的DEM，注意，分割结果要改为%Name%.tif，这样才能将分割结果存储为TIFF格式，且在文件夹中根据县名命名。

虽然在ArcMap的图层属性中，可以查看DEM的最大值、最小值、方差、均值，但是由于只能一个个查看，批量导出比较困难，在这里我选择用R语言进行批量计算。

02

R批量读取DEM

—

ArcGIS批量裁剪后的DEM如下所示：

每一块DEM都有多个包含“.tif”的文件，这个时候我们就需要使用下面的语句进行读取目录：

#读取文件名tifnames <- list.files(path = "D:/TEMP/chuandianLandscape/DEMparts", pattern = '.tif$')#为了读取栅格数据，我们需要将文件名和路径接起来：tifdir <- paste("D:/TEMP/chuandianLandscape/DEMparts/",tifnames,sep="")

我们使用list.names来提取TIFF文件名。patern中“.tif$”表示输出以“.tif”结尾的文件名，如果是“.tif”就会把所有的带“.tif”的文件全部读取。

03

R批量计算栅格参数

—

完整代码如下：

library(raster)tifnames <- list.files(path = "D:/TEMP/chuandianLandscape/DEMparts", pattern = '.tif$')tifdir <- paste("D:/TEMP/chuandianLandscape/DEMparts/",tifnames,sep="") #批量读取DEM文件allDEMS <- lapply(tifdir, raster)#计算DEM最大值Elemax <- sapply(allDEMS, maxValue)#计算DEM最小值Elemin <- sapply(allDEMS, minValue)#计算均值Elemean <- sapply(allDEMS, cellStats, stat="mean")#计算标准差Elesd<- sapply(allDEMS, cellStats, stat="sd")#将计算结果和文件名合并ElevationSummary <- cbind(tifnames, Elemin, Elemax, Elemean, Elesd)#输出计算结果write.csv(ElevationSummary, file = "DEM参数计算.csv")

有关文章与参考文献：

https://www.earthdatascience.org/courses/earth-analytics/lidar-raster-data-r/introduction-to-spatial-metadata-r/

https://stackoverflow.com/questions/52746936/how-to-efficiently-import-multiple-raster-tif-files-into-r