抚仙湖流域矢量边界与DEM高程底图数据制作全流程-育师

简介：这份资源面向从事流域分析、生态环境监测与水文地理建模的科研人员和GIS学习者，提供抚仙湖流域矢量边界及DEM高程的成套空间数据。包内共18个文件，约186.54MB，涵盖可编辑的ArcGIS MXD工程文件、标准Shapefile矢量边界、高精度DEM成图TIF，以及配套的prj投影文件、dbf属性表、tfw坐标信息和说明文档等，覆盖从数据加载到专题制图的主要环节。MXD工程预设了图层结构、投影坐标系与制图符号体系，保留矢量边界与DEM的空间配准关系，并配有图例、比例尺、指北针和标注模板，便于快速生成规范地图；Shapefile边界经空间校正与拓扑检查，采用WGS84坐标系，兼容ArcGIS、QGIS等主流软件，可用于流域分区统计与生态分区划定；DEM分辨率优于90米，支持坡度、坡向、地形起伏度计算及径流模拟、汇流路径提取等水文分析。目前已有104人学习，适合需要快速获取流域基础数据、开展空间分析与制图实践的用户。

1. 抚仙湖流域矢量边界与DEM高程：一套能直接进ArcGIS的底图数据怎么攒出来

做高原湖泊水文、面源污染或流域规划的人，多半都经历过这个场景：手头要一张抚仙湖流域的矢量边界，再配一套DEM高程，最好还能直接打开一个mxd工程文件，图层符号都调好，出图TIF拿来就能放进报告。听起来只是“找数据”，实际动手才发现——边界要么是整县行政界，要么是手工描的误差几百米；DEM要么是30米全球数据，湖岸线一带糊成一片；mxd打开全是红叹号，路径一换就崩。这套东西的核心不是“下载”，而是把矢量边界、DEM高程、mxd工程和成图TIF四样东西对齐到同一个坐标系、同一个范围、同一套符号里。适合谁？做流域水文分析的、写环评报告的、带学生做GIS课程的，以及需要快速复现一套标准底图的人。下面按我实际攒数据的顺序，把每一步的参数和坑讲清楚。

2. 抚仙湖流域矢量边界：从哪来、怎么裁、坐标系怎么定

2.1 边界数据的三个来源与取舍

抚仙湖流域的矢量边界，常见做法有三条路。第一条是水文分析法：用DEM做洼地填充、流向计算、流量累积，再在湖口位置设阈值提取汇水区。这条路最“物理”，但抚仙湖是断陷湖，周围山体陡峭，DEM精度不够时汇水区边界会跳来跳去。第二条是已有矢量裁剪：拿省市县行政界或三级流域分区，按抚仙湖位置手动修边。快，但行政界和真实分水岭不重合，湖东岸和湖西岸误差方向还不一样。第三条是遥感解译加人工修正：用高分影像或哨兵2号做湖岸线提取，再结合地形往上游推流域界。我一般用第二条做底，第三条修关键段，第一条做交叉验证。

提示：不管哪条路，最后都要落到一个面要素上，且必须是单部件、无自相交、无悬挂节点。拓扑检查不过，后面裁DEM必出碎多边形。

2.2 用QGIS或ArcGIS裁出抚仙湖流域范围

假设你手里已经有一份云南省或玉溪市的行政界，以及一份全国三级流域分区。目标是把抚仙湖所在的那个汇水单元单独提出来，再按湖岸线修边。下面以QGIS的Python控制台为例，ArcGIS里对应工具是“按位置选择”加“导出要素”。

# QGIS Python控制台：按位置提取抚仙湖所在流域面 from qgis.core import QgsVectorLayer, QgsProject, QgsSpatialIndex # 加载三级流域分区和抚仙湖中心点 basin_layer = QgsVectorLayer("data/basin.shp", "basin", "ogr") lake_point = QgsVectorLayer("data/fuxian_lake_center.shp", "lake_pt", "ogr") # 用空间索引找包含湖中心点的流域面 idx = QgsSpatialIndex(basin_layer.getFeatures()) for pt_feat in lake_point.getFeatures(): candidate_ids = idx.intersects(pt_feat.geometry().boundingBox()) for fid in candidate_ids: feat = basin_layer.getFeature(fid) if feat.geometry().contains(pt_feat.geometry()): # 选中并导出 basin_layer.select(fid) break # 导出选中要素为新shapefile from qgis.core import QgsVectorFileWriter writer = QgsVectorFileWriter.writeAsVectorFormat( basin_layer, "output/fuxian_basin_raw.shp", "UTF-8", basin_layer.crs(), "ESRI Shapefile", onlySelected=True )

这段逻辑是先定位湖中心点落在哪个流域面里，再把这个面单独导出。参数上注意：onlySelected=True只导出选中要素；编码用UTF-8避免中文属性乱码；坐标系先保持原样，下一步统一投影。裁完之后一定要打开属性表看面积，抚仙湖流域大致在1000平方公里量级，如果导出结果只有几十或上万，说明选错了层级。

2.3 坐标系选择：为什么必须用投影坐标

抚仙湖位于北纬24度附近，地理坐标系下1度经度约101公里，1度纬度约111公里。如果直接用WGS84地理坐标做面积统计和缓冲区分析，长度和面积都会变形。我一般统一到CGCS2000 3度带，中央经线102度，带号34。对应EPSG代码是4547（CGCS2000 / 3-degree Gauss-Kruger CM 102E）。在ArcGIS里就是“投影坐标系 > Gauss Kruger > CGCS2000 > 3-degree Gauss Kruger CM 102E”。所有矢量、DEM、mxd数据框全部设成这个坐标系，后面才不会出现图层错位。

# 用GDAL把边界投影到CGCS2000 3度带 ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4547 \ output/fuxian_basin_proj.shp output/fuxian_basin_raw.shp

-t_srs指定目标坐标系，源坐标系GDAL会自动从.prj文件读取。如果源文件没有.prj，需要先用-s_srs手动指定。投影后检查面积字段，和地理坐标下算出来的差异应该在合理范围内。

3. DEM高程数据：5米精度从哪来、怎么拼、怎么裁

3.1 别再到处找DEM了：三条获取路径的实测对比

热搜里那句“别再到处找dem数据了”戳中很多人。抚仙湖流域可用的DEM，常见就三类。第一类是公开全球数据：SRTM 30米、ALOS 12.5米、Copernicus 30米。免费，但湖岸陡崖处高程跳变明显，做水文分析时流向容易乱。第二类是国产高精度数据：ASTER GDEM 30米已算老将，资源三号或高分七号生成的5米DEM在部分区域可申请，但覆盖不连续。第三类是无人机或激光雷达：精度最高，但只适合小范围重点区，全流域拼接收费高、工作量大。我一般用12.5米做全流域底，重点湖滨带用5米数据替换。

注意：不同来源DEM的垂直基准可能不同，EGM96和EGM2008差几十厘米到一米。拼之前统一到同一基准，否则接边处会出现台阶。

3.2 用GDAL拼接与裁剪DEM的完整命令

假设你下载了多块12.5米DEM，覆盖抚仙湖流域。第一步拼接，第二步按边界裁剪，第三步重采样到统一分辨率。

# 1. 拼接多块DEM为一张 gdalbuildvrt -input_file_list dem_list.txt dem_mosaic.vrt gdal_translate -of GTiff -co COMPRESS=LZW -co TILED=YES \ dem_mosaic.vrt dem_mosaic.tif # 2. 按流域边界裁剪 gdalwarp -cutline output/fuxian_basin_proj.shp \ -crop_to_cutline -dstnodata -9999 \ -t_srs EPSG:4547 -tr 12.5 12.5 \ dem_mosaic.tif dem_fuxian_12m.tif # 3. 如果拿到5米数据，单独裁剪后替换重点区 gdalwarp -cutline output/fuxian_basin_proj.shp \ -crop_to_cutline -dstnodata -9999 \ -t_srs EPSG:4547 -tr 5 5 \ dem_5m.tif dem_fuxian_5m.tif

gdalbuildvrt先建虚拟拼接，不占磁盘；gdal_translate转成带压缩的GeoTIFF，COMPRESS=LZW无损压缩，TILED=YES方便后续按块读取。gdalwarp的-cutline按边界裁，-crop_to_cutline让输出范围紧贴边界，-dstnodata -9999设无效值，-tr指定输出分辨率。注意裁剪后边缘会有锯齿，因为栅格是方的、边界是弯的，这是正常现象，不是数据错误。

3.3 DEM重采样与空洞填补的取舍

如果原始DEM有空洞（云覆盖或数据缺失），常见做法是先用gdal_fillnodata.py填补，再重采样。但填补会引入平滑误差，水文分析时可能改变流向。我的习惯是：空洞面积小于1%就填，大于1%就换数据源。重采样方法上，从高精度到低精度用平均，从低精度到高精度用双线性，不要用最近邻，否则高程出现阶梯。

# 用rasterio检查DEM空洞比例 import rasterio import numpy as np with rasterio.open("dem_fuxian_12m.tif") as src: data = src.read(1) nodata = src.nodata total = data.size holes = np.sum(data == nodata) if nodata is not None else np.sum(np.isnan(data)) print(f"空洞占比: {holes/total*100:.2f}%")

这段代码读第一波段，统计等于nodata的像元数。如果nodata没定义，就用np.isnan。占比超过1%时，我会去检查是不是裁剪边界外的区域被算进来了，因为-crop_to_cutline之外的地方也是nodata，那不算空洞。

4. mxd工程文件：图层组织、符号化与路径修复

4.1 mxd里该放哪些图层、顺序怎么排

一个能直接用的抚仙湖流域mxd，图层顺序从下到上一般是：DEM高程（拉伸渲染）→ 流域边界（空心粗线）→ 湖岸线（蓝色实线）→ 河流（蓝色细线）→ 居民点（点符号）→ 注记。DEM放最底下做底图，边界压在上面才看得清。如果反过来，边界被DEM盖住，打开就是一片灰。

提示：mxd里的数据框坐标系必须和所有图层一致。如果某个图层是地理坐标，ArcGIS会动态投影，但出图时可能偏移。统一成EPSG:4547最稳。

4.2 用ArcPy批量设置符号和修复数据源

mxd最烦的是路径变了就红叹号。我一般用ArcPy在生成mxd时就写相对路径，或者用lyr文件替换数据源。下面这段是批量把图层数据源指到新路径并设符号。

# ArcPy：修复mxd数据源并设置DEM拉伸符号 import arcpy mxd = arcpy.mapping.MapDocument("output/fuxian_basin.mxd") df = arcpy.mapping.ListDataFrames(mxd)[0] df.spatialReference = arcpy.SpatialReference(4547) for lyr in arcpy.mapping.ListLayers(mxd): if lyr.name == "DEM": # 设置拉伸渲染，色带用地形 lyr.symbologyType = "RASTER_STRETCH" lyr.symbology.stretchType = "STD_DEVIATION" lyr.symbology.minimum = 1500 # 抚仙湖湖面约1720米，周边山体到2800米 lyr.symbology.maximum = 2800 elif lyr.name == "流域边界": lyr.symbologyType = "SIMPLE" lyr.symbology.symbol.outlineWidth = 2.0 lyr.symbology.symbol.outlineColor = arcpy.Color(0, 0, 0) mxd.save()

spatialReference设数据框坐标系；RASTER_STRETCH让DEM按拉伸显示，STD_DEVIATION比最小最大值更抗异常值；minimum和maximum按抚仙湖实际高程范围设，湖面约1720米，周围山顶到2800米左右。边界线宽2.0磅，黑色，出图清晰。

4.3 相对路径与打包：让mxd换台电脑也能开

ArcGIS的“地图包”功能会把数据和mxd打成一个包，但文件大。更轻的做法是设相对路径：在mxd属性里勾选“存储相对路径”，然后把mxd和data文件夹放同一级。这样换电脑只要整个文件夹拷过去就能开。如果还是红叹号，用arcpy.mapping.MapDocument的findAndReplaceWorkspacePaths批量替换。

# 批量替换mxd数据源路径 mxd = arcpy.mapping.MapDocument("output/fuxian_basin.mxd") mxd.findAndReplaceWorkspacePaths( r"\\old_server\gis\fuxian", r".\data" ) mxd.save()

findAndReplaceWorkspacePaths把旧路径替换成相对路径.\data。注意反斜杠和正斜杠在ArcPy里都能识别，但相对路径必须以点开头。

5. 标准成图TIF：出图参数、分辨率与配色避坑

5.1 出图TIF的分辨率与页面设置

成图TIF不是截图，是带地理坐标的栅格。在ArcGIS里用“导出地图”，格式选TIFF，分辨率设300dpi，勾选“写入世界文件”。页面大小按A3或A4，比例尺按流域范围定。抚仙湖流域南北约60公里，A3纸横放，比例尺1:250000左右能放下。如果设1:100000，图幅装不下，会裁掉南北两端。

# ArcPy导出成图TIF mxd = arcpy.mapping.MapDocument("output/fuxian_basin.mxd") df = arcpy.mapping.ListDataFrames(mxd)[0] df.zoomToSelectedFeatures() # 先缩放到流域边界 arcpy.mapping.ExportToTIFF( mxd, "output/fuxian_basin_map.tif", df, resolution=300, world_file=True, color_mode="24-BIT_TRUE_COLOR" )

resolution=300是印刷标准；world_file=True生成.tfw文件，让TIF带坐标；color_mode用24位真彩色，避免调色板丢失。导出后检查TIF的坐标系信息，用gdalinfo看有没有投影。

5.2 配色：DEM色带怎么选才不糊

DEM拉伸最怕两种：一是色带太花，山体纹理被颜色盖住；二是色带太灰，高低差看不出来。我一般用“地形”色带，低处绿、中间黄、高处棕白。抚仙湖湖面是平的，如果DEM没做水体掩膜，湖面会显示成低海拔绿色，和周围农田混在一起。解决办法是用湖岸线矢量把湖面区域裁掉，或者单独给湖面设一个蓝色填充。

注意：出图TIF如果用于报告印刷，色带要转成CMYK预览，否则屏幕上看是绿色，印出来偏蓝。

5.3 成图TIF的验证：坐标、范围、像素三查

导出后别急着用，先查三样。第一，gdalinfo看坐标系是不是EPSG:4547；第二，看范围是不是和流域边界一致，四角坐标有没有偏移；第三，看像素分辨率，300dpi下A3图幅大约4960×3508像素，如果只有几百像素，说明分辨率设错了。

gdalinfo output/fuxian_basin_map.tif | head -30

输出里重点看Coordinate System、Upper Left、Lower Right和Pixel Size。如果坐标系是空白的，说明导出时没写世界文件，需要回ArcGIS重新导出。

6. 避坑与排查：抚仙湖底图数据最常见的5个翻车现场

6.1 边界和DEM对不上，湖岸线偏移几百米

现象：mxd里边界压在DEM上，湖岸线位置明显偏西或偏东。原因：边界和DEM用了不同坐标系，或者边界是从地理坐标直接裁的没投影。解决：用arcpy.Project_management或ogr2ogr -t_srs统一到EPSG:4547，再重新裁DEM。检查方法是在湖岸线上取几个点，看DEM高程是不是接近1720米。

6.2 mxd打开全是红叹号，数据源找不到

现象：换电脑后mxd里所有图层前面有红色感叹号。原因：mxd存的是绝对路径，原路径不存在。解决：用findAndReplaceWorkspacePaths批量替换，或者提前设相对路径。如果数据在移动硬盘上，盘符变了也会红叹号，所以尽量用相对路径。

6.3 DEM裁剪后边缘出现锯齿状空洞

现象：裁剪后的DEM边缘像被狗啃过，有三角形空洞。原因：gdalwarp裁剪时边界外的像元被设成nodata，但边界是弯的，栅格是方的，边缘像元部分覆盖。解决：这是正常现象，出图时用边界线盖住边缘即可。如果做水文分析，先用gdal_fillnodata.py填边缘，但会轻微改变流向。

6.4 成图TIF没有坐标，放进报告对不上底图

现象：TIF在PS里能看，但拖进ArcGIS没有坐标，和边界叠不上。原因：导出时没勾选“写入世界文件”，或者用了截图工具。解决：回ArcGIS用ExportToTIFF重新导出，world_file=True。如果已经导出了，可以用gdal_translate补写坐标，但需要知道四角坐标。

6.5 5米DEM和12.5米DEM拼接处出现台阶

现象：重点区用5米数据，周边用12.5米，接边处高程突然跳变。原因：两个数据源垂直基准不同，或者分辨率不同导致重采样误差。解决：拼接前统一垂直基准到EGM2008，接边处用gdal_merge.py做羽化过渡，或者直接用12.5米全流域，5米数据只做参考不拼接。

7. 进阶：用Python脚本一键生成抚仙湖流域底图包

如果你经常要重复这套流程，写一个脚本把边界提取、DEM裁剪、mxd生成、TIF导出串起来，比每次手动点工具快得多。下面这个脚本框架用geopandas加rasterio加arcpy（ArcPy需在ArcGIS自带Python环境跑），核心是把前面几步的参数固化下来。

# 一键生成抚仙湖流域底图包（框架） import geopandas as gpd import rasterio from rasterio.mask import mask import arcpy import os # 参数区：改这里就能换流域 BASIN_SHP = "data/basin.shp" LAKE_PT = "data/fuxian_lake_center.shp" DEM_IN = "data/dem_mosaic.tif" OUT_DIR = "output" CRS = "EPSG:4547" # 1. 提取流域边界 basin = gpd.read_file(BASIN_SHP).to_crs(CRS) lake = gpd.read_file(LAKE_PT).to_crs(CRS) target = basin[basin.contains(lake.geometry.iloc[0])] target.to_file(f"{OUT_DIR}/fuxian_basin.shp", encoding="UTF-8") # 2. 裁剪DEM with rasterio.open(DEM_IN) as src: src_crs = src.crs geom = [target.to_crs(src_crs).geometry.iloc[0].__geo_interface__] out_image, out_transform = mask(src, geom, crop=True, nodata=-9999) out_meta = src.meta.copy() out_meta.update({ "driver": "GTiff", "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform, "nodata": -9999, "compress": "lzw" }) with rasterio.open(f"{OUT_DIR}/dem_fuxian.tif", "w", **out_meta) as dest: dest.write(out_image) # 3. 生成mxd（需ArcGIS环境） mxd = arcpy.mapping.MapDocument("CURRENT") # 或新建 df = arcpy.mapping.ListDataFrames(mxd)[0] df.spatialReference = arcpy.SpatialReference(4547) # 添加图层、设符号、导出TIF，略 mxd.saveACopy(f"{OUT_DIR}/fuxian_basin.mxd")

参数区里BASIN_SHP和LAKE_PT换成你自己的数据路径，CRS固定4547。第一步用contains找包含湖中心点的流域面，比空间索引更直接。第二步用rasterio.mask裁剪，crop=True让输出范围紧贴边界，nodata=-9999。第三步的mxd生成依赖ArcPy，如果没装ArcGIS，可以用qgis的Python API替代，但符号化设置更麻烦。