1.背景介紹
我網盤里有12.5m的全球分辨率數據，但是沒有分縣級單位裁剪的DEM數據。每一次使用都非常麻煩。
2.簡介
數字高程模型（Digital Elevation Model)，簡稱DEM，是通過有限的地形高程數據實現對地面地形的數字化模擬。
目前能獲取到的最高的DEM數據是12.5m分辨率的ALOS衛星數據。原始數據下載網址為:https://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個小時內回復。