我网盘里有12.5m的全球分辨率数据，但是没有分县级单位裁剪的DEM数据。每一次使用都非常麻烦。
数字高程模型（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:
geo = shpdata.geometry[j]
data_shp_name=shpdata.全称[j]
data_filepath=str(data_shp_name)
feature = [geo.__geo_interface__]
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:
裁剪后的影像单个县会出现多张影像，因此需要进行镶嵌。
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+"\\"
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,
outputType=gdalconst.GDT_Int16)
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+"\\"
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中:
1.使用gdal、geopandas可以很方面地使用矢量裁剪栅格。
2.此次镶嵌的范围只是基于县的，没有超出gdal默认的镶嵌范围，因此不需要设置输出影像的范围大小。如果输出范围过大，是需要设置输出范围的。
3.数据镶嵌后，需要设置一个函数，批量删除文件夹中的中间文件，参考节3.3；
4.针对全国的DEM数据分县的裁剪，算法不是问题，电脑的算力才是主要问题。
6.数据分享
需要哪个县的DEM数据，请在评论区留言，我会在12个小时内回复。
- 165
-
AirtestProject
Python
