简介浙江省绍兴市30米分辨率DEM数字高程数据包面向GIS从业者、城市规划与地理科研人员可用于地形分析、坡度计算、洪水模拟、视域分析等场景。压缩包共12个文件总大小约25.19MB核心为tif格式的高程栅格配套保存投影信息的tfw与金字塔优化文件ovr同时包含绍兴市行政边界shp矢量文件及其dbf属性表、prj坐标系、sbn/sbx空间索引、xml元数据等附件构成了一套可直接在ArcGIS、QGIS等软件中读取分析的标准地理数据集合。数据已覆盖绍兴市全域地形30米精度足以清晰呈现山丘、河谷等主要地貌特征辅助识别地质灾害风险区域、优化交通线路与城市空间布局。已有785人学习适合需要开展区域地形建模、环境评价或生成三维地形底图的相关人员使用。1. 拿到绍兴市DEM数据包先别急着拖进ArcGIS一个名为“浙江省绍兴市DEM数字高程数据含区域范围shp文件.zip”的压缩包拆开后是DEM栅格和绍兴市行政边界shp看起来数据齐了但直接拖进ArcGIS往往不是黑屏就是范围对不上。原因多半出在坐标系不一致、NoData值被当成0、shp带了多余的属性或裁剪后数据量暴涨。这篇内容以这个zip为样本讲一条用GDAL命令行和Python脚本把“能看的DEM”变成“能用的DEM”的路径覆盖数据体检、边界裁剪、重投影和批量处理。适合做GIS数据处理、遥感制图、WebGIS服务发布的工程师也适合刚接触栅格数据的人照着步骤走一遍。2. 拆解DEM与shp数据包里的两个主角2.1 DEM数字高程数据栅格里的高度表DEMDigital Elevation Model本质是单波段栅格每个像元的值代表该位置的地面高程单位通常是米。绍兴市地处浙东丘陵和宁绍平原交界南侧会稽山、北侧平原高程范围从几米到几百米这决定了DEM数据里像元值分布跨度很大。常见的dem文件格式有GeoTIFF、IMG、ASCII Grid等这个zip里如果是.tif那就可以直接用GDAL处理。拿到dem文件后第一步不是打开看而是查元数据。用下面这个命令快速读出栅格的行列数、波段数量、坐标系和NoData值gdalinfo shengsha_dem.tif输出里重点看三层信息Size is后面的行列数决定数据量大小Coordinate System is后面是投影或地理坐标系统NoData Value决定后续计算时哪些像元应该被忽略。比如如果NoData是-9999那在做坡度或填洼时就必须带上这个值否则会把无效区域当成平地。2.2 shp文件用来“圈地”的矢量边界shp是ESRI Shapefile的简称但它不是一个文件而是至少由.shp、.shx、.dbf三个文件组成的矢量数据集。绍兴市区域范围shp可能是一个包含多个面要素的图层比如越城区、柯桥区、上虞区、诸暨市、嵊州市和新昌县等。这意味着裁剪前得先搞清楚这个shp是“一个整体范围”还是“多个分区”。用ogrinfo可以快速查看shp的几何类型、要素数量和属性表结构ogrinfo -al -so shengsha_boundary.shp选项-al表示列出所有图层-so表示只输出概要信息。结果中可以看到Feature Count、Extent、Geometry Column以及字段列表。如果Extent的坐标值在120-121左右、纬度在29-30左右那大概率是WGS84经纬度如果坐标值是几十万到三百万量级那就是投影坐标系比如UTM 51N或CGCS2000 3度带。2.3 坐标系和范围检查别让两个文件各说各话DEM和shp必须处在同一套坐标系下才能正确叠加。如果DEM是WGS84经纬度shp却是CGCS2000投影直接裁剪时GDAL会报错“Transform failed”或者输出一个空栅格。检查坐标系的简单方式是分别用gdalinfo和ogrinfo抓取信息然后对比。数据典型坐标系Extent示例DEMWGS84 / UTM 51N517xxx.xxx, 331xxxx.xxxshpCGCS2000 / 3-degree GK zone 3945xxx.xxx, 33xxxxx.xxx如果两者不一致常见做法是把矢量重投影到栅格坐标系而不是反过来。因为重投影栅格会引入像元重采样误差。重投影shp用ogr2ogr命令如下ogr2ogr -t_srs EPSG:32651 shengsha_boundary_utm.shp shengsha_boundary.shp这里-t_srs EPSG:32651指定目标坐标系为WGS84 UTM 51N。执行完再用ogrinfo确认Extent是否和DEM的对上。注意如果有多个要素这步不会改变要素个数只会改坐标值。检查完坐标系后再进入真正的裁剪环节。3. 用GDAL让shp“咬住”DEM裁剪与重投影实战3.1 解压和目录规划别让中文路径绊一跤Windows下解压这个zip时如果目录名包含中文GDAL在读写shp的.prj文件时偶尔会出编码问题。虽然GDAL 3.x对中文路径支持好了很多但为了保险建议把数据解压到纯英文路径比如D:\gisdata\shaoxing。解压后先看目录里有什么ls -la D:/gisdata/shaoxing通常能看到*.tif、*.shp、*.shx、*.dbf、*.prj。如果缺少.prj文件shp的坐标系信息就丢了后面裁剪时必须手工指定坐标系麻烦得多。所以先文件体检再动工具。3.2 用gdalinfo确认DEM尺寸和投影参数在裁剪前把DEM值和shp范围打印出来确认数据包和实际行政区是否匹配。一个更详细的命令gdalinfo -stats shengsha_dem.tif加上-stats会扫描像元值并输出最小、最大、平均值和标准差。这个信息能快速判断DEM有没有异常值。比如如果最小值是-9999那说明NoData没被识别如果最大值是9999也说明存在填充值。对应shp的范围则用ogrinfo -so -al shengsha_boundary.shp | grep -E Extent把输出的Extent和gdalinfo里的Upper Left、Lower Right对比。如果shp范围明显大于DEM裁剪后会有大片空白区域如果小于那就正常。确认无误后进入下一步。3.3 用gdalwarp执行带shp边界的裁剪最常见的裁剪方式是gdalwarp配合-cutline参数它把shp作为裁剪范围直接作用于栅格。命令可以这样写gdalwarp -cutline shengsha_boundary.shp -crop_to_cutline -dstnodata 0 -of GTiff shengsha_dem.tif shaoxing_dem_cut.tif这条命令的内部逻辑是读取原始DEM用shp的边界多边形做掩膜裁掉边界外的所有像元边界内部保留原始值。-crop_to_cutline表示把输出栅格的范围压到shp外接矩形上而不是保留原始DEM的完整范围。-dstnodata 0把裁剪后边界外的无效区域设成0这样后续查看时不会出现黑洞异常值。执行完后再跑一次gdalinfo查看新文件gdalinfo -stats shaoxing_dem_cut.tif这时Size的行列数应该远小于原始DEM而且NoData Value显示为0。如果输出文件大小反而变大多半是-tr分辨率参数没有跟随裁剪范围调整导致插值点变多后面会解释。4. 裁剪参数和常见坑从分辨率到NoData4.1 五个必调参数-cutline、-crop_to_cutline、-tr、-r、-dstnodatagdalwarp的裁剪参数看起来简单实际组合起来能改出四种不同效果。下表是实际动手时最常调的五个参数参数作用建议值坑点-cutline指定矢量裁剪边界shp路径不支持中文路径时读取失败-crop_to_cutline把输出窗口裁剪到shp外接矩形裁剪时必须加不加会保留原始范围文件巨大-tr设置输出像元分辨率与原始DEM一致设太大会丢细节设太小会内存爆炸-r重采样算法near或bilinear用nearest保真用bilinear平滑-dstnodata设置NoData值常用0或-9999不设会继承源值造成统计错误这里特别说一下-tr。如果DEM原始分辨率是12.5米裁剪后shp范围比原始DEM小那么输出行列数会随范围缩小而减少但如果-tr设置了比原始还小的值比如5米输出像元数量会激增四倍以上裁剪花的时间也成倍上涨。我一般会先用gdalinfo查原始像素大小然后让裁剪结果沿用这个值只在需要统一多源数据时才改-tr。4.2 坐标系不匹配导致输出全黑裁剪完成后输出文件存在但打开全是黑色或没有数据往往是因为shp和DEM坐标系不一致。gdalwarp内部默认会做坐标变换但如果shp缺少.prj文件它会假定shp和DEM坐标相同这时裁剪窗口可能落在数据范围之外。排查方法是先看裁剪输出文件的Extentgdalinfo shaoxing_dem_cut.tif | grep -E Lower Left|Upper Right如果Extent的坐标单位是米而原始DEM是经纬度那肯定不对。解决办法是先给shp设置正确的坐标系用ogr2ogr -s_srs EPSG:4326 -a_srs EPSG:32651重新赋值或者转换然后再执行裁剪。注意.prj文件缺失时ogrinfo会输出一个警告这个警告不要忽略。4.3 NoData值把山脉变成筛子很多DEM数据在和shp叠加裁剪时NoData会被保留下来。比如原始DEM的NoData是-9999裁剪后没设-dstnodata输出文件就会保留这个值。如果后续用gdal_calc或者Python做分析没有排除-9999会出现“山脉峰顶全部消失”的诡异结果。一个实际验证方法是统计裁剪后栅格的直方图gdalinfo -hist shaoxing_dem_cut.tif直方图里如果某个像元值出现特别多且正好等于-9999那就说明NoData没被正确改写。处理方式是在裁剪时明确指定-dstnodata -9999并且后续所有处理都带上-ignore选项或者干脆把NoData重设成0前提是真实高程不可能为0。4.4 处理shp里多个行政区划按属性分别裁剪绍兴市shp可能包含多个区县面要素如果想把整个shp裁成一整块前面命令已经够用。但如果要按每个区县分别输出DEM文件就需要按字段拆分。常见做法是用ogr2ogr先按属性导出单个shp再循环执行gdalwarp。比如按name字段提取柯桥区ogr2ogr -where name柯桥区 keqiao.shp shengsha_boundary.shp gdalwarp -cutline keqiao.shp -crop_to_cutline -dstnodata 0 shengsha_dem.tif keqiao_dem.tif注意-where后面的SQL表达式里中文字段值要用单引号包住。如果.dbf字段是GBK编码Windows下要加-lco ENCODINGUTF-8避免乱码。逐个手工写命令效率低所以下一步会给出批量脚本。5. 进阶用几何脚本批量出图和验证数据质量5.1 批量裁剪绍兴市各区县DEM既然shp里有多个区县边界就可以写一个Python脚本遍历要素并逐个裁剪。以下脚本用subprocess调用GDAL命令无需安装额外库import subprocess import ogr shp_path shengsha_boundary.shp ds ogr.Open(shp_path) layer ds.GetLayer() dem_file shengsha_dem.tif for feature in layer: name feature.GetField(name) out_shp f{name}.shp out_tif f{name}_dem.tif subprocess.run([ ogr2ogr, -where, fname{name}, out_shp, shp_path ], checkTrue) subprocess.run([ gdalwarp, -cutline, out_shp, -crop_to_cutline, -dstnodata, 0, dem_file, out_tif ], checkTrue)这里用了Python内置的ogr库遍历要素实际项目里如果没用过这个库可以用subprocess.run调用ogrinfo解析输出更简单但不够灵活。跑这个脚本前要确认Python环境已经安装了GDAL绑定的包比如wheel里的gdal版本。5.2 从DEM生成等高线图和坡度图裁剪不是终点很多人是为了做等高线或坡度分析。生成等高线可以用gdal_contour命令行像这样gdal_contour -a elev -interval 10 shaoxing_dem_cut.tif contour.shp-a elev是往输出shp里添加一个名为elev的属性字段存放每条等高线的高程值。-interval 10表示每隔10米画一条线。绍兴市北边平原地区如果也用10米间隔等高线会非常稀疏这时候可以改成5米或者3米。坡度图则用gdaldem生成gdaldem slope shaoxing_dem_cut.tif slope.tif -p-p表示输出坡度百分比而不是度。做这一步前务必确认DEM没有坑洞否则坡度图上会出现明显的条状噪声。5.3 验证裁剪结果让数字说话处理好数据后用gdalinfo -stats检查每个区县的DEM最小值和最大值就能判断是否裁偏了。比如柯桥区北部是平原最低海拔应该在3到10米左右南部靠近会稽山的地方最高海拔可能在500米以上。如果某个区县输出高程范围异常大或全为0就回头检查裁剪边界和NoData设置。另一个实用的验证方式是直接读一个像素值确认边界外是0、边界内是真实高程python -c from osgeo import gdal; dsgdal.Open(keqiao_dem.tif); print(ds.GetRasterBand(1).ReadAsArray(0,0,1,1)[0][0])这里读取左上角第一个像素如果输出是0说明落在了边界外如果是一个合理的高程值说明裁剪正常。这样五个区县逐一检查数据就可以放心交出去了。本文还有配套的精品资源点击获取
用GDAL和Python处理DEM数据：绍兴市高程栅格裁剪与坐标转换实战
简介浙江省绍兴市30米分辨率DEM数字高程数据包面向GIS从业者、城市规划与地理科研人员可用于地形分析、坡度计算、洪水模拟、视域分析等场景。压缩包共12个文件总大小约25.19MB核心为tif格式的高程栅格配套保存投影信息的tfw与金字塔优化文件ovr同时包含绍兴市行政边界shp矢量文件及其dbf属性表、prj坐标系、sbn/sbx空间索引、xml元数据等附件构成了一套可直接在ArcGIS、QGIS等软件中读取分析的标准地理数据集合。数据已覆盖绍兴市全域地形30米精度足以清晰呈现山丘、河谷等主要地貌特征辅助识别地质灾害风险区域、优化交通线路与城市空间布局。已有785人学习适合需要开展区域地形建模、环境评价或生成三维地形底图的相关人员使用。1. 拿到绍兴市DEM数据包先别急着拖进ArcGIS一个名为“浙江省绍兴市DEM数字高程数据含区域范围shp文件.zip”的压缩包拆开后是DEM栅格和绍兴市行政边界shp看起来数据齐了但直接拖进ArcGIS往往不是黑屏就是范围对不上。原因多半出在坐标系不一致、NoData值被当成0、shp带了多余的属性或裁剪后数据量暴涨。这篇内容以这个zip为样本讲一条用GDAL命令行和Python脚本把“能看的DEM”变成“能用的DEM”的路径覆盖数据体检、边界裁剪、重投影和批量处理。适合做GIS数据处理、遥感制图、WebGIS服务发布的工程师也适合刚接触栅格数据的人照着步骤走一遍。2. 拆解DEM与shp数据包里的两个主角2.1 DEM数字高程数据栅格里的高度表DEMDigital Elevation Model本质是单波段栅格每个像元的值代表该位置的地面高程单位通常是米。绍兴市地处浙东丘陵和宁绍平原交界南侧会稽山、北侧平原高程范围从几米到几百米这决定了DEM数据里像元值分布跨度很大。常见的dem文件格式有GeoTIFF、IMG、ASCII Grid等这个zip里如果是.tif那就可以直接用GDAL处理。拿到dem文件后第一步不是打开看而是查元数据。用下面这个命令快速读出栅格的行列数、波段数量、坐标系和NoData值gdalinfo shengsha_dem.tif输出里重点看三层信息Size is后面的行列数决定数据量大小Coordinate System is后面是投影或地理坐标系统NoData Value决定后续计算时哪些像元应该被忽略。比如如果NoData是-9999那在做坡度或填洼时就必须带上这个值否则会把无效区域当成平地。2.2 shp文件用来“圈地”的矢量边界shp是ESRI Shapefile的简称但它不是一个文件而是至少由.shp、.shx、.dbf三个文件组成的矢量数据集。绍兴市区域范围shp可能是一个包含多个面要素的图层比如越城区、柯桥区、上虞区、诸暨市、嵊州市和新昌县等。这意味着裁剪前得先搞清楚这个shp是“一个整体范围”还是“多个分区”。用ogrinfo可以快速查看shp的几何类型、要素数量和属性表结构ogrinfo -al -so shengsha_boundary.shp选项-al表示列出所有图层-so表示只输出概要信息。结果中可以看到Feature Count、Extent、Geometry Column以及字段列表。如果Extent的坐标值在120-121左右、纬度在29-30左右那大概率是WGS84经纬度如果坐标值是几十万到三百万量级那就是投影坐标系比如UTM 51N或CGCS2000 3度带。2.3 坐标系和范围检查别让两个文件各说各话DEM和shp必须处在同一套坐标系下才能正确叠加。如果DEM是WGS84经纬度shp却是CGCS2000投影直接裁剪时GDAL会报错“Transform failed”或者输出一个空栅格。检查坐标系的简单方式是分别用gdalinfo和ogrinfo抓取信息然后对比。数据典型坐标系Extent示例DEMWGS84 / UTM 51N517xxx.xxx, 331xxxx.xxxshpCGCS2000 / 3-degree GK zone 3945xxx.xxx, 33xxxxx.xxx如果两者不一致常见做法是把矢量重投影到栅格坐标系而不是反过来。因为重投影栅格会引入像元重采样误差。重投影shp用ogr2ogr命令如下ogr2ogr -t_srs EPSG:32651 shengsha_boundary_utm.shp shengsha_boundary.shp这里-t_srs EPSG:32651指定目标坐标系为WGS84 UTM 51N。执行完再用ogrinfo确认Extent是否和DEM的对上。注意如果有多个要素这步不会改变要素个数只会改坐标值。检查完坐标系后再进入真正的裁剪环节。3. 用GDAL让shp“咬住”DEM裁剪与重投影实战3.1 解压和目录规划别让中文路径绊一跤Windows下解压这个zip时如果目录名包含中文GDAL在读写shp的.prj文件时偶尔会出编码问题。虽然GDAL 3.x对中文路径支持好了很多但为了保险建议把数据解压到纯英文路径比如D:\gisdata\shaoxing。解压后先看目录里有什么ls -la D:/gisdata/shaoxing通常能看到*.tif、*.shp、*.shx、*.dbf、*.prj。如果缺少.prj文件shp的坐标系信息就丢了后面裁剪时必须手工指定坐标系麻烦得多。所以先文件体检再动工具。3.2 用gdalinfo确认DEM尺寸和投影参数在裁剪前把DEM值和shp范围打印出来确认数据包和实际行政区是否匹配。一个更详细的命令gdalinfo -stats shengsha_dem.tif加上-stats会扫描像元值并输出最小、最大、平均值和标准差。这个信息能快速判断DEM有没有异常值。比如如果最小值是-9999那说明NoData没被识别如果最大值是9999也说明存在填充值。对应shp的范围则用ogrinfo -so -al shengsha_boundary.shp | grep -E Extent把输出的Extent和gdalinfo里的Upper Left、Lower Right对比。如果shp范围明显大于DEM裁剪后会有大片空白区域如果小于那就正常。确认无误后进入下一步。3.3 用gdalwarp执行带shp边界的裁剪最常见的裁剪方式是gdalwarp配合-cutline参数它把shp作为裁剪范围直接作用于栅格。命令可以这样写gdalwarp -cutline shengsha_boundary.shp -crop_to_cutline -dstnodata 0 -of GTiff shengsha_dem.tif shaoxing_dem_cut.tif这条命令的内部逻辑是读取原始DEM用shp的边界多边形做掩膜裁掉边界外的所有像元边界内部保留原始值。-crop_to_cutline表示把输出栅格的范围压到shp外接矩形上而不是保留原始DEM的完整范围。-dstnodata 0把裁剪后边界外的无效区域设成0这样后续查看时不会出现黑洞异常值。执行完后再跑一次gdalinfo查看新文件gdalinfo -stats shaoxing_dem_cut.tif这时Size的行列数应该远小于原始DEM而且NoData Value显示为0。如果输出文件大小反而变大多半是-tr分辨率参数没有跟随裁剪范围调整导致插值点变多后面会解释。4. 裁剪参数和常见坑从分辨率到NoData4.1 五个必调参数-cutline、-crop_to_cutline、-tr、-r、-dstnodatagdalwarp的裁剪参数看起来简单实际组合起来能改出四种不同效果。下表是实际动手时最常调的五个参数参数作用建议值坑点-cutline指定矢量裁剪边界shp路径不支持中文路径时读取失败-crop_to_cutline把输出窗口裁剪到shp外接矩形裁剪时必须加不加会保留原始范围文件巨大-tr设置输出像元分辨率与原始DEM一致设太大会丢细节设太小会内存爆炸-r重采样算法near或bilinear用nearest保真用bilinear平滑-dstnodata设置NoData值常用0或-9999不设会继承源值造成统计错误这里特别说一下-tr。如果DEM原始分辨率是12.5米裁剪后shp范围比原始DEM小那么输出行列数会随范围缩小而减少但如果-tr设置了比原始还小的值比如5米输出像元数量会激增四倍以上裁剪花的时间也成倍上涨。我一般会先用gdalinfo查原始像素大小然后让裁剪结果沿用这个值只在需要统一多源数据时才改-tr。4.2 坐标系不匹配导致输出全黑裁剪完成后输出文件存在但打开全是黑色或没有数据往往是因为shp和DEM坐标系不一致。gdalwarp内部默认会做坐标变换但如果shp缺少.prj文件它会假定shp和DEM坐标相同这时裁剪窗口可能落在数据范围之外。排查方法是先看裁剪输出文件的Extentgdalinfo shaoxing_dem_cut.tif | grep -E Lower Left|Upper Right如果Extent的坐标单位是米而原始DEM是经纬度那肯定不对。解决办法是先给shp设置正确的坐标系用ogr2ogr -s_srs EPSG:4326 -a_srs EPSG:32651重新赋值或者转换然后再执行裁剪。注意.prj文件缺失时ogrinfo会输出一个警告这个警告不要忽略。4.3 NoData值把山脉变成筛子很多DEM数据在和shp叠加裁剪时NoData会被保留下来。比如原始DEM的NoData是-9999裁剪后没设-dstnodata输出文件就会保留这个值。如果后续用gdal_calc或者Python做分析没有排除-9999会出现“山脉峰顶全部消失”的诡异结果。一个实际验证方法是统计裁剪后栅格的直方图gdalinfo -hist shaoxing_dem_cut.tif直方图里如果某个像元值出现特别多且正好等于-9999那就说明NoData没被正确改写。处理方式是在裁剪时明确指定-dstnodata -9999并且后续所有处理都带上-ignore选项或者干脆把NoData重设成0前提是真实高程不可能为0。4.4 处理shp里多个行政区划按属性分别裁剪绍兴市shp可能包含多个区县面要素如果想把整个shp裁成一整块前面命令已经够用。但如果要按每个区县分别输出DEM文件就需要按字段拆分。常见做法是用ogr2ogr先按属性导出单个shp再循环执行gdalwarp。比如按name字段提取柯桥区ogr2ogr -where name柯桥区 keqiao.shp shengsha_boundary.shp gdalwarp -cutline keqiao.shp -crop_to_cutline -dstnodata 0 shengsha_dem.tif keqiao_dem.tif注意-where后面的SQL表达式里中文字段值要用单引号包住。如果.dbf字段是GBK编码Windows下要加-lco ENCODINGUTF-8避免乱码。逐个手工写命令效率低所以下一步会给出批量脚本。5. 进阶用几何脚本批量出图和验证数据质量5.1 批量裁剪绍兴市各区县DEM既然shp里有多个区县边界就可以写一个Python脚本遍历要素并逐个裁剪。以下脚本用subprocess调用GDAL命令无需安装额外库import subprocess import ogr shp_path shengsha_boundary.shp ds ogr.Open(shp_path) layer ds.GetLayer() dem_file shengsha_dem.tif for feature in layer: name feature.GetField(name) out_shp f{name}.shp out_tif f{name}_dem.tif subprocess.run([ ogr2ogr, -where, fname{name}, out_shp, shp_path ], checkTrue) subprocess.run([ gdalwarp, -cutline, out_shp, -crop_to_cutline, -dstnodata, 0, dem_file, out_tif ], checkTrue)这里用了Python内置的ogr库遍历要素实际项目里如果没用过这个库可以用subprocess.run调用ogrinfo解析输出更简单但不够灵活。跑这个脚本前要确认Python环境已经安装了GDAL绑定的包比如wheel里的gdal版本。5.2 从DEM生成等高线图和坡度图裁剪不是终点很多人是为了做等高线或坡度分析。生成等高线可以用gdal_contour命令行像这样gdal_contour -a elev -interval 10 shaoxing_dem_cut.tif contour.shp-a elev是往输出shp里添加一个名为elev的属性字段存放每条等高线的高程值。-interval 10表示每隔10米画一条线。绍兴市北边平原地区如果也用10米间隔等高线会非常稀疏这时候可以改成5米或者3米。坡度图则用gdaldem生成gdaldem slope shaoxing_dem_cut.tif slope.tif -p-p表示输出坡度百分比而不是度。做这一步前务必确认DEM没有坑洞否则坡度图上会出现明显的条状噪声。5.3 验证裁剪结果让数字说话处理好数据后用gdalinfo -stats检查每个区县的DEM最小值和最大值就能判断是否裁偏了。比如柯桥区北部是平原最低海拔应该在3到10米左右南部靠近会稽山的地方最高海拔可能在500米以上。如果某个区县输出高程范围异常大或全为0就回头检查裁剪边界和NoData设置。另一个实用的验证方式是直接读一个像素值确认边界外是0、边界内是真实高程python -c from osgeo import gdal; dsgdal.Open(keqiao_dem.tif); print(ds.GetRasterBand(1).ReadAsArray(0,0,1,1)[0][0])这里读取左上角第一个像素如果输出是0说明落在了边界外如果是一个合理的高程值说明裁剪正常。这样五个区县逐一检查数据就可以放心交出去了。本文还有配套的精品资源点击获取
