简介:全球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 | 局部条纹噪声与异常尖峰 |
| AW3D30 | 30m | 全球 | EGM96 | 山区细节较好,裸岩区偶有伪等高线 |
| GLO-30 | 30m | 全球 | 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 # 读取shp,rows=1000表示先只抽样1000行做体检 gdf_small = gpd.read_file("global_dem.shp", rows=1000) # 查看几何类型,结果是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") # 把高程字段转成浮点,errors="coerce"将无法解析的文本变NaN elev = pd.to_numeric(gdf["elevation"], errors="coerce") # 大范围经验截断:只保留-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), method="linear", 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存坐标系。压缩包里如果漏放.prj,GIS软件就不知道这些坐标是经纬度还是米制投影坐标,通常会按默认的未知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"], errors="coerce")。关键是转换后马上用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,意味着这份数据的垂直基准可能被改动过,后续坡度计算要格外小心。
第三层体检,验证抽稀密度是否足够支撑你的应用。一个实用办法是随机抽出一半点重新插值,和原始全量插值结果做交叉验证;两版在同一位置的高程差如果超过你业务要求的垂直精度,说明点密度不够,需要寻找更密的版本或者降级只用于宏观展示。
这三层体检全过,这份数据才能被当作可信的宏观地形参考。说实话,这类“不大很好用”的数据包,恰恰是项目里最容易因为体积小而跳过质量验证的环节,可它一旦有问题,影响会顺着下游每个环节扩散,等分析都快做完了才发作,那真是连后悔药都没得吃。现在你来接手这份数据,先用半小时做一轮体检,后面省下的是成倍的调试时间。希望帮到你。
本文还有配套的精品资源,点击获取