全球DEM数据shp包：小而轻的地形分析利器

简介全球DEM数据(shp文件不大很好用)是一份面向ArcGIS用户的数字高程模型资源覆盖全球陆地范围适用于地形渲染、坡度坡向分析、水文模拟和区域规划等基础GIS工作尤其适合需要快速获取全球尺度地形数据的初学者与从业者。压缩包共4个文件主数据为GeoTIFF格式的DEM栅格配套XML文件记录坐标系与元数据OVR金字塔文件加速分层显示整体大小仅18.82MB加载轻便且兼容性强。目前已有1622人学习/下载数据已做投影处理用户无需额外设置即可直接导入ArcGIS开展分析与制图。虽然分辨率并非极高但足以支撑教学演示、趋势判断和初步科研是解决全球DEM数据获取与预处理问题的实用资源。1. 全球DEM数据的这个shp包为什么“不大”反而是优点解压“全球DEM数据(shp文件不大很好用).rar”之后你得到的不是那种动辄几十GB的栅格瓦片而是一份能用QGIS或ArcGIS直接拖进去的矢量shp。标题里特意强调“不大挺好用”说白了就是这份数据已经被人处理过一轮把全球高程按点或简化多边形存了下来绕开了栅格数据在低配电脑上打开即卡死的毛病。它解决的是“我要一张全球尺度的地形图做宏观分析但不想要一台工作站”的需求适合做区域对比、课堂演示、插图底图和气候模型前期试探。严格说它不适合当工程级地形依据但作为一份全球DEM的轻量入口价值就在“看得见、拉得动、查得了”。2. 一份“不大”的全球DEM是怎么来的数据源血统与体积换算2.1 先判断血统这份shp来自哪个高程数据集拿到这份shp文件第一步不是急着打开而是确认它的原始血统。全球公开的DEM数据集就那么几个体系SRTM系列、ASTER系列、AW3D30、GLO-30这几类它们的原始分辨率都到30米左右差异在覆盖范围、空洞处理和局部噪声特征。数据集代号原始分辨率覆盖范围高程基准常见已知问题SRTM系列30m / 90m北纬60°到南纬56°之间EGM96大地水准面水体处理粗糙山区空洞ASTER系列30m全球EGM96局部条纹噪声与异常尖峰AW3D3030m全球EGM96山区细节较好裸岩区偶有伪等高线GLO-3030m全球EGM96综合质量最高但原始体量也最大判断这份shp血统的常规做法不是靠猜而是看三处压缩包里有没有附带的说明文档属性表里高程字段叫什么名字常见命名有height、elevation、alt、z_wgs84以及shp记录的空间坐标范围如果纬度最大到60°N左右那大多是SRTM派生数据如果覆盖全全球且负值异常少那多半来自ASTER或GLO-30的重采样版本。我一般会先看说明文档没有文档时看字段类型。有些打包者会把处理过程写进一个txt里包括原始分辨率、采用的简化阈值和投影参数这些信息直接决定了后续使用方式。如果纯粹只有几个shp文件那就只能靠属性表里的取值范围反推。2.2 shapefile存高程的三种组织方式点、面、线的取舍同样是“shp格式的DEM”实际几何类型可能完全不同。常见做法是三种组织方式点、面和等值线。点类型最常见每个点自身坐标加一个高程属性值用来做插值和统计最灵活面类型是一个个矩形网格多边形每个格网带一个高程值适合直接做专题填色但文件体积最大因为每个矩形至少4个顶点坐标跨越多边形边界时顶点数量还会翻倍等值线类型则不是用来做数值分析的它偏向制图表达属性表里通常是高程级别编号或具体高度。判断这份数据的几何类型很简单在QGIS图层列表里看“Geometry Type”或者在ArcGIS的图层属性里看要素类型就能确认。三种方式里我个人最推荐点类型因为用点做后续的重采样、插值和栅格化容错率最高。面类型最大的坑在于很多打包者对格网做了简化把相邻等高值合并成一个多边形这时候属性表里的面积字段可能与真实地形面积严重不符想恢复网格还得重新栅格化。2.3 “不大”的代价30米重采样到1公里体量缩到什么程度标题里“不大”两个字背后是一笔实打实的换算账。原始30米分辨率栅格一个像元覆盖约900平方米全球陆地面积约1.5亿平方公里意味着原始栅格有上百亿个有效像元这就是为什么完整版全球DEM拿普通电脑根本打不开。如果按每1公里取一个抽稀点来重采样每个点代表1平方公里全球陆地大约1.5亿个点。看起来还是不少但若按每10公里取一个点全球就只剩约150万个点配合shapefile这种紧凑的二进制存储体积就能压到百MB级甚至几十MB级。这个密度用来画一张A4大小的全球地形图、做大陆尺度的地势起伏统计已经够用但想看清楚某一个山谷的坡向和坡度就完全不行。所以说“不大”未必是坏事关键是你要清楚这份数据到底抽稀到了什么密度。打开属性表看要素总数再和上面这个换算对照一下就能反推出它的有效空间分辨率。3. 把这份全球小DEM跑起来QGIS快速体检与Python读取统计3.1 在QGIS里完成几何类型、坐标系和高程范围的快速检查拿到shp后第一件事不是做分析而是做体检。打开QGIS通过菜单“图层 添加图层 添加矢量图层”选中解压目录下的.shp主文件加载后做三步检查。第一步看“图层属性 源”确认坐标参考系。如果显示的是EPSG:4326说明是经纬度坐标这是全球DEM最常见的组织方式如果显示未知CRS说明压缩包里缺了配套的投影文件要手动补指定。第二步看几何类型和要素数量。打开属性表看左下角的总记录数。同时看几何是点还是面。如果属性表里高程字段排序后是一个个文本形态的数字比如“450”“1200”“80”那么字段类型极可能是文本型需要先转数值。第三步看高程值的分布范围。右键图层选择“属性 符号化”按高程字段做渐变配色观察图层是否呈现合理的从海平面到高原的过渡。如果大片区域颜色突变呈椒盐状说明原数据含有不少异常尖峰值后面必须做截断清洗。3.2 用GeoPandas快速读取并统计全局高程分布在Python里用GeoPandas读取这个shp比在GIS软件里点鼠标更可控。下面的代码把数据读进来先看几何类型和CRS再做基础统计。import geopandas as gpd # 读取shprows1000表示先只抽样1000行做体检 gdf_small gpd.read_file(global_dem.shp, rows1000) # 查看几何类型结果是Point、Polygon或LineString print(gdf_small.geom_type.unique()) # 查看字段列表定位高程字段 print(gdf_small.columns.tolist()) # 查看坐标参考系期望是EPSG:4326 print(gdf_small.crs)这里用rows参数做抽样读取是针对大文件的一个实用技巧。它不会真把全球所有要素读进内存而是先拿1000个样本判断数据形态避免一次性加载导致内存溢出。geom_type的输出决定了后续操作方向如果是Point直接插值和栅格化如果是Polygon需要先提取每个面的中心点或几何质心。确认结构后再把全量数据读入注意把高程字段强制转为浮点型。很多打包者在导出dbf时会把高程写成字符串不转换的话后面做统计和插值都白搭。import geopandas as gpd # 正式读取全量数据 gdf gpd.read_file(global_dem.shp) # 把高程字段转成浮点errorscoerce将无法解析的文本变NaN elev pd.to_numeric(gdf[elevation], errorscoerce) # 大范围经验截断只保留-500到9000米的记录 # 全球陆地最高点约8848米海底最深处约-11000米 # 超过这个范围的数值基本是异常噪声应直接过滤 valid elev[(elev -500) (elev 9000)] print(f有效记录数: {valid.shape[0]}) print(valid.describe())pd.to_numeric是处理“字段是字符串但内容像数字”的通用办法比直接astype更稳妥因为它遇到解析失败会产生NaN而不是抛异常。经验截断区间-500到9000基于全球地形的物理合理范围超过这个区间的值在绝大多数情况下都是传感器的坏值或重采样产生的边缘伪影。3.3 把点DEM插值回栅格线性插值最小命令如果你的最终目标是在某一区域画等值线或者做山体阴影渲染那需要先把这份点数据重新栅格化。对全球尺度的稀疏点最简单的做法是线性插值。import numpy as np from scipy.interpolate import griddata # 从GeoDataFrame中取出坐标和高程 coords np.array([(pt.x, pt.y) for pt in gdf.geometry]) values gdf[elevation].astype(float).values # 生成目标网格步长0.5度全球范围经度360列纬度180行 grid_lon, grid_lat np.meshgrid( np.arange(-180, 180, 0.5), np.arange(-90, 90, 0.5) ) # 线性插值适合宏观展示fill_value填-9999标记海洋无数据区 grid_elev griddata( coords, values, (grid_lon, grid_lat), methodlinear, fill_value-9999 ) # 把结果写出成GeoTIFF方便后续在GIS中叠加 from osgeo import gdal driver gdal.GetDriverByName(GTiff) ds driver.Create(global_dem_0.5deg.tif, 360, 180, 1, gdal.GDT_Float32) ds.GetRasterBand(1).WriteArray(grid_elev.astype(np.float32)) ds.GetRasterBand(1).SetNoDataValue(-9999) ds.SetGeoTransform((-180, 0.5, 0, 90, 0, -0.5)) ds.SetProjection(GCS_WGS_1984) ds.FlushCache()griddata是SciPy里做散点插值的标准函数method参数有三个选择nearest速度最快但结果呈块状linear平滑适中cubic连续性好但计算量大且容易在稀疏区域产生超调。对全球0.5度网格这种尺度linear是绝对不会出大错的选择。目标网格的步长不要小于点数据的实际密度否则插值出来的细节全是算法脑补的假象这一点特别容易被人忽略。4. 全球DEM shp数据的避坑清单坐标系、负值和字段类型的4次翻车4.1 坑一压缩包缺了.prj坐标系被误判为未知或WGS84现象把这个shp加载进QGIS后图层始终对不上底图要么整体偏移几十公里要么经纬度显示正常但面积变形离谱。属性表里坐标范围看似正常却无法和同一经纬度的地块匹配。原因shapefile是一个多文件格式.shp存几何.dbf存属性.prj存坐标系。压缩包里如果漏放.prjGIS软件就不知道这些坐标是经纬度还是米制投影坐标通常会按默认的未知CRS处理导致图层在页面上乱跑。解决先确认坐标值范围。如果X在-180到180之间、Y在-90到90之间基本可以断定是经纬度坐标手动给图层指定EPSG:4326即可。如果坐标是六到七位数的大数值那可能是原始栅格投影转换后残留的米制坐标需要找原始说明文档里的投影参数。最稳妥的办法是自建一个同名.prj文件右键图层 导出 保存要素为选择“与图层相同”的CRS导出后把生成的.prj文件拷回原目录。4.2 坑二属性表里高程字段看着是数字实际是字符串现象对高程字段做统计分级生成的图例顺序完全不对数值从小到大排列竟然按“1000、1100、1200、200、300”这种字典序排。画出来的剖面图毛刺成片像被人打乱过。原因很多转换工具把dbf字段类型定义成了字符型C而不是数值型N或F。字段名本身看不出区别但排序方式暴露了一切。字符串排序按首字符比较所以“200”会排在“1000”前面。解决在QGIS里打开字段计算器新建一个浮点型字段表达式写to_real(elevation)把转换结果放进去。用Python更简单读入后执行df[elevation] pd.to_numeric(df[elevation], errorscoerce)。关键是转换后马上用df.dtypes检查字段类型是否变成float。4.3 坑三负值到底是NoData还是海平面以下现象统计高程分布时发现大量负值数量多到不像海洋底部而且低到-9999、-32768这种整数值。如果直接按负数渲染陆地与海洋的边界会糊成一片。原因不同数据集处理NoData的约定不同。GLO-30用-32768标记无效区域ASTER系列常用-9999某些重采样工具用-32767而另一些版本用0表示陆地、负值表示海洋。打包者在做简化时如果没有统一掩膜两种标记会混在一起。解决先做一次频数统计如果某个负数值的数量比其他负值高出几个量级它极大概率是NoData掩膜。先过滤掉明显的NoData值再把余下负值视作海洋深度处理。最省心的做法是配合一份全球陆地边界多边形做空间相交把落在陆地多边形内的负值一律视为异常值剔除。4.4 坑四全球范围选择导致跨180度经线的要素“撕裂”现象做全球渲染时最东侧和最西侧的陆地被甩到图框两端中间隔着一道空隙典型的如太平洋区域亚欧大陆东岸和北美西岸在图上被割开。用GeoPandas做缓冲区或邻近分析时结果还会跨过整个半球。原因WGS84经纬度坐标的X范围是-180到180跨越180度经线的几何对象在数学上被判为“从经度179度跳到-179度”两个点的直线距离被错误地算成358度而不是2度。解决最稳妥的办法是按经度区间拆成两半处理。以0到360度西向为正则域将X值小于0的要素加360度得到0到360坐标范围再用适合全球的伪圆柱投影展示。另一种做法是专门用跨越太平洋为中心的投影坐标系QGIS里可以自定义一个以180度中央经线的等距圆柱投影专门用来承接这类数据集。5. 让这份数据更可信一个三层质量体检技巧接手的这类全球小DEM做一个三层体检再投入使用是我最近两年形成的习惯确实比直接出图省了后续返工的力气。第一层体检先看山体阴影中的异常尖峰。把抽稀后的点插值成栅格在QGIS里生成山体阴影渲染那些比周围高出几百米又孤立存在的亮点基本就是原始数据里残存的噪声点。处理办法是设置高程截断阈值将偏离局部中位数超过3倍标准差的点直接判为坏点在重采样阶段剔除。第二层体检做一次与参考数据的区域差值分析。这里不需要全量对比选择几个有代表性的区域比如一个高原区、一个沿海平原区、一个山区分别用这份shp与原始分辨率更高的DEM产品做高程差统计。差值中位数如果接近0说明整体无系统偏移如果中位数明显偏离0意味着这份数据的垂直基准可能被改动过后续坡度计算要格外小心。第三层体检验证抽稀密度是否足够支撑你的应用。一个实用办法是随机抽出一半点重新插值和原始全量插值结果做交叉验证两版在同一位置的高程差如果超过你业务要求的垂直精度说明点密度不够需要寻找更密的版本或者降级只用于宏观展示。这三层体检全过这份数据才能被当作可信的宏观地形参考。说实话这类“不大很好用”的数据包恰恰是项目里最容易因为体积小而跳过质量验证的环节可它一旦有问题影响会顺着下游每个环节扩散等分析都快做完了才发作那真是连后悔药都没得吃。现在你来接手这份数据先用半小时做一轮体检后面省下的是成倍的调试时间。希望帮到你。本文还有配套的精品资源点击获取