今天把栅格数据的空间分析整个过了一遍。说实话栅格分析和矢量分析虽然都叫空间分析但思维方式完全不一样。矢量玩的是“对象”一个点一条线一块面属性挂在对象身上栅格玩的是“场”整个研究区被切成一堆排列整齐的小方块所有信息都写在每个像元里。复习完最大的感触是很多人的栅格分析做不好根本不是工具不熟而是没把“像元”这件事想透。这篇文章算是我自己的复习笔记也是给正在学GIS、准备考试或者刚接触空间分析的同行看的。我不打算按教科书那样从定义抄起而是把我实际复习中用到的理解方式、踩过的坑、反复验证过的参数设置习惯都梳理一遍。内容偏实操适合那种“原理知道一点、但一操作就出问题”的人。计算机基础好一些的读者可以跟着推演计算逻辑基础薄一点的直接照步骤做也能跑通。1. 复习栅格分析得先重新认识“像元”1.1 栅格本质是一张按规则网格切分的大表栅格数据听起来很抽象但拆开了其实就是一层网格。你可以把研究区想象成一张棋盘棋盘上的每个格子就是一个像元格子的边长就是分辨率也叫像元大小。比如一副30米分辨率的DEM意思是地面每30米乘30米的区域被一个像元代表这个像元里只存一个数值。我第一次带学生做栅格分析时不少人会问一个问题既然一个像元只能存一个值那这30米见方的地方如果既有山地又有平地怎么办答案是看这个像元的取值规则。如果栅格采集时遵循“中心点取值”那就看像元中心落在什么地物上如果遵循“面积占优”那就看哪种地物占这块面积比例最大。这个细节平时不太起眼但在做土地利用分类、植被覆盖度估算时影响很大。同一块区域用两种规则生成的栅格面积统计结果能差出好几个百分点。栅格数据的第二个特性是波段。航拍影像或者卫星影像一般是多波段的比如红、绿、蓝三个波段合成一张真彩色影像再加近红外波段就能参与植被指数计算。而DEM、坡度、人口密度这些派生产品通常是单波段的每个像元只存一个数值。复习时我习惯把波段理解为“同一套网格下存的多张属性表”波段之间靠行列号严格对应不会错位。1.2 数据类型和NoData是栅格分析的两条隐形底线很多栅格操作结果不对问题不在操作本身而在数据类型的处理上。栅格数据按存储方式分为整型整数和浮点型小数。重分类、区域统计这类操作要求输入栅格是整型NDVI计算、坡度计算这种涉及运算的操作通常生成浮点型。ArcGIS里有一步“转为整型”的工具QGIS的栅格计算器也能用int()函数转换但转换时要考虑精度损失和信息取舍。更隐形的是NoData。栅格的空白区域并不是0而是NoData意思是“这个地方没有有效值”。区分0和NoData太重要了。计算NDVI时水体区域的像元可能被掩膜成NoData如果你把它当成0去参与后续统计平均值会被拉低一大截。我自己的习惯是在分析开始前做一次NoData审查。QGIS里可以用“属性”面板查看栅格数据集的统计信息和无效值设置ArcGIS里右键图层属性在源选项卡里能看到像元大小、行列数和NoData值。如果发现影像边缘一圈是NoData尽量用边界裁剪掉不然邻域统计、距离分析都会在边缘产生奇怪的误差。2. 栅格空间分析的四种基本模式一次讲透2.1 局部分析一个像元跟一个像元算局部分析是栅格分析里最基础也最常用的一类学术上叫Local Operation或Per-Cell Operation。它只看每个像元自己的值不关心周围像元有点类似Excel里对两列同行数据做公式运算。最典型的局部分析是NDVI计算公式是(NIR - Red) / (NIR Red)对应到栅格计算器就是两个波段逐像元相减再相除。还有重分类按数值区间把原来的像元值映射成新值比如把高程分成平原、丘陵、山地三个等级。这个操作本质上也是逐像元处理。局部分析最大优势是逻辑简单、运算效率高因为每个像元之间没有依赖关系可以并行计算。但正因为只看单像元它对“碎斑”问题无能为力。比如用阈值法从温度栅格里提取高温区会出现大量零星的单像元高温点这就是典型椒盐噪声。解决椒盐噪声靠局部分析自己不行得依赖接下来要说的邻域分析做平滑。2.2 邻域分析计算单位和邻居绑定在一起邻域分析在GIS里常被称为焦点统计Focal Statistics或滤波。它的核心是对每个像元开一个滑动窗口计算窗口内所有像元的值再把这个统计结果赋给中心像元。窗口常见的有3×3、5×5也可以按半径做圆形邻域还可以选环形、楔形。我用一个具体例子解释它的意义。假设有一副30米DEM你要提取山脊线。山脊线本质上是一条条相对两侧都高的线性体。如果单看一个像元你根本判断不出它是不是山脊因为你不知道旁边像元是比它高还是比它低。只有把周围3×3甚至5×5范围内的像元都拉进来比较才能判定“这个像元是局部高值”。这就是邻域分析的价值——它让栅格分析从“看独点”升级到“看局部结构”。坡度坡向计算是邻域分析的另一个高频场景。表面分析工具计算坡度时就是用一个3×3窗口去拟合局部地表平面进而求出水平面与这个拟合平面的夹角。ArcGIS里的“焦点统计”工具也可以自选统计类型比如算窗口内的平均值、最大值、标准差。做平滑滤波时我一般先试3×3均值窗口如果噪声还明显再上5×5但不能无脑加大窗口否则会丢失真实的地形细节。2.3 区域分析用矢量分区去统计栅格信息区域分析和邻域分析有点像但区别在于“分组依据”。邻域分析的滑动窗口形状固定是靠空间距离决定的区域分析的统计范围和位置完全由另一个图层决定这个图层可以是矢量面也可以是另一个栅格分区。典型操作是分区统计Zonal Statistics。举个例子你手头有一个省级行政区划面文件又有一副全国年降水量栅格你想知道每个省的平均降水量那就用分区统计。输入分区矢量选择降水量栅格统计类型选“平均值”输出的表格就是每个省的平均降水量。这个过程看起来简单但背后有个细节矢量分区在参与分区统计前必须先和栅格数据在空间上对齐。如果区划边界划分到了栅格像元的中间系统会用“像元中心点落在哪个区里”来确定归属边界位置的像元归属可能会有几个像元的偏差。另外还有个工具叫区域直方图它会统计每个分区内不同栅格值的像元数量分布。做土地利用转移矩阵时这个功能非常有用用两期土地利用栅格做分区统计第一期的每个类型作为分区看第二期中各类型像元数量就能得到转移矩阵。2.4 全局分析距离和连通性绕不开它全局分析不考虑“单个像元”而是把整个栅格视为一个连通体。最常见的全局分析是距离分析。简单欧氏距离很容易理解每个像元到最近一个源像元的直线距离。ArcGIS里“欧氏距离”工具输出结果同时会给两个栅格一个是距离栅格一个是方向栅格。这里有个隐藏细节欧氏距离用的是像元中心点之间的距离计算过程中对角方向的距离是乘根号2的这就是为什么图中看起来相对整齐的格状距离线会在对角方向产生一圈一圈的圆弧伪影。比欧氏距离更贴近实际的是成本距离Cost Distance。欧氏距离只看物理远近成本距离还考虑“走这段路要付出多少代价”。比如你从营地出发要穿过一片沼泽沼泽区域的成本栅格值很高即使绕远路过平坦草地总成本反而更低。成本距离计算基于像元间的移动代价做路径规划、保护区廊道分析时经常用到。它计算量比欧氏距离大得多因为它要迭代求解ArcGIS底层用的算法类似于最短路径的栅格版实现。水文分析也是全局分析的经典代表。从DEM提取流向D8算法、计算汇流累积量、提取河网、划分集水区每一步都涉及“水往哪流”的全局逻辑。填洼通常是第一步因为DEM里常有一些虚假的凹陷点如果不填掉水流会被困在这些坑里流向就会断掉。填洼本身是个有争议的操作真实地理环境里确实有天然洼地、溶蚀坑等模拟洪水淹没场景时填洼可能不合理。复习时要把这个“为什么填、什么时候不该填”拎出来想清楚面试和实操都有用处。3. 高频实操拆解从数据准备到地形因子计算3.1 数据预处理的第一步永远是统一投影和分辨率我见过太多人跳过这一步直接把数据喂给工具结果用两个不同投影的栅格做代数运算输出的结果完全是乱的。栅格空间分析对数据空间参考的要求比矢量更严格。矢量分析允许动态投影软件实时对齐就行栅格是像素级运算一旦投影不一致像元中心对不齐算出来的逐像元结果就没有物理意义。栅格分析的预处理我总结三步走把所有参与分析的栅格统一到同一个投影坐标系建议用投影坐标系不用地理坐标系因为地理坐标系的度单位对距离量算不友好。检查像元大小是否一致。如果不一致以最终成果需要的分辨率作为基准用重采样工具统一。设置环境里的“捕捉栅格”让所有输出栅格的像元边界和基准栅格严格对齐。第3步最容易忽略但特别关键。ArcGIS里在【环境设置】-【栅格分析】-【捕捉栅格】中指定基准栅格QGIS的处理框架里则是在“运行”前找到“重投影图层”或“对齐栅格”来处理。如果不设置捕捉栅格两个像元大小相同的栅格也可能因为原点坐标差半个像元而产生错位做差值运算时结果会出现一圈圈的条纹噪声。重采样方式也有讲究。离散数据土地利用类型用“最邻近法”保留原始像元值连续数据气温、降水用“双线性插值”或“三次卷积插值”能平滑过渡但不能用于离散分类数据否则会生成不存在的中间类型值。3.2 地形因子计算的原理与实操参数地形因子是栅格空间分析里最容易出成就感的内容因为坡度、坡向、曲率这些成果很直观看一眼就明白分析是干什么的。但它们的计算原理值得多花点时间我复习时把坡度、坡向和曲率放在一起对比理解效率很高。坡度是地面在某点的倾斜程度ArcGIS里用Horn算法计算。它把每个像元及其周围8个像元组成3×3窗口分别计算东西方向和南北方向的高程变化率dz/dx、dz/dy再用反正切函数得出坡度角。这个过程对分辨率很敏感在30米DEM上算出的坡度和10米DEM上算出的坡度数值分布明显不同。因为地形细节在低分辨率下被平滑了陡坡被“抹平”。所以做坡度分级时一定要先确认DEM分辨率和成果尺度匹配不能拿30米DEM去算精细到小地块的坡度。坡向表示地形面对的方向通常分成平地、北、东北、东、东南、南、西南、西、西北等类别。计算时也是从3×3窗口推导出的dz/dx和dz/dy入手用反正切换算成方位角。平地像元的坡向一般定义为-1这在制图符号化时要注意不要让-1被拉伸渲染成奇怪的颜色。曲率则反映地表地形的凹凸程度正值代表凸形坡面负值代表凹形坡面。做土壤侵蚀分析时曲率与坡度坡向一起使用效果更好。实际操作中我建议用10米DEM做坡度分析时窗口不必调大直接用默认3×3。但若是高分辨率DEM1米级别做城市地表分析反而可以考虑先做一次低通滤波因为1米数据里的建筑细节、噪声点太多直接算坡度会得到大量无意义的破碎陡坡。3.3 栅格计算器看懂表达式比点工具更重要局部分析最核心的工具就是栅格计算器。ArcGIS叫“地图代数”QGIS叫“栅格计算器”GDAL里用gdal_calc.py命令本质都一样用代数表达式把多个栅格输入变成新的栅格输出。很多人觉得栅格计算器就是个“加减乘除计算器”其实它还有条件运算、逻辑运算和数学函数。最常用的条件表达式是Con函数。ArcGIS语法是Con(condition_raster, true_value, false_value)QGIS里更接近Python布尔索引写法。我举个例子如果用NDVI提取植被覆盖区阈值取0.3那么大于0.3的地方赋值1否则赋值0。表达式可以写成Con(ndvi 0.3, 1, 0)。这个操作在QGIS栅格计算器里要写成ndvi1 0.3。条件成立的栅格输出1不成立的输出0。注意不同软件的条件表达式写法不一样千万别把ArcGIS的语法直接搬到QGIS里用。栅格计算器还有个容易踩的坑浮点型运算的NoData传播。如果输入栅格中有一个像元是NoData默认情况下输出栅格的这个像元也是NoData。做NDVI时只要红外影像里有个别像元因为传感器故障是NoData计算结果就凭空多出一个空洞。解决方案是在表达式里对NoData做专门处理或者先用焦点的邻域统计把NoData填补掉。再有就是“地图代数”中的名称引用。ArcGIS中调试栅格计算器时图层名称带空格或横线会导致表达式语法错误最好把图层重命名成不带空格的简洁名称。QGIS里使用“波段号”的写法相对灵活但也要求输入图层是在画布中加载的栅格图层不是文件路径。4. 工具选型与参数细节不同场景下的取舍经验4.1 三套主流程怎么选栅格分析主流工具无外乎三套ArcGIS带Spatial Analyst扩展、QGIS带SAGA/GRASS插件、GDAL命令行/脚本。我三套都用过感受是各有不可替代的地方。ArcGIS的好处是生态完整、工具链闭环。从工程管理、符号化出图到地理数据库存储一系列流程都很顺尤其Spatial Analyst扩展的“栅格分析环境设置”对新手特别友好。代价是商业授权贵而且批处理操作不够灵活遇到成百上千个栅格文件手工点工具能点到手酸。QGIS免费开源SAGA和GRASS这两个插件补充了大量专业工具。比如SAGA的多种地形分析工具、水文分析工具比ArcGIS原生工具更细GRASS的r.terraflow做水文填洼对大DEM的性能比ArcGIS更好。但QGIS的界面参数命名和ArcGIS不一样像“像元大小”有时候叫Cell Size有时候叫Grid Spacing刚切换工具时容易发懵。GDAL是命令行和Python的核心适合做批量、自动化、服务器端处理。gdal_calc.py处理栅格代数表达式、gdalwarp做重投影与重采样、gdal_polygonize做栅格转矢量这三条命令我几乎每次批量任务都用。缺点是没有可视化界面地形因子计算需要进一步调用其他库或编译工具。我的选型建议是教学和复习阶段优先ArcGIS或QGIS因为你能直接看到每一步中间结果方便理解原理工程批量任务用GDAL脚本如果做水文分析和地形分析还有更冷门的需求可以单独试试SAGA。工具之间没有绝对优劣能解决当前问题的就是合适的。4.2 三个我反复栽过跟头的参数第一个是“像元大小”设置。做分区统计时输出表里有一项是“面积”或“像元计数”它由像元大小决定。如果你在分区统计时没有注意输入栅格的像元大小统计出来的面积可能偏大或偏小。尤其两个分辨率相差很大的栅格叠在一起时用3×3窗口和9×9窗口得到的邻域水文参数差别非常大。第二个是“输出范围”设置。有些工具默认输出范围等于输入栅格范围而有些则是两个输入的交集。比如用两期栅格做差值如果它们的范围不完全一致输出会变成两者重叠区域。这些重叠区域之外的NoData常常让统计结果出现整片空白。我踩过一次用两期NDVI做变化检测结果输出栅格边缘缺了一块后来查发现是其中一期的影像剪切范围短了一截默认输出用了交集范围。所以做范围不同的栅格运算前一定要在环境设置里手动指定研究的完整范围。第三个是“像元统计类型”。做重分类或“以表格显示分区几何”时软件会提示你设置栅格统计类型整型还是浮点型。如果原始数据是浮点型但你想做分类计数需要先转成整型或重新分类。浮点型数据做直方图统计时像元值是连续分布的分类很难得到清晰边界转成整型后每一类会被归并到整数值统计才有意义。我把常见的栅格数据类型、可参与的命令类型整理成一个速查表方便复习时照着检查数据类型典型来源可使用的主要操作注意事项整型栅格土地利用分类、重分类结果区域统计、频率统计、转矢量转矢量前建议先“栅格清理”去除碎斑浮点型栅格高程DEM、气温插值、NDVI坡度坡向、栅格计算器、表面分析做“转为整型”时注意信息损失多波段栅格遥感影像、航拍DOM波段合成、指数计算、主成分分析参与代数运算前确认波段号和顺序二值掩膜栅格提取结果(0/1)掩膜提取、条件函数、距离分析0和NoData不要混为一谈5. 常见问题与排查技巧实录5.1 输出结果一大片NoData怎么办这应该是最常见的栅格分析问题。排查顺序我建议先确认原始范围打开属性里的“范围”和“像素宽度”看看是不是两个输入栅格的范围没有重叠。再看NoData设置如果输入栅格的NoData被设成0或带了背景值且没有正确识别分析时会把背景值当成有效值参与运算。还要注意栅格计算器里的空值传播规则。ArcGIS默认输出遇到NoData就输出NoDataQGIS里则可以通过“如果条件不成立时”选项把NoData替换成一个指定值。排查时可以先计算“IsNull”类型的判断把NoData区域可视化出来确认NoData分布是否符合预期。如果是单幅影像本身内部的NoData斑点造成的结果空洞可以尝试邻域填补在QGIS里用“焦点统计”将NoData视为缺失值用邻域有效像元的均值填补ArcGIS里等价操作是“焦点统计”后用Con函数把空值填上。这个方法对DEM孔洞比较有效但要注意填补值是光滑估计不代表真实地物。5.2 分辨率不一致引起的“对不齐”问题两个栅格的分辨率不同时不能直接做逐像元运算。比如一个30米分辨率的坡度和一个10米分辨率的NDVI它们虽然覆盖同一区域但像元网格错位直接相减属于“坏数据操作”。解决办法是先重采样让目标对齐到同一分辨率。这里我坚持一个原则重采样方向优先“粗转细”还是“细转粗”取决于你的应用场景。做水文学分析通常以地形数据的分辨率为主其他辅助数据采样到同一网格做植被变化检测则用统一的较粗分辨率降低噪声影响避免把多源传感器差异当成真实变化。检查是否对齐有一个实用技巧加载两个栅格把其中一个设为半透明放大到像元级别观察网格是否重合。如果发现网格偏移半个像元说明不是分辨率问题而是原点坐标不一致此时要调整环境设置里的“捕捉栅格”或使用“对齐栅格”工具。5.3 分类结果出现大量碎斑的清理技巧重分类或多值提取的结果经常会看到密密麻麻的细小斑块原因是我们只用单像元阈值判断完全忽略了空间连续性。解决思路有三类先用邻域统计做平滑做一次均值滤波或众数滤波把孤立碎斑“融化”进周边主导类型。使用“多数滤波”工具Majority Filter它专门针对分类栅格把零星像元替换成周围多数像元的值ArcGIS和SAGA里都有。做栅格转矢量后的碎面清理转到矢量后按面积阈值筛选掉过小的面片再融合相邻面转回栅格或直接用矢量成果。我在复习时特别强调一个点碎斑清理要控制尺度别为了图面美观把真实的地物细节也抹掉了。比如农业区的小块旱地本身就只有几百平方米如果阈值设得太大会把它们全删掉产出结果就会失去真实性。经验上阈值大小最好先用直方图或面积分布曲线看看碎斑的数量构成再决定砍到哪一档。5.4 栅格和矢量的来回转换转换坑得心里有数栅格分析经常需要和矢量数据配合。矢量转栅格时最关键的是字段选择转出去的是属性表中哪个字段转好之后是分值还是类别。ArcGIS里“面转栅格”工具里要选值字段和像元大小QGIS里是“栅格化”工具。如果选错字段后续所有统计都会变成错误的分类。栅格转矢量则要注意两个问题一是浮点型栅格直接转矢量会对每个栅格值产生独立多边形而浮点栅格里的噪声像元会产生海量碎面二是NoData区域转出来的多边形往往是不规则空洞如果只是做制图可以不考虑但做面积统计一定要先把NoData区域排除或填充。我遇到过最抓狂的一个问题是栅格转完矢量后的面积和原栅格的像元计数对不上。原因是矢量化边界是沿着像元边缘走的按理说面积应该等于像元数量乘以单像元面积。但如果栅格有NoData洞、或者在边缘有半像元对齐问题面积就会差一点。后来我总结了凡是做面积统计直接在栅格里用“区域统计”或“像元计数”计算不要在栅格转矢量后再统计面积避免二次转换引入误差。6. 复习建议和最后的经验这次完整复习下来我对栅格分析的理解比刚学时深了很多。最大的转变是我不再把栅格分析当成一个“按钮集合”我拿到数据时会先想这次分析是局部、邻域、区域还是全局模式是只处理单像元的值还是需要看像元之间的相对关系还是直接看全图的连通结构想清楚这一步工具选择就自然浮现出来了。我给准备系统学习或考试复习的朋友一个建议不要死记工具位置要按“数据结构→运算模式→工具实现”的链路来学习。先理解栅格像元和NoData这两个底层概念再掌握局部/邻域/区域/全局这四种运算模式最后落到ArcGIS和QGIS的实操里验证。每做一个练习都问自己一句“这个工具如果不用默认参数我改哪个参数会对结果产生什么影响”这个问题能帮你把实操经验真正变成原理理解。最后一个实用技巧。复习栅格空间分析时与其拿现成数据反复做我建议自己动手生成一份“人造栅格”用QGIS的栅格计算器生成一个像元值从1到100线性变化的栅格再生成几个圆形高值区。你自己知道真实情况长什么样做坡度、焦点统计、欧氏距离时就能直观判断工具算得对不对。这个方法看着笨但比刷十遍文档都管用。做过一次人造数据实验你对栅格分析的基本功就真正扎实了。