企业官网建设流程全解析

简介：这份资源是四川省自贡市30米分辨率的DEM数字高程数据包，面向地理信息系统学习者、测绘与城市规划从业者以及开展地形分析的研究人员，可用于地形建模、洪水模拟、地质灾害评估和城市空间分析等教学与实践场景。压缩包共12个文件，约14.72MB，核心为自贡市dem.tif高程栅格，辅以ovr金字塔与tfw坐标文件便于快速浏览和配准；自贡市范围.shp、shx、dbf、prj、sbn、sbx构成完整Shapefile边界数据，另有多个xml元数据说明数据结构与属性。数据覆盖自贡市行政区域并延伸至部分邻近地区，便于处理边界效应或扩展研究范围。目前已有421人学习下载，适合需要真实区域高程底图、希望快速搭建GIS分析环境的初中级用户参考使用。

1. 自贡市 30m DEM 到手之后：先别急着往 GIS 里拖

拿到「四川省自贡市DEM数字高程数据30m（含本市级范围shp文件）.zip」这类数据包，很多人的第一反应是解压、拖进 QGIS 或 ArcGIS，然后直接出图。我见过太多人卡在这一步：要么坐标对不上，要么范围裁歪了，要么高程值全是异常。自贡这个地方地形很有代表性——西南丘陵区，海拔落差从两百多米到一千米出头，沱江穿城而过，低山、丘陵、平坝交错。30m 分辨率的 DEM 在这个尺度上刚好够用，既能看清丘陵起伏，又不会让文件大到跑不动。这个数据包的核心价值在于「本市级范围 shp 文件」——它给了你一个现成的行政边界，省去了自己去抠边界的麻烦。但边界文件和高程栅格能不能严丝合缝地对上，取决于坐标系、裁剪方式和 NoData 处理。这篇笔记就按我实际处理这类数据的顺序，把从解压到出可用成果的整条链路讲清楚，适合做水文分析、选址评估、地形可视化的从业者照着复现。

2. 解压后先看清三样东西：栅格、边界、坐标系

2.1 压缩包里通常有什么，怎么快速判断能不能直接用

这类数据包解压后一般包含三类文件：DEM 栅格文件（常见 .tif 或 .img 格式）、本市级行政边界 shp 文件（.shp/.shx/.dbf/.prj 一套）、可能还有一个说明文档或元数据文件。先别打开 GIS 软件，用命令行快速扫一遍文件结构和大小，能省不少事。

# 解压后进入目录，先看文件清单和大小 unzip -l "四川省自贡市DEM数字高程数据30m（含本市级范围shp文件）.zip" # 解压到指定目录 unzip "四川省自贡市DEM数字高程数据30m（含本市级范围shp文件）.zip" -d ./zigong_dem # 查看解压后的文件结构 find ./zigong_dem -type f | head -50 ls -lh ./zigong_dem/

逻辑说明：unzip -l只列出压缩包内容不解压，用来确认里面有没有你需要的 shp 和 tif。find加head是为了防止文件太多刷屏。ls -lh看文件大小，30m 分辨率下自贡全市范围（约 4381 平方公里）的 DEM 栅格通常在几十到一百多 MB 之间，如果只有几 MB，可能是被压缩过或者范围不全。

参数说明：-d指定解压目录，避免污染当前工作目录。如果压缩包内文件名含中文，Linux 下可能需要加-O GBK或-O CP936参数解决乱码，Windows 下用 7-Zip 或 Bandizip 更省心。

2.2 坐标系不核对，后面全白干

这是血泪经验里排第一的坑。DEM 栅格和 shp 边界文件如果坐标系不一致，你在 GIS 里看到的要么是边界飘到几千公里外，要么是栅格和边界错位。自贡市常用的坐标系有两种：地理坐标系 CGCS2000（经纬度，单位度）和投影坐标系 CGCS2000 3 度带高斯-克吕格投影（单位米，自贡大概在 104°E 到 105°E 之间，对应 3 度带第 35 带，中央经线 105°E）。

# 用 gdalinfo 查看 DEM 栅格的坐标系和范围 gdalinfo ./zigong_dem/zigong_dem_30m.tif | grep -E "Coordinate|Origin|Pixel|Upper|Lower" # 用 ogrinfo 查看 shp 文件的坐标系 ogrinfo -al -so ./zigong_dem/zigong_boundary.shp | grep -E "Extent|Layer|Geometry|Feature"

逻辑说明：gdalinfo输出的Coordinate System字段告诉你栅格用的什么坐标系，Origin和Pixel Size告诉你左上角坐标和像元大小。ogrinfo -al -so只输出摘要信息不打印所有要素，Extent字段显示边界的范围。把两者的范围对比一下，如果数量级差了几十倍，基本可以确定一个是经纬度一个是投影坐标。

参数说明：-so是 summary only 的意思，数据量大时必加，否则终端会被要素信息刷爆。如果gdalinfo输出里Coordinate System显示GCS_China_Geodetic_Coordinate_System_2000，说明是地理坐标系；如果显示CGCS2000_3_Degree_GK_Zone_35或类似带GK字样的，就是投影坐标系。

提示：如果两者坐标系不一致，不要直接改栅格或 shp 的坐标系定义（那是自欺欺人），要用gdalwarp或ogr2ogr做真正的坐标转换。

2.3 用 gdalwarp 统一坐标系的最小命令

假设 DEM 是地理坐标系，shp 是投影坐标系，或者反过来，统一到投影坐标系对后续面积、距离计算更友好。

# 将 DEM 从地理坐标系转换到 CGCS2000 3 度带 35 带投影 gdalwarp -t_srs EPSG:4544 -r bilinear -of GTiff \ ./zigong_dem/zigong_dem_30m.tif \ ./zigong_dem/zigong_dem_30m_proj.tif # 如果 shp 是地理坐标系，同样转换 ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4544 \ ./zigong_dem/zigong_boundary_proj.shp \ ./zigong_dem/zigong_boundary.shp

逻辑说明：gdalwarp的-t_srs指定目标坐标系，-r bilinear是重采样方法，DEM 连续数据用双线性插值比最近邻更平滑。ogr2ogr的-t_srs对矢量做投影转换。EPSG:4544 对应 CGCS2000 3 度带 35 带（中央经线 105°E），自贡市全域都落在这个带内。

参数说明：-r可选near（最近邻，适合分类数据）、bilinear（双线性，适合连续数据）、cubic（三次卷积，更平滑但可能过冲）。DEM 用bilinear是常见做法。如果转换后栅格出现条带或异常值，检查源数据是否有 NoData 值没设对。

3. 按行政边界裁剪 DEM：gdalwarp 和按掩膜提取怎么选

3.1 两种裁剪方式的本质区别

裁剪 DEM 到自贡市范围，常见做法有两种：一是用gdalwarp -cutline直接按 shp 边界裁，二是先用gdal_translate裁矩形范围再用gdalwarp按掩膜提取。前者一步到位但边界外像元会被设为 NoData，后者多一步但可控性更强。

我一般会先用-cutline试一次，如果边界复杂（自贡下辖四区两县，边界有飞地或狭长区域），再考虑分步处理。关键参数是-crop_to_cutline，不加这个参数的话，输出栅格范围还是原始范围，只是边界外变成 NoData，文件大小没减。

# 方式一：一步裁剪到自贡市边界 gdalwarp -cutline ./zigong_dem/zigong_boundary_proj.shp \ -crop_to_cutline -dstnodata -9999 \ -r bilinear -of GTiff \ ./zigong_dem/zigong_dem_30m_proj.tif \ ./zigong_dem/zigong_dem_30m_clip.tif # 方式二：先裁矩形再按掩膜提取（适合边界复杂或需要分县处理） gdal_translate -projwin 104.5 29.5 105.5 28.8 \ -of GTiff \ ./zigong_dem/zigong_dem_30m_proj.tif \ ./zigong_dem/zigong_dem_30m_rect.tif gdalwarp -cutline ./zigong_dem/zigong_boundary_proj.shp \ -crop_to_cutline -dstnodata -9999 \ ./zigong_dem/zigong_dem_30m_rect.tif \ ./zigong_dem/zigong_dem_30m_clip.tif

逻辑说明：-cutline指定裁剪边界 shp，-crop_to_cutline让输出范围贴合边界外接矩形，-dstnodata -9999把边界外像元设为 -9999，方便后续识别。方式二的-projwin参数顺序是xmin ymax xmax ymin，注意是「左上右下」不是「左下右上」，这个顺序搞反了裁出来是空的。

参数说明：-dstnodata的值要和后续分析工具兼容，ArcGIS 里常用 -9999，QGIS 里也可以用 -9999 或 0，但 0 在 DEM 里可能是真实高程（自贡最低点约 240m），所以别用 0。-projwin的坐标要跟栅格坐标系一致，投影坐标下单位是米，地理坐标下是度。

3.2 裁剪后必做的三项检查

裁完不是就完事了，至少检查三样：范围对不对、NoData 设没设对、高程值有没有异常。

# 检查裁剪后栅格的范围和 NoData 值 gdalinfo ./zigong_dem/zigong_dem_30m_clip.tif | grep -E "Upper|Lower|NoData|Size" # 统计高程值分布（排除 NoData） gdalinfo -stats ./zigong_dem/zigong_dem_30m_clip.tif | grep -A 20 "STATISTICS" # 用 Python 快速检查异常值 python3 -c " from osgeo import gdal import numpy as np ds = gdal.Open('./zigong_dem/zigong_dem_30m_clip.tif') band = ds.GetRasterBand(1) arr = band.ReadAsArray() nodata = band.GetNoDataValue() valid = arr[arr != nodata] print(f'有效像元数: {valid.size}') print(f'高程范围: {valid.min():.1f} ~ {valid.max():.1f} 米') print(f'均值: {valid.mean():.1f} 米') print(f'NoData 像元数: {arr.size - valid.size}') "

逻辑说明：gdalinfo -stats会计算栅格的统计信息，包括最小值、最大值、均值、标准差。自贡市高程范围大概在 240m 到 1000m 出头，如果统计出来最小值是 -9999 或最大值是 9999，说明 NoData 没设对或者有异常值。Python 脚本用numpy过滤 NoData 后统计，更直观。

参数说明：GetNoDataValue()返回栅格设置的 NoData 值，如果返回None说明没设，需要手动处理。arr != nodata做布尔索引，注意如果 nodata 是浮点数，直接用!=比较可能有精度问题，稳妥做法是用np.isclose或设一个容差范围。

注意：如果裁剪后有效像元数远小于预期（自贡全市约 4381 平方公里，30m 像元约 487 万个），检查是不是-crop_to_cutline没加，或者 shp 边界本身有问题（比如坐标系不对导致边界跑到别处）。

4. 避坑与排查：自贡 DEM 处理中最容易翻车的五个地方

4.1 坑一：shp 边界和 DEM 范围对不上，差了几十公里

现象：在 QGIS 里同时加载 DEM 和 shp，发现边界飘在栅格外面，或者只覆盖了栅格的一个角。

原因：最常见的是坐标系不一致。DEM 是地理坐标系（经纬度），shp 是投影坐标系（米），或者反过来。另一个可能是 shp 文件本身的范围就是错的，比如从某个在线地图下载的边界，坐标系标的是 WGS84 但实际是 GCJ02 偏移过的。

解决：先用gdalinfo和ogrinfo分别看两者的坐标系和范围。如果坐标系不一致，用gdalwarp和ogr2ogr统一。如果坐标系标称一致但范围还是对不上，用 QGIS 的「缩放到图层」功能分别看两个图层的实际位置，确认是不是数据本身有偏移。自贡地区如果用到从某些在线地图获取的边界，注意 GCJ02 偏移问题，需要做坐标纠偏。

4.2 坑二：裁剪后栅格全是 NoData 或者只有一条边有数据

现象：执行gdalwarp -cutline后，输出栅格大部分是 NoData，只有边缘一小条有数据。

原因：-cutline的 shp 坐标系和输入栅格坐标系不一致，gdalwarp 按 shp 的坐标去裁栅格，但两者不在一个空间参考下，导致裁剪区域错位。另一个可能是 shp 的几何类型有问题，比如是线而不是面，或者面有自相交。

解决：确认 shp 和栅格坐标系一致后再裁。用ogrinfo -al -so看 shp 的Geometry字段，必须是Polygon或MultiPolygon。如果是线，需要用ogr2ogr或 QGIS 转成面。面有自相交的话，用 QGIS 的「修复几何」工具处理。

4.3 坑三：高程值出现负值或异常大值

现象：统计高程时发现最小值是 -9999 或 -32768，或者最大值是 9999。

原因：-9999 通常是 NoData 值没被正确识别，-32768 是某些格式（如 ERDAS IMG）的默认 NoData。9999 可能是原始数据里的填充值或错误值。

解决：用gdalwarp的-dstnodata重新指定 NoData 值，或者在 Python 里用numpy过滤。如果原始数据本身就有异常值，需要用gdal_calc或 Python 做条件替换。

# 用 Python 将异常值替换为 NoData from osgeo import gdal import numpy as np ds = gdal.Open('./zigong_dem/zigong_dem_30m_clip.tif', gdal.GA_Update) band = ds.GetRasterBand(1) arr = band.ReadAsArray() # 将小于 0 或大于 2000 的值设为 NoData arr[(arr < 0) | (arr > 2000)] = -9999 band.SetNoDataValue(-9999) band.WriteArray(arr) ds = None

逻辑说明：自贡市真实高程不会低于 0 米也不会高于 2000 米，用这个范围做过滤是合理的。GA_Update以可写模式打开，SetNoDataValue设置 NoData 值，WriteArray写回。最后ds = None关闭数据集，确保数据落盘。

参数说明：阈值 0 和 2000 是根据自贡实际地形定的，如果换到其他地区要调整。写回前建议先备份原始文件，避免误操作覆盖。

4.4 坑四：裁剪后文件太大，跑不动

现象：自贡全市 30m DEM 裁剪后文件还有几百 MB，在 QGIS 里缩放卡顿。

原因：GeoTIFF 默认不压缩，30m 分辨率下像元多，文件自然大。另外如果 NoData 区域没有用掩膜或压缩，存储效率低。

解决：用gdal_translate加压缩参数重新输出。

# 用 LZW 压缩和金字塔重采样减小文件 gdal_translate -co COMPRESS=LZW -co PREDICTOR=2 \ -co TILED=YES -co BIGTIFF=IF_SAFER \ ./zigong_dem/zigong_dem_30m_clip.tif \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif # 添加金字塔（加速缩放显示） gdaladdo -r average ./zigong_dem/zigong_dem_30m_clip_compressed.tif 2 4 8 16

逻辑说明：COMPRESS=LZW是无损压缩，PREDICTOR=2对连续数据（如 DEM）压缩率更好，TILED=YES分块存储加速读取，BIGTIFF=IF_SAFER在文件可能超过 4GB 时自动用 BigTIFF 格式。gdaladdo添加金字塔，QGIS 和 ArcGIS 缩放时会自动调用，显示更流畅。

参数说明：PREDICTOR可选 1（不预测）、2（水平差分）、3（浮点预测），DEM 用 2 或 3 都行。金字塔层级2 4 8 16表示 1/2、1/4、1/8、1/16 分辨率，一般加到 16 或 32 就够。

4.5 坑五：用 ArcGIS 裁剪时结果和 gdalwarp 不一致

现象：同样的 shp 和 DEM，ArcGIS 的「按掩膜提取」和 gdalwarp 裁出来边界处像元值不一样。

原因：两者的重采样默认方法和 NoData 处理逻辑不同。ArcGIS 默认可能用最近邻，gdalwarp 默认用最近邻但可以指定双线性。边界处像元如果跨在裁剪线上，不同方法取值不同。

解决：统一重采样方法。gdalwarp 加-r bilinear，ArcGIS 里在环境设置中把「重采样技术」改为「双线性」或「三次卷积」。另外确认两者的 NoData 值设置一致。如果做定量分析，建议全程用同一套工具链，别混用。

5. 从 DEM 到可用成果：坡度坡向提取与水文分析的参数怎么设

5.1 坡度坡向提取：gdal 和 ArcGIS 的参数差异

DEM 最常用的衍生成果是坡度和坡向。自贡丘陵区坡度分析对农业选址、水土保持很有价值。用gdaldem命令行提取坡度，参数设置直接影响结果。

# 提取坡度（度为单位） gdaldem slope -of GTiff -compute_edges \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_slope.tif # 提取坡向 gdaldem aspect -of GTiff -compute_edges \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_aspect.tif # 提取山体阴影（可视化用） gdaldem hillshade -of GTiff -z 2 -az 315 -alt 45 \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_hillshade.tif

逻辑说明：gdaldem slope默认输出度，加-p输出百分比坡度。-compute_edges让边缘像元也参与计算，不加的话边缘一圈是 NoData。hillshade的-z 2是垂直夸张系数，自贡地形起伏不大，用 2 到 3 比较合适，-az 315是光源方位角（西北方向），-alt 45是光源高度角。

参数说明：坡度提取的算法基于 Horn 方法（3x3 窗口），ArcGIS 的「坡度」工具默认也是 Horn 方法，但 ArcGIS 会先做边缘填充。如果两者结果在边缘处不一致，是正常现象。坡向输出 0-360 度，0 为正北，90 为正东，ArcGIS 的坡向输出范围是 -1 到 360，-1 表示平地。

5.2 水文分析：填洼、流向、流量累积的关键参数

自贡有沱江及其支流，做水文分析前必须填洼，否则流向计算会断。

# 用 gdal_fillnodata 填洼（简单场景） gdal_fillnodata.py -md 10 -si 0 \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_dem_filled.tif # 用 Python + RichDEM 做更专业的填洼和流向分析 python3 -c " import richdem as rd dem = rd.LoadGDAL('./zigong_dem/zigong_dem_30m_clip_compressed.tif') dem_filled = rd.FillDepressions(dem, epsilon=True, in_place=False) flow_accum = rd.FlowAccumulation(dem_filled, method='D8') rd.SaveGDAL('./zigong_dem/zigong_flow_accum.tif', flow_accum) "

逻辑说明：gdal_fillnodata.py的-md 10是最大搜索距离（像元数），-si 0是搜索步长。RichDEM 的FillDepressions用epsilon=True做微填洼，避免大范围平坦区域导致流向不确定。FlowAccumulation用 D8 算法（八方向），适合丘陵区。

参数说明：填洼的-md参数根据洼地大小调，自贡丘陵区一般 10 到 20 够用。RichDEM 的epsilon填洼会在平坦区加微小梯度，保证流向唯一。流量累积结果中，高值对应河道，可以用来提取河网，阈值一般设 1000 到 5000 个像元（30m 分辨率下约 0.9 到 4.5 平方公里汇水面积）。

提示：如果做正式水文分析，建议用 ArcGIS 的 Hydrology 工具箱或 WhiteboxTools，RichDEM 适合快速验证。不同工具的填洼算法有差异，结果会有细微不同，选一套用到底就行。

5.3 用自贡 shp 边界做分区统计的实操

拿到坡度、坡向后，常需要按行政区统计。用zonal工具或 Python 的rasterstats库。

# 用 rasterstats 按自贡市边界统计坡度 from rasterstats import zonal_stats import geopandas as gpd boundary = gpd.read_file('./zigong_dem/zigong_boundary_proj.shp') stats = zonal_stats(boundary, './zigong_dem/zigong_slope.tif', stats=['min', 'max', 'mean', 'median', 'std'], nodata=-9999) print(stats)

逻辑说明：zonal_stats接受矢量边界和栅格路径，stats指定要计算的统计量。nodata参数要和栅格实际 NoData 值一致，否则会把 NoData 算进去。输出是一个列表，每个要素对应一个字典。

参数说明：如果边界有多个要素（比如自贡下辖的区县），结果会按要素顺序返回。all_touched=True参数会让边界接触到的所有像元都参与统计，默认只统计中心点在边界内的像元。做精确面积统计时用all_touched=True更合理。

6. 一个容易被忽略的技巧：用 DEM 做自贡市地形起伏度分级

地形起伏度是描述地表切割程度的指标，定义为特定窗口内高程最大值与最小值之差。自贡丘陵区用这个指标做地貌分级比单纯看坡度更直观。我一般用 Python 加scipy做滑动窗口计算，窗口大小取 3x3 到 15x15 像元，对应 90m 到 450m 的空间尺度。

from osgeo import gdal import numpy as np from scipy.ndimage import maximum_filter, minimum_filter # 读取 DEM ds = gdal.Open('./zigong_dem/zigong_dem_30m_clip_compressed.tif') band = ds.GetRasterBand(1) dem = band.ReadAsArray().astype(np.float32) nodata = band.GetNoDataValue() # 将 NoData 设为 NaN 避免影响计算 dem[dem == nodata] = np.nan # 计算 5x5 窗口（150m 尺度）的地形起伏度 window_size = 5 max_dem = maximum_filter(dem, size=window_size, mode='nearest') min_dem = minimum_filter(dem, size=window_size, mode='nearest') relief = max_dem - min_dem # 按起伏度分级（参考自贡实际地形） # < 30m 平坝，30-70m 浅丘，70-150m 中丘，> 150m 深丘/低山 classified = np.zeros_like(relief, dtype=np.uint8) classified[relief < 30] = 1 classified[(relief >= 30) & (relief < 70)] = 2 classified[(relief >= 70) & (relief < 150)] = 3 classified[relief >= 150] = 4 classified[np.isnan(relief)] = 0 # 输出分级结果 driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create('./zigong_dem/zigong_relief_class.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(classified) out_band.SetNoDataValue(0) out_ds = None

逻辑说明：maximum_filter和minimum_filter是scipy.ndimage提供的滑动窗口极值滤波，size=5表示 5x5 窗口，mode='nearest'让边缘像元用最近的有效值填充。relief就是起伏度。分级阈值 30m、70m、150m 是根据自贡实际地形定的——平坝区起伏度通常小于 30m，浅丘 30 到 70m，中丘 70 到 150m，深丘和低山大于 150m。

参数说明：窗口大小决定分析尺度，5x5 对应 150m 见方，适合村级规划；15x15 对应 450m 见方，适合乡镇级。分级阈值可以根据具体项目调整，但建议先统计relief的分位数（比如 25%、50%、75% 分位）再定阈值，比拍脑袋更靠谱。输出用GDT_Byte节省空间，NoData 设为 0。

这个方法的实用价值在于：自贡很多区域坡度看着不大，但起伏度很高，说明是破碎丘陵，修路和选址成本比平坝高得多。把起伏度分级图和 shp 边界叠加，能快速识别哪些乡镇以平坝为主、哪些以深丘为主。我一般会把这个结果和坡度图交叉，坡度大于 15 度且起伏度大于 70m 的区域标记为「建设难度高」，给规划部门做参考。

做这类分析最大的教训是：别一上来就追求花哨的算法，先把坐标系、NoData、裁剪范围这三样核对清楚，后面所有分析都顺。我见过太多人卡在坐标系上折腾一整天，最后发现只是 shp 的 .prj 文件缺失。拿到数据先gdalinfo和ogrinfo扫一遍，花不了五分钟，能省几小时。希望帮到你。

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