最近我想清楚了一件事用GIS做项目这么久栅格数据的空间分析始终是我的短板。不是不会点工具而是每次干活都靠“现学现卖”——今天要算坡度就打开坡度工具明天要提河网再临时刷一遍水文分析流程。工具倒是点得越来越熟可底层的原理一直模模糊糊。直到前阵子跟一位做规划的同行聊天他随口问我“你填洼的时候怎么判断要填到多深”我当场愣住。就是那一瞬间我决定专门用一个完整的日子把栅格数据的空间分析从头到尾认真复习一遍。这篇博文就是这个复习日整理出来的完整笔记。这篇内容不是零散技巧的堆积而是按照一个GIS从业者复习栅格空间分析时最顺手的组织方式来写的先搞清楚栅格与矢量在思维上的差异再把像元、分辨率、NoData这些地基概念重建一遍接着区分本地、邻域、区域三类运算随后深入到DEM表面分析与水文分析最后把实操里最容易翻车的几个坑和一条完整的汇水区划定案例全盘捋清。适合正在学GIS的学生、天天用ArcGIS/QGIS但没系统梳理过原理的从业者以及准备面试或技术方案前想临时补基础的人。1. 复习日的学习地图栅格空间分析到底在分析什么1.1 先搞懂栅格与矢量在思维上的差异我复习时做的第一件事不是打开软件而是把栅格和矢量两种数据模型的差异重新想了一遍。矢量的本质是“边界与路径”——用点、线、面描述离散的实体每个要素有独立的属性记录。栅格的本质则是“铺满空间的规则网格”——把一片区域切成无数个正方形像元每个像元只记录一个数值。所以栅格适合表达连续变化的场比如高程、气温、降水、土壤属性矢量适合表达有确定边界的实体比如地块、道路、行政区。这件事看似基础却直接决定了分析思路。比如要算一个地块的坡度用矢量数据很难处理因为坡度本质上是一个“邻域场”的概念必须借助DEM这类栅格通过相邻像元的高差来推算。反过来要统计某条河穿过哪些地块矢量叠加会比栅格逐像元分类更直观。很多从矢量入门的人总觉得栅格“反人类”其实只是因为没适应这种“没有边界只有网格”的思维方式。1.2 一份可以照着抄的复习清单既然是复习日就得按一条科学的路径走不能东一榔头西一棒子。我给自己列的清单是这样的数据准备与预处理投影统一、坐标范围裁剪、重采样、NoData处理、像元对齐。本地运算栅格计算器、数学变换、指数计算如NDVI、条件赋值。邻域运算焦点统计、坡度/坡向/山体阴影、滤波平滑。区域统计Zonal Statistics按分区聚合。专门分析表面分析、距离分析、水文分析。综合建模加权叠加、适宜性评价。这个顺序不是随便排的背后的逻辑是“每一步都是下一步的输入”。投影和重采样没做好后续所有叠加都会错位NoData没处理干净任何运算都会把边界污染掉不懂本地运算和邻域运算的区别就无法理解坡度工具为什么会输出那样一张图。所以复习必须从地基往上修而不是直接去看高阶工具。1.3 先跑算例、后补原理记得最牢我的复习方式有一个原则先动手跑出一个直观的结果再回头补原理。比如我先在一片DEM上做了填洼看到水流累积图之后再去查手册才真正明白填洼是为了消除DEM里的伪坑。如果反过来先看长篇大论的功能介绍再动手操作大概率是看完就忘。所以后面几节里的知识点我都会尽量结合具体的操作结果来讲而不是干巴巴地背定义。2. 像元、分辨率与NoData把“地基”重新打三遍2.1 分辨率不是越高越好分辨率就是单个像元对应的实际地面尺寸。30米分辨率的DEM一个像元代表地面上30米乘30米的区域2米分辨率的DSM则是一个2米乘2米的小格。很多人容易陷入“分辨率越高越精细”的误区。但分辨率越高像元数量按平方级增长处理时间和文件体积都会迅速膨胀。一片10公里乘10公里的区域30米分辨率大约11万个像元换成2米分辨率就是2500万个像元中间差了二十多倍计算耗时可不是线性增长那么简单。选择分辨率的核心依据是分析尺度。做全国尺度的地形粗略分析30米DEM完全够用做一个小流域的精细水文建模可能需要更高分辨率但也要考虑数据源本身的精度。如果你手里只有30米原始数据硬重采样成5米并不能增加真实信息只是把同一片区域切得更碎而已。这个道理就像数码照片插值放大不会真的补出细节。2.2 NoData是“非数值”不是0栅格数据里的NoData表示“此处无测量值”它可能来自传感器覆盖范围之外、云遮挡或数据拼接的空白区。问题在于不少初学者会不自觉地把它当成0来处理结果得到一堆诡异的结果。我举一个真实例子一份降雨栅格只覆盖了流域东部另一份DEM覆盖整个流域。如果直接用栅格计算器做某个公式运算边界外的NoData一旦被当成0计算结果的边缘会突然塌陷出现一条纵向的“零值带”。你拿到那张图很难看出到底是数据本身的问题还是公式出了问题排查起来非常痛苦。正确做法是在运算前先检查NoData的范围用条件表达式显式处理。比如需要把NoData替换成固定值时在ArcGIS里可以写Con(IsNull(dem), -9999, dem)需要把无效值直接排除时则用SetNull。处理完以后再看一遍栅格的属性统计确认最大值、最小值都在合理区间再进入下一步。2.3 数据类型与位深选择栅格的值类型也很关键。整型栅格适合存分类编号比如土地利用类型码、分区编号浮点型栅格适合存连续值比如高程、气温、NDVI。位深决定了文件大小和值域范围8位无符号只能存0到25516位有符号可以存-32768到3276732位浮点则适合带小数的高精度数据。双精度64位栅格很少用因为文件会大得离谱。做重分类时尤其要注意如果输出是浮点型分类结果会带着一堆无意义的半整数后续做面积统计时反而麻烦。所以我习惯在重分类工具的输出设置里强制指定整型后面的Zonal统计会更干净。3. 三大运算类型本地运算、邻域运算、区域统计的取舍这是整个复习日里价值最高的一部分。栅格空间分析的操作成百上千但底层的运算逻辑其实只有三大类本地运算、邻域运算、区域统计。把这三类搞清楚了以后再遇到什么新工具你都能快速归位。3.1 本地运算逐像元做“一对一或一对多”的函数映射本地运算指的是输出像元值只由同一个位置的输入像元值决定完全不考虑周围环境。最典型的应用就是栅格计算器里的四则运算、三角函数、条件判断以及NDVI这类指数计算。NDVI在QGIS栅格计算器里通常这么写(NIR - Red) / (NIR Red)它逐像元读取近红外波段和红波段的反射率算出一个在-1到1之间的值用来反映植被覆盖状况。因为整个过程不依赖像元邻居这就是标准的本地运算。本地运算最常见的技术要求是输入栅格必须完全对齐——坐标范围一致、分辨率一致、像元网格的边严格重合。如果没有对齐两个栅格在同一位置像元对应的地面范围并不一致虽然工具通常不会报错但结果已经悄悄失真了。所以我每次做多波段叠加之前都会把数据统一重采样到同一分辨率并设定一个基准栅格做对齐参考。3.2 邻域运算像元加上邻居一起参与计算邻域运算的核心思想是输出值不仅由中心像元决定还由它周围的邻居像元共同决定。最典型的就是焦点统计用一个3×3、5×5或者更大的窗口扫过整个栅格对窗口内的像元求平均值、最大值、众数或标准偏差。坡度、坡向、山体阴影这些地形表面分析工具本质上也都是邻域运算。坡度就是通过中心像元与相邻像元的高差拟合出一个局部斜面然后计算这个斜面相对水平面的倾角。公式上可以简化为坡度度 atan(sqrt(dz/dx^2 dz/dy^2)) * 180 / π其中dz/dx和dz/dy就是通过3×3窗口推算的南北方向和东西方向的高程变化率。邻域运算非常吃性能。窗口每扩大一倍参与计算的像元数会大幅增加大范围的高分辨率影像做一次9×9滤波都能卡到怀疑人生。如果只是想去噪平滑先试试3×3窗口追求更大范围的平滑效果更好的思路是先重采样降分辨率或者做分块并行计算而不是盲目把窗口调大。3.3 区域统计在指定分区内做聚合区域统计Zonal Statistics解决的问题是在某个空间划分单元内另一个值栅格的平均值、总和、极值、标准差是多少。最简单的应用一个vector面图层划定了若干流域边界一张栅格是多年平均降雨量你不需要把每个像元单独导出来再求均值直接用Zonal Statistics工具就能得到每个流域内的平均降雨量、降雨总量。这为后续的水资源总量估算提供了直接输入。做区域统计时要留意分区栅格和值栅格最好在同一分辨率下。如果分区是矢量面建议在工具参数里指定“以栅格模拟为准”或者把矢量先转成栅格否则工具会按像元中心点落在哪个面里来判断归属边缘像元的归属可能会和你想的不一样。3.4 先判断再选择一张对比表说清楚差异复习到这一步我把三类运算拉了一张对比表随时可以查阅考察角度本地运算邻域运算区域统计输出值依赖同一位置的输入值中心像元及周边邻域区域内所有像元聚合典型用途NDVI、条件赋值、数学变换坡度、坡向、滤波流域平均降雨、分区统计常配工具栅格计算器、Map Algebra焦点统计、坡度工具Zonal Statistics常见误区输入栅格没有对齐窗口大小拍脑袋分区与值栅格分辨率不一致我的口诀是先问结果看重“单个像元、周边像元、还是区域集合”再把问题归类到三类运算中选工具基本不会跑偏。4. DEM表面分析与水文分析最常用也最容易做错的模块4.1 表面分析工具背后的物理含义表面分析是我平时用得最多的模块其中坡度、坡向、山体阴影、曲率是最基本的四个。坡度的输出有两种单位度数和百分比。前者适合表达斜坡倾角后者适合工程中的坡率计算比如“坡度5%”意味着水平距离每100米抬升5米。工具里默认输出什么单位不同软件不一样出图前一定要确认。坡向输出0到360度的方向角表示坡面朝向哪个方位北为0度顺时针递增。它能用来做日照分析、干燥度评估。山体阴影则依赖太阳方位角和太阳高度角两个参数参数的调整会彻底改变画面起伏的视觉。很多人把它当成普通渲染其实它是用邻域运算模拟太阳照射下的阴影是可视化成果里最出效果的工具之一。4.2 水文分析第一步为什么一定是“填洼”水文分析是一套固定流程填洼Fill、流向Flow Direction、流量累积Flow Accumulation然后再提取河网和集水区。填洼这一步经常被跳过或者敷衍对待。DEM中的洼地是指那些高程低于周围但实际并不蓄水的像元可能是原始数据噪声或插值造出来的伪坑。如果直接计算流向水流在这些伪坑面前无法找到更低的下游像元流向就会中断流量累积图里会出现大量不连贯的零值区域。Fill工具做的工作就是把洼地抬高到能让水流继续外溢的高度。但填洼并不是越彻底越好。在喀斯特地貌或者人工水库区域真实洼地是存在且应该保留的过度填洼会把这些真实地形抹掉导致生成的水系在图纸上“凭空穿过”一座水库或汇水洼地。所以在自动填洼之后最好对比原始DEM检查填挖深度分布用阈值控制填洼范围。4.3 D8与D-infinity流向算法的两种思路流向计算的核心是判断每个像元的水流方向。最常用的D8算法只允许水流流向8个邻域中坡度最陡的那一个像元。它简单、稳定流域边界划定这类应用基本都够用但在相对平缓的地区D8会生成大量彼此平行的直线河道看起来不太自然。另一种D-infinity算法会把水流按照坡向比例分配到多个邻域像元更接近实际扩散的径流路径。如果做的是地形湿度指数、污染扩散模拟这类需要考虑水流分散的模型D-infinity更合适如果只是划定汇水区和流域边界D8仍然是最稳妥的选择。4.4 河网提取阈值不要照抄默认值大部分教程会告诉你流量累积超过某个阈值就能当作河流。但这个阈值不是固定的它取决于分辨率、地形和当地径流特征。流域尺度大、地形起伏明显阈值可以设高一些地形平缓、研究区域又小阈值就得调低。我的经验是直接用试错法在流量累积图上依次尝试阈值100、500、1000、5000把每次生成的河网矢量和卫星影像或地形图上的实际水系叠加比对看哪一版和实际沟道最接近。这个方法虽然看起来“笨”但比查任何网上的默认参数都可靠。阈值设得太小河网密得像毛细血管太大会把实际的季节性冲沟全部丢掉后续汇水区面积就会偏小。5. 栅格实操里的翻车现场投影、分辨率、NoData与内存复习日最有价值的部分其实是复盘踩过的坑。下面这些坑我敢说大部分经常做栅格分析的人都遇到过。5.1 坡度计算前先把投影统一有一次我拿到一份WGS84经纬度坐标系的DEM没多想就直接做了坡度计算结果出来的坡度值大到离谱整片山坡忽高忽低完全无法使用。原因很简单经纬度坐标以度为单位而高程以米为单位两者根本不在同一量纲上坡度公式在混合单位下会计算出错。所以现在我的习惯是只要涉及坡度、坡向、距离、面积的栅格分析一律先把数据转换到合适的投影坐标系比如UTM、Albers等轴等积投影再做计算。经纬度数据只适合做定位和范围可视化不适合做量算。5.2 分辨率不匹配与像元未对齐把0.2米分辨率的无人机正射影像和10米分辨率的DEM叠加时因为两者的像元网格完全不在同一套“棋盘”上叠加结果会出现大量局部拉伸后的脏数据。软件通常不会给你红色报错但结果图一眼就能看出来不对。解决方法是统一重采样并设置一个基准栅格Snap Raster。在ArcGIS的环境设置里指定好基准栅格后后续所有分析输出都会自动对齐到同一套像元网格。这一步看着不起眼却能避免后来的一大堆麻烦。重采样方法也要选对分类数据土地利用类型码必须用最近邻法因为插值会产生不存在的类别值连续数据高程、NDVI可以用双线性或三次卷积画面更平滑。这个原则我在第6章案例里会再验证一遍。5.3 边界出现诡异的零值带先查NoData如果你发现计算结果边缘有一片不自然的低值区先别急着怀疑公式检查一下NoData是不是被当成了0参与运算。我处理夜光遥感数据和气象栅格时经常踩这个坑。排查方法很简单用识别工具点开边界像元看它的值到底是什么再用栅格统计里的像元总数和被有效值覆盖的数量做比较。处理方式也很直接先判断NoDataSetNull(IsNull(dem), dem)或者把它替换成一个明显不可能出现的哨兵值比如-9999后续再统一处理。5.4 大数据栅格的内存不足分块是正经出路一次对一个5万乘5万像元的较大影像做3×3均值滤波我直接把整个栅格读进内存结果机器卡到几乎死机。之后再做大范围栅格处理我都改用分块思路把栅格按512×512或1024×1024的窗口切块逐块计算再拼接结果。用Rasterio可以很轻便地做到import rasterio from rasterio.windows import Window with rasterio.open(big_dem.tif) as src: for row in range(0, src.height, 512): for col in range(0, src.width, 512): window Window(col, row, 512, 512) data src.read(1, windowwindow) # 在这里对分块数据做滤波或其他处理如果你的场景只是预览效果也可以先建立金字塔只看概览层或者先重采样到较低分辨率做试算确认参数无误后再跑全分辨率省下的时间非常可观。5.5 重分类阈值不能拍脑袋做土地适宜性评价时我见过有人把坡度“0到2度”设为最适宜、“2到6度”设为适宜问为什么这么切回答是“感觉差不多”。这种拍脑袋的阈值会直接影响后续加权叠加的结果。正确做法是先看坡度的直方图或分位统计了解数据实际分布再结合行业规范或研究文献来定断点。分类的本质是切分决策空间断点应该反映数据本身的结构或业务规则而不是随便凑一个整数。6. 一次完整复盘从DEM到汇水区划定的全流程复习日的最后我用一个真实做过的小流域案例把整个流程串了一遍。虽然案例本身不复杂但每一步都能踩到前面提到的知识点对系统性理解很有帮助。6.1 案例背景与数据检查任务是要在某个丘陵地区划定一片面积约15平方公里的小流域边界。手头数据是5米分辨率的机载LiDAR DEM覆盖范围约10公里乘8公里投影是UTM 50N。拿到数据的第一件事不是直接跑工具而是先检查投影、分辨率和NoData情况顺便看一眼高程直方图确认是否存在明显的高程异常值。6.2 统一投影与重采样这个DEM已经是投影坐标系所以省掉了投影转换。但为了后续和矢量边界叠加不产生错位我把所有参与分析的数据都重采样到同一个5米网格并设置了基准栅格确保所有中间结果严格对齐。这一步做完后面的工具链跑起来才会干净。6.3 水文分析参数和检查点流程是标准的四步Fill填洼先自动填洼再对比填挖前后的高程变化检查有没有把喀斯特坑洞或人工水库周围的地形过度“抹平”。Flow Direction流向用D8算法输出流向栅格。Flow Accumulation流量累积输出每个像元汇入的累积像元数这个栅格在视觉上已经能看出水网的骨架。河网提取经过多次试阈值把阈值定为600对比卫星影像上的实际沟道后确认河网密度合理。随后把栅格河网转为矢量线再以流域出口点为汇水点生成集水区。每一步都有检查点流向栅格要确保没有大片无值区流量累积图要能隐约看到真实冲沟的形态河网矢量和影像上的季节沟道基本重合。如果中间某一步不对我的原则是倒回去看前一步而不是继续往后跑。6.4 结果验证用面积估算倒推合理性生成的集水区面积为15.2平方公里。为了判断这个结果是否合理我用区域统计算了流域内的平均年降水量和多年平均径流系数再简单估算一下理论年径流量如果该流域多年平均径流深度约600毫米则理论年径流量约为15.2平方公里×600毫米折算后大约是912万立方米。把这个估算值和附近水文站的实测数据对比如果落在同一量级就说明集水区边界基本可靠如果差很多就要检查是不是某个汇水点选错了位置。这个“面积倒推水量”的验证方法并不需要多复杂的模型却能在几分钟内给你一个强有力的参考信号比对着屏幕反复看边界要靠谱得多。复习日过完我最深的感受是栅格分析看起来工具繁多真正决定成败的往往不在工具本身而在数据准备阶段的那些枯燥细节——投影对不对、分辨率齐不齐、NoData处理干净没有。还有一点很实用的小习惯把常用栅格分析做成一页纸检查清单每次出图前按顺序过一遍比临时翻文档高效多了。