企业官网建设流程全解析

简介：汕头市30米分辨率DEM数字高程数据包，面向GIS分析、城市规划、灾害评估等场景，提供完整的地形栅格与行政边界矢量数据，可直接用于坡度计算、可视域分析、径流模拟、环境评估等专业应用。压缩包内共12个文件，涵盖TIFF栅格、Shapefile矢量、DBF属性表、投影定义文件等核心类型，并附空间索引、坐标参考和金字塔概览文件，保障数据在ArcGIS、QGIS等平台中的准确加载与快速浏览，整体大小约4.94MB。目前已有456人学习下载。数据以汕头市全域为范围，附带.shp边界矢量，且因以市界裁切也可能包含周边过渡区域，用户无需再自行裁剪和配准；配合30米地面分辨率，足以支撑中尺度地形分析与工程前期调研，尤其适合地理信息相关专业学生、科研人员和城乡规划从业者使用。

1. 拿到这份30m汕头DEM，先别急着解压

第一次看到“广东省汕头市DEM数字高程30m（含区域范围shp文件）.zip”时，大多数人只在意里面的高程栅格，直到在 ArcGIS 里打开才发现矩形范围比汕头市辖区大一圈，裁剪、统计、出图都得多做一步。这个包好就好在把汕头市行政边界拆成了标准 Shapefile，和 DEM 一起交付，省掉了从全国矢量库里抠边界的步骤。30 米分辨率意味着每个栅格覆盖约 900 平方米，在城市规划和中小流域分析里已经能看出地形起伏趋势，又比 12.5 米数据体量小，加载和计算都快得多。适合刚接触 DEM 的 GIS 从业者，也适合需要快速复现完整地形分析流程的人。

2. 拆包之前先认全文件：Shapefile 与 DEM 的完整拼图

2.1 压缩包里的十二个文件各管什么

解压后看到一堆后缀，新手容易只盯着汕头市DEM.tif，其余文件似乎无关紧要。实际上，没有汕头市范围.shp那组文件，你连准确的市域边界都要重新找；没有.tfw，栅格在高斯平面里的位置可能错位一整格。下面这张表把每个文件的角色列清楚。

shp实际上不是单个文件，而是shp + shx + dbf的最小组合，加上prj才能正确投影。这里还多了sbn/sbx，说明这个边界很可能是从 ArcGIS 环境里导出的。给这些文件做备份时，最好保留原压缩包，因为.sbn丢了以后 ArcGIS 会重新生成，但.shp一旦缺了dbf，属性表就保不住了。

2.2 .tfw、.aux.xml 和 .ovr：地理配准的底层约定

.tfw是 ESRI 世界文件，只有六行文本，分别代表 X 方向像元宽度、旋转项、旋转项、Y 方向像元高度、左上角 X、左上角 Y。GDAL 读取 TIF 时会优先使用内嵌的 GeoTIFF 标签，如果没有内嵌标签再读.tfw。这个包同时保留.tif和.tfw，说明数据可能在某个环节做过一次格式归一化。

.tif.aux.xml里的STATISTICS_MINIMUM和STATISTICS_MAXIMUM能提前告诉你高程范围，可以先用文本编辑器打开看一眼。.ovr金字塔只影响显示性能，不影响原始分析结果；如果觉得文件太大，删掉后 QGIS 会重新生成。常见做法是保留它，因为 30m 的 DEM 覆盖汕头全市后，缩放到县级视角时没有金字塔会明显卡顿。

2.3 用 gdalinfo 和 ogrinfo 检查数据质量

解压后建议做的第一件事不是直接拖进 ArcGIS，而是用 GDAL 命令行验证数据完整性。我一般会先进入解压目录，执行下面两条命令。

cd /path/to/shantou_dem gdalinfo 汕头市DEM.tif ogrinfo -so -al 汕头市范围.shp

gdalinfo输出里重点看四段内容：Size表示栅格行列数，Pixel Size是像元分辨率，Coordinate System是投影描述，NoData Value用于区分真实高程和无效区。ogrinfo -so -al中的-so是 summary only，只汇总不打印全要素；-al表示读取全部图层。输出中的Extent就是边界的左上右下坐标，和 DEM 的范围比对后，能确认边界是否完全落在 DEM 覆盖范围里。

如果发现中文文件名在终端里乱码，不是文件坏了，是字符编码问题。Windows 终端建议先chcp 65001切换 UTF-8，或者用英文文件夹名重新解压。GDAL 的 Python 绑定对中文路径的支持比命令行更友好，后面第 4 章会用到。

3. 用边界 shp 裁剪 30m DEM，并生成坡度、坡向与等高线

3.1 裁剪前先确认两个数据集的空间参考一致

直接用gdalwarp裁剪前，最好把上一步拿到的两个Extent摆在一起看。DEM 如果是经纬度坐标，Extent 单位是度，像元大小大约是 0.000277778；shp 如果是高斯投影，Extent 单位是米。这两者坐标系不同，不能直接叠加。用ogrinfo返回的Geometry: Polygon和GEOGCS判断，或者统一交给gdalwarp的-t_srs参数处理。

如果 shp 的坐标系和 DEM 不一致，常见做法是先给 shp 做投影转换，再裁剪。不过gdalwarp本身支持-cutline和-crop_to_cutline，它会在裁剪的同时完成重投影，所以直接执行下面的命令也能得到正确结果。

gdalwarp -cutline 汕头市范围.shp -crop_to_cutline -dstnodata -9999 -overwrite 汕头市DEM.tif shantou_clip.tif

-cutline指定边界矢量，-crop_to_cutline让输出范围严格按面边界生成而不是矩形包络，-dstnodata -9999把边界外的区域统一写成 -9999，后续做坡度、三维渲染时不会把空值当 0 米处理。如果不加这个参数，GDAL 默认继承原文件的 NoData 设置，但一旦原文件没有明确的 NoData，边界外很可能出现黑色斑块。-overwrite用于覆盖已有输出文件，避免重复运行报错。

3.2 坡度与坡向：单位选择比算法更关键

坡度计算是地形分析最常用的步骤，gdaldem提供了现成实现，不需要自己写滑动窗口。执行下面的命令：

gdaldem slope shantou_clip.tif shantou_slope.tif -p -s 111120 -of GTiff gdaldem aspect shantou_clip.tif shantou_aspect.tif -of GTiff

-p表示坡度输出为百分比而不是角度，如果你后续要按水利规范分级，通常用百分比阈值；不写-p时输出角度。-s 111120是水平和垂直方向的尺度比，只有 DEM 坐标系为度时才需要。它的含义是纬度方向 1 度约等于 111120 米，gdaldem slope需要让水平和垂直方向统一到同一单位才能算出真实坡度角。如果 DEM 已经是投影坐标系，就不要加-s，否则会把真实坡度缩小一半以上。

gdaldem aspect输出的坡向以正北为 0 度，按顺时针到 360 度。坡向本身没有单位争议，但要注意 NoData 区也会参与计算导致边界一圈出现异常坡向。解决办法是在计算前把 -9999 设置为 NoData，GDAL 在-dstnodata -9999时已经写入元数据，这里一般不会出错。

3.3 等高线生成与 shp 转 KML、TXT

等值线是拿 DEM 做图件时的“快速交付物”。gdaldem contour可以直接从裁剪后的栅格生成矢量等高线。

gdaldem contour -i 10 shantou_clip.tif shantou_contour.shp

-i 10表示等高线间隔 10 米，适合汕头这种低海拔地区；如果是丘陵区，间隔取 20 或 25 米更合适。生成出来的 shp 包含ELEV属性，记录每条线的高程值。注意gdaldem contour输出的 shp 会默认带上该字段，不需要手工拼接。

如果要把等高线分享给外业同事或放入其他 GIS 平台，常见做法是转成 KML 或者带表头的 TXT。ogr2ogr可以一条命令完成：

ogr2ogr -f KML shantou_contour.kml shantou_contour.shp ogr2ogr -f CSV shantou_contour.csv shantou_contour.shp

CSV 实际上就是“shp 转 txt”的一种变体：属性表被扁平化成行列格式，几何字段会丢。如果你需要每个顶点坐标，建议用ogr2ogr的-segmentize参数把曲线打散成点后再导出。这里给一个常用对比表。

这一套流程走完后，你会得到一个按汕头市边界裁剪的 DEM、配套的坡度、坡向和等高线。下一步可以做更有业务味道的分析，比如视域、汇水，但在此之前我通常先做一次高程统计，判断数据里是否存在不合理负值或空洞。

4. 用 Python 批量处理 30m DEM：高程统计、边界叠加与点位采样

4.1 为什么选择 rasterio + geopandas

GDAL 命令行适合单张处理，但当你需要跑几十个站点的节点高程，或者把 DEM 和道路 shp、行政区 shp 叠加做批量出图时，Python 更合适。rasterio负责读写栅格，geopandas负责矢量 io 和投影转换。这两者底层分别是 GDAL 和 Fiona，所以和前面命令行看到的数据行为完全一致。

汕头市的 DEM 大约几千乘几千像元，直接read(1)读入内存没问题。如果是 12.5m 或更高分辨率，建议用window参数按块读取，避免内存溢出。下面的代码先读取 DEM 基本信息并统计高程范围。

import numpy as np import rasterio with rasterio.open("汕头市DEM.tif") as src: elev = src.read(1) nodata = src.nodata res = src.res bounds = src.bounds if nodata is not None: elev = np.where(elev == nodata, np.nan, elev) print("分辨率: {}x{}, 像素数: {}".format(res[0], res[1], elev.shape)) print("高程范围: {:.2f} 到 {:.2f} 米".format(np.nanmin(elev), np.nanmax(elev)))

这里src.res返回 (dx, dy)，dy 通常是负值，但方向不影响分析；src.bounds是栅格矩形范围。把nodata替换为np.nan很关键：否则最小值会被 -9999 这类占位符污染。如果输出出现 -9999 作为最小值，说明 NoData 没有被正确识别。

4.2 裁剪到 shp 边界并叠加行政边界图

有了 DEM 数组，接下来用rasterio.mask做出和gdalwarp一样的结果，但保留在 Python 内存里，方便继续画图和采样。

import geopandas as gpd from rasterio.mask import mask with rasterio.open("汕头市DEM.tif") as src: shp = gpd.read_file("汕头市范围.shp", encoding="utf-8") if shp.crs.to_string() != src.crs.to_string(): shp = shp.to_crs(src.crs) clip_img, clip_transform = mask(src, shp.geometry, crop=True, nodata=-9999) clip_img = np.squeeze(clip_img) clip_img = np.where(clip_img == -9999, np.nan, clip_img)

shp.crs.to_string()和src.crs.to_string()比较的是 PROJ 字符串，如果一个是 EPSG:4490，一个是 EPSG:4618，字符不同但都是 CGCS2000 下的经纬度坐标，严格判断会触发投影转换。这里多一步转换并不会引入明显误差。mask(src, shp.geometry, crop=True, nodata=-9999)返回裁剪后的数组和新的仿射变换矩阵；裁剪结果是从原图行列中抽出的窗口，并不丢失空间参考。

把边界直接叠加在栅格图上，用matplotlib的extent参数把数组行列坐标映射到地理坐标即可。

import matplotlib.pyplot as plt from matplotlib import cm fig, ax = plt.subplots(figsize=(12, 10)) im = ax.imshow(clip_img, cmap=cm.terrain, extent=[clip_transform[2], clip_transform[2] + clip_transform[0] * clip_img.shape[1], clip_transform[5] + clip_transform[4] * clip_img.shape[0], clip_transform[5]]) gdf = shp.boundary.plot(ax=ax, color="black", linewidth=1.2) plt.colorbar(im, label="Elevation (m)") plt.savefig("shantou_dem_boundary.png", dpi=300)

extent四元组顺序是 左、右、下、上。注意clip_transform[5]是左上角 Y，clip_transform[4]是行方向像元高度，二者相乘后再加上首值，得到的是右下角 Y。如果直接拿bounds做 extent，在旋转栅格里会出错，但这里 30m DEM 的旋转项一般为 0，所以两种写法结果一致。

4.3 从 DEM 提取站点高程：shp 转 txt 的常见场景

地形分析常遇到“有一批站点坐标，要查这些位置的高程”的需求。用rasterio.sample可以在不逐像元读取的情况下完成采样。假设你有一份stations.csv，包含 lon 和 lat 两列。

import pandas as pd stations = pd.read_csv("stations.csv") points = list(zip(stations["lon"], stations["lat"])) with rasterio.open("汕头市DEM.tif") as src: vals = [x[0] for x in src.sample(points)] if nodata is not None: vals = [np.nan if v == nodata else v for v in vals] stations["elev"] = vals stations.to_csv("stations_elev.txt", sep="\t", index=False)

src.sample(points)需要点坐标和栅格坐标一致，如果你的站点是经纬度而 DEM 是高斯投影，必须先投影。最稳妥的方式是用geopandas.GeoSeries创建点图层，然后to_crs(src.crs)。这段代码导出的是制表符分隔的stations_elev.txt，既方便 Excel 打开，也和 GIS 里的“shp 转 txt”数据交换习惯一致。

4.4 批量处理多幅 DEM 的策略

如果后续从多个来源收集了相邻区域的 DEM 瓦片，比如从 OpenTopography 下载的 12.5m 数据，需要先做图幅拼接再裁剪。GDAL 提供一个vrt中间格式，不实际合并文件，只建立懒加载的索引视图。

gdalbuildvrt mosaic.vrt tile1.tif tile2.tif tile3.tif gdalwarp -cutline 汕头市范围.shp -crop_to_cutline mosaic.vrt shantou_merge_clip.tif

mosaic.vrt是 XML 文件，记录各瓦片的范围和分辨率。把多张 DEM 拼成 VRT 后，再用前面第 3 章的gdalwarp命令裁剪，这样不用生成中间大文件，节省磁盘空间。需要注意的是，不同来源的 DEM 高程基准可能差几米，拼接前用gdalinfo -stats检查重叠带的均值差，若差值大于 1 米，说明两个产品的高程参考不一致，不能直接混用。

5. 解压与打开 DEM 时的排错：zip 损坏、EOCD 缺失和高程异常的判断

5.1 zip 完整性测试与 invalid zip archive 修复

从网盘或旧硬盘拷出来的 zip 包，最常见的问题是下载截断，导致解压时报could not find EOCD。这个错误的本质是解压器在文件末尾找不到中央目录记录，文件不完整或结构被破坏。先不要急着用各种工具反复点“修复”，第一步是做完整性测试。

unzip -t 广东省汕头市DEM数字高程30m.zip

-t参数逐文件测试 CRC 校验，如果输出里有bad CRC或mismatch，说明某个文件已经损坏。此时用zip -FF把能恢复的部分抢救出来会损坏内部数据，不建议作为常规手段。最常见的可靠做法是重新下载原包，下载后先比对文件大小和源页面标注的字节数，一致再解压。

如果 zip 能解压但 GDAL 打开 TIF 报error read zip archive，先确认.tif是否被安全软件拦截或占用了别名为folded的随机盘符。移动文件到纯英文路径下再试一次，通常能解决 Windows 长路径或中文路径导致的文件访问异常。

5.2 DEM 显示全黑或范围不对时，检查 tfw 和 aux.xml

裁完的 DEM 在 QGIS 里显示全黑，先不要怀疑数据坏了，打开图层样式把渲染类型改成“单波段伪彩色”并重新计算最值。如果显示范围偏离真值，大概率是.tfw和.tif内嵌的 GeoTIFF 标签不一致。可以用gdalinfo查看输出的 Corner Coordinates，如果左上角坐标和 shp 的范围相距很远，就手动读取.tfw的前两行和后两行，确认像元尺寸和原点。

以下命令可以帮助重建丢失的.tfw，但前提是你知道栅格左上角坐标和像元大小。

gdal_translate -a_srs EPSG:4490 -a_ullr 116.2 23.6 117.2 23.0 汕头市DEM.tif corrected.tif

-a_ullr后面依次是左上经度、左上纬度、右下经度、右下纬度，这是应急方案，不能凭记忆填写，必须从辅助元数据或原始来源中获取，否则坐标系信息会被错误覆盖。

5.3 高程异常检测：用gdalinfo -stats或 Python 快速定位

高程最小值如果是 -9999、0 或负几万，说明 NoData 没有被正确屏蔽。先用gdalinfo -stats 汕头市DEM.tif查看统计信息和STATISTICS_VALID_PERCENT。如果有效像元占比低于 95%，说明大量区域是空值，边界裁剪、坡度计算都会产生异常。

Python 侧的判断更快：

import numpy as np import rasterio with rasterio.open("汕头市DEM.tif") as src: data = src.read(1) valid = data[data > -1000] print("有效像元: {}%".format(100 * valid.size / data.size)) print("5%-95%分位: {} {}".format(np.percentile(valid, 5), np.percentile(valid, 95)))

超过 5% 像元落在 -1000 以下，基本可以判定 NoData 掩膜丢失；5%-95%分位比最大值最小值更稳健，适合发现极端孤立异常点。对异常点位可以配合站点高程 CSV 输出做矩阵定位，再用 shp 边界叠加看是否落在汕头市辖区外。这一步做完，数据就能放心进入后续的填洼、视域或洪水淹没分析了。

本文还有配套的精品资源，点击获取