1.背景介绍

我网盘里有12.5m的全球分辨率数据，但是没有分县级单位裁剪的DEM数据。每一次使用都非常麻烦。

2.简介

数字高程模型（Digital Elevation Model)，简称DEM，是通过有限的地形高程数据实现对地面地形的数字化模拟。 目前能获取到的最高的DEM数据是12.5m分辨率的ALOS卫星数据。原始数据下载网址为:urs.earthdata.nasa.gov.

3.数据制作

3.1 流程图

总体思路是：全国12.5m的DEM数据筛选、按县级单位进行裁剪、镶嵌，结果数据后处理。

3.1数据准备工作

首先是下载数据，筛选出中国区的范围。

3.2 DEM按省份、地级市、县裁剪

3.2.1 脚本进行裁剪

使用python的RasterIO模块进行单景的DEM读取，使用Geooandas模块操作省份矢量裁剪。这里暂时只放重要的裁剪脚本：

#裁剪函数
def clip(pathDir,shpdata,rasterfile):
    for i in tqdm(range(len(pathDir))):
        # 读入栅格文件
        rasterfile = files_path+"\\"+pathDir[i]
        rasterdata = rio.open(rasterfile)
        #获取栅格信息
        profile = rasterdata.profile
        #标识符
        note = pathDir[i]
        # 投影变换，使矢量数据与栅格数据投影参数一致
        shpdata = shpdata.to_crs(rasterdata.crs)
        # 按照所有矢量进行循环裁剪
        for j in range(0, len(shpdata)):
        try:
                # 获取矢量数据的features
                geo = shpdata.geometry[j]
                #获取该要素的属性信息
                data_shp_name=shpdata.全称[j]
                #文件保存位置的文件夹 各省
                data_filepath=str(data_shp_name)
                feature = [geo.__geo_interface__]
                # 通过feature裁剪栅格影像
                out_image, out_transform = rio.mask.mask(rasterdata, feature, all_touched=True, crop=True, nodata=0)
                profile.update(
                    height=out_image.shape[1],
                    width=out_image.shape[2],
                    shape=(out_image.shape[1],out_image.shape[2]),
                    nodata=0,
                    bounds=[],
                    transform=out_transform,
                    )
                # 定义要创建的目录
                mkpath = "目录名"
                # 检测目录是否存在
                mkdir(mkpath)
                # 文件名字
                name="文件名"
                with rasterio.open(name, mode='w', **profile) as dst:
                    dst.write(out_image)
        except:
                pass
                

裁剪后的影像单个县会出现多张影像，因此需要进行镶嵌。

3.2 DEM按县级单位进行镶嵌

遍历所有地区的文件夹，镶嵌使用gdal库的Warp函数。

#镶嵌一个文件夹中的所有文件
import os
from osgeo import gdal, gdalconst
import rasterio as rio
import rasterio.mask
import rasterio
from tqdm import tqdm

tifPath = '待镶嵌文件的目录'  # 待融合的图像所在的文件夹

tifPaths_folder = os.listdir(tifPath)

# 循环目录
for path in tqdm(tifPaths_folder):
    DEM_SMALL_PATH = os.path.join(tifPath, path)
    for (pathname, dirs, files) in os.walk(DEM_SMALL_PATH):
        son_Paths_file=files
    # 如果影像数量大于一
    if len(son_Paths_file) >= 2:
        DEM_SMALL_PATH2=DEM_SMALL_PATH+"\\"
        # 循环子目录，进行镶嵌
        #循环同一个文件下的tif文件
        inputFiles = []
        for path_small in son_Paths_file:
            # 每一个栅格的路径
            son_Paths_PATH = os.path.join(DEM_SMALL_PATH2, path_small)
            #读取影像
            inputrasfile = gdal.Open(son_Paths_PATH, gdal.GA_ReadOnly)  # 读取影像

            inputProj = inputrasfile.GetProjection()  # 获取坐标系


            inputFiles.append(inputrasfile)  # 推入列表

            options = gdal.WarpOptions(srcSRS=inputProj,  # 输入坐标系
                                       dstSRS=inputProj,  # 输出坐标系
                                       format='GTiff',  # 图像格式
                                       resampleAlg=gdalconst.GRIORA_NearestNeighbour,  # 重采样算法，这里是双线性内插
                                       # dstNodata=Nodata,  # 缺省值
                                       #    cutlineLayer=outline,  # 输出范围，这里可以是一个外轮廓shp数据
                                       outputType=gdalconst.GDT_Int16)  # 数据类型，这里是有符号32位整型

        #输出文件名
        outputfilePath =DEM_SMALL_PATH+"\\"+path+"_DEM"+"_12.5分辨率_ALOS数据"+".tif"
        #写栅格
        gdal.Warp(outputfilePath, inputFiles, options=options)  # 图像镶嵌

3.3 数据后处理

主要是使用python脚本，对每一个省份文件夹中的结果影像进行重命名，并删除多余的切片文件。这一步需要添加一个判定函数，判定是否为镶嵌文件， 是则保留，不是则删除。

# 删除非融合的原始文件  并将单个影像重新命名
import os
# 待融合的图像所在的文件夹
tifPath = '存储县级DEM的目录'  # 待融合的图像所在的文件夹

tifPaths_folder = os.listdir(tifPath)

# 循环目录
for path in tifPaths_folder:
    DEM_SMALL_PATH = os.path.join(tifPath, path)
    for (pathname, dirs, files) in os.walk(DEM_SMALL_PATH):
        son_Paths_file=files
        # 如果影像数量大于一
        if len(son_Paths_file) >= 2:
            DEM_SMALL_PATH2=DEM_SMALL_PATH+"\\"
            # 循环子目录，进行镶嵌
            #循环同一个文件下的tif文件
            inputFiles = []
            for path_small in son_Paths_file:
                #判定不是镶嵌影像，删除
                if path_small.replace("_DEM_30m分辨率_ASTER数据.tif","")!=path:
                    # 每一个栅格的路径
                    son_Paths_PATH = os.path.join(DEM_SMALL_PATH2, path_small)
                    os.remove(son_Paths_PATH)
        else:
            for (pathname, dirs, files) in os.walk(DEM_SMALL_PATH):
                son_Paths_file = files
            for path_small in son_Paths_file:
                DEM_SMALL_PATH2 = DEM_SMALL_PATH + "\\"
                son_Paths_PATH = os.path.join(DEM_SMALL_PATH2, path_small)
                new_name=DEM_SMALL_PATH2+path+"_DEM_12.5分辨率_ALOS数据.tif"
                os.rename(son_Paths_PATH,new_name)

4.结果展示

裁剪得到了全国2700多个县12.5m的DEM影像。以加载四川省-资阳市-乐至县的DEM数据为例。

首先选择对应的省份、市，然后找到县，下载数据，加载到arcgis中:

5.总结

1.使用gdal、geopandas可以很方面地使用矢量裁剪栅格。 2.此次镶嵌的范围只是基于县的，没有超出gdal默认的镶嵌范围，因此不需要设置输出影像的范围大小。如果输出范围过大，是需要设置输出范围的。 3.数据镶嵌后，需要设置一个函数，批量删除文件夹中的中间文件，参考节3.3； 4.针对全国的DEM数据分县的裁剪，算法不是问题，电脑的算力才是主要问题。

6.数据分享

需要哪个县的DEM数据，请在评论区留言，我会在12个小时内回复。