企业官网建设流程全解析

简介：这份广东省30米分辨率DEM高程数据面向GIS从业者、地理科研人员及规划决策者，用于地形分析、坡度坡向计算、排水方向模拟与自然灾害风险评估等场景。数据以Shapefile矢量格式组织，同时附带GeoTIFF栅格影像，可导入ArcGIS、QGIS等主流平台直接使用。压缩包共26个文件，涵盖shp、dbf、shx、prj等Shapefile核心组件，以及tif、tfw、ovr、aux等栅格与辅助文件，另有adf、xml元数据与投影信息，整体约53.97MB，结构完整、便于直接加载。目前已有1642人学习下载，适合需要广东省高精度地形底图的中高级用户，可支撑土地利用规划、基础设施选址与空间叠加分析等工作。

1. 广东省30米高清DEM高程数据：从瓦片拼接到地形分析的完整落地路径

做华南地区水文或选址分析的人多半遇到过这个场景：拿到一份区域边界，想算坡度、汇流或者淹没范围，第一步就卡在DEM上。公开的全球30米数据覆盖广东没问题，但直接下载下来往往是分幅的、带投影差异的、甚至高程基准不统一的碎片。广东省30米高清DEM高程数据这个方向，核心要解决的就是把覆盖广东全省的30米格网高程整合成一份可直接进GIS和脚本的连续栅格，再支撑坡度坡向、流域提取、通视分析这类具体任务。它适合做区域规划、遥感预处理、水文建模的从业者，也适合想用真实地形数据练手的学生。下面按数据获取、拼接裁剪、质量核查、地形派生、避坑、进阶验证的顺序讲透。

2. 数据获取与坐标系统一：先搞清楚你手里的是什么

2.1 30米DEM的常见来源与选型理由

覆盖广东的30米级DEM，从业者常用的来源有几类：一类是全球尺度的免费数据集，分辨率约1弧秒，在广东纬度上地面间距接近30米；一类是科研机构发布的区域产品，做过国内基准的校正；还有一类是测绘部门的基础地理信息成果，精度更高但获取门槛也高。选型时不要只看分辨率数字，要盯三个指标：水平基准、高程基准、垂直精度。水平基准决定你叠加矢量边界时会不会整体偏移，高程基准决定高程值和当地水准是否对得上，垂直精度决定你能不能做厘米级判断——30米数据做米级地形分析够用，做精细土方就不够。

我一般会先确认数据是地理坐标（经纬度）还是投影坐标（米）。广东跨了多个投影带，如果拿到的是分带投影数据，拼接前必须统一到同一个坐标系，否则接边处会出现错位。常见做法是统一到地理坐标WGS84或CGCS2000，再做后续处理；如果要做面积和坡度计算，再投影到适合广东的投影坐标系，比如UTM 49N/50N或高斯克吕格带。

2.2 用命令行批量下载与整理分幅数据

假设你已经从某个数据源拿到了分幅的GeoTIFF列表，第一步是批量下载并检查文件完整性。下面用Python写一个下载与校验脚本，逻辑是读取URL清单、断点续传、校验文件头。

import os import requests from osgeo import gdal # 分幅数据URL清单，实际使用时替换为你的数据源地址 url_list = [ "https://example.com/dem/tile_01.tif", "https://example.com/dem/tile_02.tif", ] save_dir = "./dem_tiles" os.makedirs(save_dir, exist_ok=True) for url in url_list: fname = os.path.join(save_dir, os.path.basename(url)) if os.path.exists(fname): # 已存在则跳过，避免重复下载 continue try: r = requests.get(url, stream=True, timeout=60) r.raise_for_status() with open(fname, "wb") as f: for chunk in r.iter_content(chunk_size=8192): f.write(chunk) except Exception as e: print(f"下载失败 {url}: {e}") continue # 用GDAL逐个打开，确认能被识别且不是空文件 for f in os.listdir(save_dir): path = os.path.join(save_dir, f) ds = gdal.Open(path) if ds is None: print(f"无法打开，可能损坏: {f}") else: print(f"{f} 尺寸: {ds.RasterXSize}x{ds.RasterYSize} 投影: {ds.GetProjection()[:60]}")

这段代码的关键点：stream=True避免大文件一次性读进内存；os.path.exists做断点续传的简化版；最后用GDAL打开每个文件，确认尺寸和投影信息。如果某个文件gdal.Open返回None，说明下载不完整或格式不对，需要重新下载。参数上，timeout建议设60秒以上，华南地区访问境外数据源可能较慢；chunk_size用8192是通用值，网络差可以降到4096。

2.3 坐标系统一与重采样参数

分幅数据如果投影不一致，拼接前要统一。用gdalwarp做批量重投影，目标坐标系选EPSG:4326（WGS84地理坐标）或EPSG:4490（CGCS2000地理坐标）。重采样方法选bilinear还是cubic？高程数据我一般用bilinear，它在连续表面上更平滑，不会像nearest那样产生阶梯；cubic在边缘可能过冲，产生异常高值。命令如下：

# 批量重投影到WGS84地理坐标，重采样用双线性 for f in ./dem_tiles/*.tif; do gdalwarp -t_srs EPSG:4326 -r bilinear -of GTiff \ "$f" "./dem_wgs84/$(basename "$f")" done

-t_srs指定目标坐标系，-r bilinear指定重采样，-of GTiff指定输出格式。注意：如果原始数据本身是地理坐标但基准不同（比如WGS84和CGCS2000），差异在米级，做区域分析可以接受，做高精度叠加就需要做基准转换，这需要七参数，一般从业者拿不到，所以选数据时尽量选基准一致的。

3. 拼接、裁剪与无效值处理：把碎片变成一份能用的栅格

3.1 用GDAL拼接分幅并处理接边

分幅数据统一坐标系后，用gdal_merge.py或gdalbuildvrt拼接。gdalbuildvrt更快，它先建虚拟栅格，不实际写文件，适合先预览；gdal_merge.py直接输出拼接结果。我一般先用VRT确认范围，再转成GeoTIFF。

# 第一步：建立虚拟栅格，把所有分幅按顺序列出 gdalbuildvrt ./dem_mosaic.vrt ./dem_wgs84/*.tif # 第二步：查看VRT信息，确认范围和分辨率 gdalinfo ./dem_mosaic.vrt | head -30 # 第三步：转成实际GeoTIFF，压缩用LZW节省空间 gdal_translate -of GTiff -co COMPRESS=LZW -co TILED=YES \ ./dem_mosaic.vrt ./dem_mosaic.tif

gdalbuildvrt会自动处理接边，但如果分幅之间有重叠且高程值不一致，重叠区会取最后一个文件的值。所以拼接前要检查相邻分幅在重叠区的高程差，差太多说明数据源有问题。-co COMPRESS=LZW是无损压缩，-co TILED=YES让后续按块读取更快，这两个参数对大数据量很实用。

3.2 按广东行政边界裁剪的两种方式

拼接完的全省数据可能包含周边省份，需要按广东边界裁剪。常见做法有两种：用gdalwarp -cutline直接裁剪，或者用rasterio的mask函数。前者适合命令行批量，后者适合Python流程。

# 用矢量边界裁剪，-crop_to_cutline精确到边界 gdalwarp -cutline ./guangdong_boundary.shp -crop_to_cutline \ -of GTiff -co COMPRESS=LZW \ ./dem_mosaic.tif ./dem_guangdong.tif

-cutline指定矢量边界，-crop_to_cutline让输出范围严格贴合边界外接矩形。注意：如果边界是经纬度坐标，要和DEM坐标系一致，否则裁剪结果会偏移。裁剪后边界外的区域会被设为无效值，通常是0或-9999，下一步要处理。

3.3 无效值与异常高程的清洗

裁剪后边界外是无效值，数据内部也可能有异常值，比如负高程（广东陆地最低点接近海平面，但不会大面积负值）或突变高值。用numpy加rasterio做清洗：

import rasterio import numpy as np with rasterio.open("./dem_guangdong.tif") as src: dem = src.read(1) profile = src.profile nodata = src.nodata # 把无效值、负值、超过2000米的值（广东最高约1900米）设为NaN mask = (dem == nodata) | (dem < -50) | (dem > 2000) dem_clean = dem.astype(np.float32) dem_clean[mask] = np.nan # 用3x3中值滤波去除孤立噪点，只对NaN以外的区域 from scipy.ndimage import median_filter valid = ~np.isnan(dem_clean) filled = dem_clean.copy() filled[valid] = median_filter(dem_clean[valid], size=3) # 写回新文件，NaN用-9999表示 filled[np.isnan(filled)] = -9999 profile.update(nodata=-9999, dtype=rasterio.float32) with rasterio.open("./dem_guangdong_clean.tif", "w", **profile) as dst: dst.write(filled, 1)

逻辑说明：先按阈值掩膜，广东最高峰约1900米，设2000米上限能滤掉异常高值；负值下限设-50米，因为广东有部分低于海平面的区域但不会太低。中值滤波只对有效区域做，避免把NaN扩散。参数上，size=3是3x3窗口，适合去除孤立噪点；如果噪点成片，要改用更大的窗口或先做插值。写回时nodata=-9999，方便后续GIS软件识别。

4. 从DEM到地形因子：坡度、坡向与汇流提取的实操

4.1 坡度与坡向计算的参数选择

DEM最常用的派生是坡度和坡向。用gdaldem可以直接算，参数里最关键的是-alg（算法）和-s（比例因子）。广东纬度约20到25度，地理坐标下1度纬度约111公里，1度经度约100到105公里，所以经度和纬度的地面距离不同，计算坡度时要考虑。常见做法是先投影到米制坐标再算，或者用gdaldem的-s参数指定比例。

# 先投影到UTM 49N（广东西部）或50N（广东东部），这里以49N为例 gdalwarp -t_srs EPSG:32649 -r bilinear \ ./dem_guangdong_clean.tif ./dem_utm49.tif # 计算坡度，单位度 gdaldem slope ./dem_utm49.tif ./slope.tif -alg Horn -s 1.0 # 计算坡向 gdaldem aspect ./dem_utm49.tif ./aspect.tif -alg Horn

-alg Horn是Horn算法，适合大多数地形；-s 1.0在投影坐标下比例是1，因为单位已经是米。如果在地理坐标下算，-s要设成纬度对应的米每度，比如-s 111120。坡向输出是0到360度，0为正北，顺时针增加。注意：平坦区域坡向是-9999或0，取决于实现，后续分析要排除。

4.2 汇流累积与河网提取的最小流程

水文分析常用汇流累积。用rasterio加richdem或pysheds可以做，但更稳的是用GDAL加TauDEM或WhiteboxTools。这里给一个用pysheds的最小示例：

import numpy as np import rasterio from pysheds.grid import Grid grid = Grid.from_raster("./dem_utm49.tif") dem = grid.read_raster("./dem_utm49.tif") # 填洼，去除内部凹陷 pit_filled = grid.fill_pits(dem) # 填平，处理平坦区域 flooded = grid.fill_depressions(pit_filled) # 计算流向，D8算法 flow_dir = grid.flowdir(flooded) # 计算汇流累积 acc = grid.accumulation(flow_dir) # 按阈值提取河网，阈值1000个栅格单元 streams = acc > 1000 with rasterio.open("./flow_acc.tif", "w", **grid.profile) as dst: dst.write(acc, 1)

逻辑说明：fill_pits和fill_depressions是水文分析的标准预处理，不填洼会导致流向中断；flowdir用D8算法，每个栅格流向8个邻居之一；accumulation统计每个栅格上游汇水单元数。阈值1000是经验值，30米数据下约0.9平方公里，广东湿润区河网密，阈值可以降到500。参数上，pysheds的flowdir默认D8，也可以选D-infinity，但D8更常用。

4.3 用坡度分级做选址筛选

坡度算出来后，按分级做筛选是常见落地。比如光伏选址要求坡度小于15度，可以用numpy做分级统计：

import rasterio import numpy as np with rasterio.open("./slope.tif") as src: slope = src.read(1) profile = src.profile # 分级：0-5, 5-15, 15-25, 25-35, >35 bins = [0, 5, 15, 25, 35, 90] labels = [1, 2, 3, 4, 5] slope_class = np.digitize(slope, bins) - 1 slope_class = np.clip(slope_class, 0, 4) # 统计各等级面积占比 total = slope_class.size for i, label in enumerate(labels): count = np.sum(slope_class == i) print(f"等级{label}: {count/total*100:.2f}%") profile.update(dtype=rasterio.uint8, nodata=255) with rasterio.open("./slope_class.tif", "w", **profile) as dst: dst.write(slope_class.astype(np.uint8), 1)

np.digitize按bins分级，np.clip确保索引在0到4。输出是uint8，节省空间。统计面积占比能快速判断区域地形特征，比如广东珠三角平原区0-5度占比高，粤北山区25度以上占比高。

5. 避坑与排查：30米DEM处理中最容易翻车的5个点

5.1 拼接后接边处出现明显台阶

现象：拼接后的DEM在分幅接边处有线性台阶，高程突变几米到几十米。原因：相邻分幅来自不同数据源或不同基准，高程值不一致。解决：拼接前检查重叠区高程差，超过5米就要做接边改正，常见做法是用重叠区均值做线性调整，或者只选同一来源的数据。

5.2 裁剪后边界外无效值被当成0参与计算

现象：算坡度时边界外出现大片0度或异常值。原因：裁剪后无效值设为0，而0在坡度计算中被当成有效高程。解决：裁剪时用-dstnodata -9999指定无效值，后续计算前用numpy把-9999设为NaN，所有派生计算都排除NaN。

5.3 地理坐标下算坡度结果偏大或偏小

现象：直接用经纬度DEM算坡度，结果和实际不符。原因：经纬度下x和y方向地面距离不同，gdaldem默认比例是1，导致坡度计算错误。解决：先投影到米制坐标再算，或者用-s参数指定纬度对应的米每度，广东约111120。

5.4 填洼后河网位置偏移

现象：填洼后提取的河网和实际河流位置对不上。原因：填洼改变了局部高程，尤其是平原区，导致流向改变。解决：填洼前先用高分辨率河流数据做引导，或者用breach算法代替填洼，pysheds支持breach_depressions，它只打通出口不改变内部高程。

5.5 大文件处理内存不足

现象：全省30米DEM拼接后文件几十GB，numpy读入直接内存溢出。原因：一次性读整个栅格。解决：用rasterio的窗口读取，分块处理，或者用gdal_translate先降分辨率做预览，确认流程后再全分辨率跑。分块大小建议512x512或1024x1024。

6. 进阶验证：用等高线回套与剖面检查DEM质量

6.1 用等高线回套验证高程一致性

DEM处理完，怎么确认没被改坏？我一般会生成等高线，和原始地形图或已知高程点对比。用gdal_contour从DEM生成等高线，再叠加到原始数据上目视检查。

# 从清洗后的DEM生成10米间隔等高线 gdal_contour -a elev -i 10 ./dem_guangdong_clean.tif ./contour.shp

-a elev指定高程字段名，-i 10是等高距。生成后可以在QGIS里叠加原始分幅或影像，看等高线是否贴合地形。如果某区域等高线密集但影像上是平地，说明该区域DEM有异常。

6.2 剖面检查：沿一条线看高程变化

另一个验证方法是取剖面。用rasterio沿一条线采样，画高程曲线，看是否有突变。

import rasterio import numpy as np import matplotlib.pyplot as plt with rasterio.open("./dem_guangdong_clean.tif") as src: # 定义剖面线起点和终点（行列号） row_start, col_start = 1000, 500 row_end, col_end = 1000, 2000 # 沿行方向采样 line = src.read(1, window=((row_start, row_end+1), (col_start, col_end+1))) profile = line[0, :] plt.plot(profile) plt.xlabel("采样点") plt.ylabel("高程 (米)") plt.title("DEM剖面检查") plt.savefig("./profile_check.png")

如果剖面出现垂直台阶或尖刺，说明该区域有拼接缝或噪点。正常地形剖面应该是连续变化的，除非有陡崖。参数上，采样线要跨过分幅接边，才能检查接边质量。

6.3 一个我常犯的错误：忽略垂直基准

最后说个血泪经验。早期做广东某区域淹没分析，DEM高程直接用了，结果算出来的淹没范围和实际差了几十米。后来才发现数据是椭球高，不是正常高，两者在广东差约10到20米。这个差异在平坦地区影响巨大。所以拿到DEM第一件事是确认高程基准，如果是椭球高，要用高程异常模型改正到正常高。没有改正模型就别做淹没分析，这是后悔药都买不到的坑。

验证方法上，找几个已知水准点，对比DEM高程和实际高程，差多少心里有数。如果差值是常数，可能是基准问题；如果是随机分布，可能是数据精度问题。这个检查花不了多少时间，但能避免后面所有分析白做。

希望帮到你。

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