1. 栅格空间分析第一课先把“像素思维”扭转过来做GIS的人绕不开栅格数据。哪怕你今天只处理矢量数据只要你做坡度分析、做流域提取、做城市扩张监测、做地下水位模拟最终都会被拉回栅格的世界。原因很简单——很多自然现象是连续分布的而矢量结构本质上描述的是离散边界你没法用一个个多边形把一个光滑的水位面、温度场或者人口密度面管理好这时候必须靠规则格网。我在过去的项目里经常遇到这种场景甲方给了一个全国地下水位观测点shp要求做区域水位插值插完还要分区分级出图或者给你一个城市形态栅格数据集让你算不同城市建成区的破碎度、连片度指标。这些任务如果直接把矢量点和面拿来算你会发现操作链非常冗长而一旦转入栅格语言很多指标几行工具就能跑通。所以我打算把栅格数据空间分析做成一个系列从最常用、最容易踩坑的基础操作开始讲第一篇先说清楚三个问题栅格空间分析的底层逻辑是什么、拿到数据后怎么把底子打牢、以及最常用到的几个核心分析工具应该怎么用。这一篇适合谁适合刚入门的GISer也适合那些已经把ArcGIS或者QGIS工具栏翻过一遍但对栅格一直“知其一不知其二”的朋友。你放心我不会直接丢出一堆功能菜单让你无脑点我会把每个操作背后的原理讲清楚这样以后无论换什么软件、什么版本的Python库你都能很快迁移过去。2. 栅格数据的本质和空间分析思维2.1 一个栅格背后藏着哪些信息栅格数据说白了就是一块规则铺开的棋盘每个格子叫像元Cell或者像素Pixel每个像素里填一个数值这个数值可以代表高程、气温、降雨、地下水位、土地利用类型编码等任何空间属性。你可以把它理解成一个带地理坐标的Excel表格每一行每一列都有实际的位置信息而单元格里的数字就是我们要分析的对象。从这个角度看栅格有一个矢量无法替代的优势空间连续。比如地下水位今天你去井里测只能得到离散点但真实的水位面是一个连续曲面你必须借助栅格插值把它表达成一个个像元上的估计值。同一个逻辑也适用于城市形态分析——建筑物密度、透水面率、夜间灯光亮度这些指标本身就适合按格网统计因为城市边界不是一条干净利落的线而是慢慢过渡的场。栅格的“分辨率”是一个极其关键的概念它指每个像元在地面上代表的实际尺寸。比如一个30米分辨率的栅格意味着每个格子代表地面30乘30米的区域。分辨率越高数据量越大分析精度也越高但噪声和存储压力同时上升。做空间分析时我每次拿到栅格数据都会先看一眼像元大小因为它直接决定后续邻域分析、坡度计算时的尺度感。分辨率也决定了你“敢不敢”做某些分析。比如用30米DEM算坡度坡面逻辑是成立的如果你用一公里分辨率的DEM去算坡度那算出来的东西只能叫“宏观地势起伏趋势”不能叫工程意义上的坡度。很多新手拿到全国栅格就一通百通地分析最后结果图看起来挺漂亮实际已经失真了。这一点要刻在脑门上。2.2 栅格空间分析的三大思维模式做栅格分析前我建议你先建立三种思维模式。第一种是“字段思维”——像元上的数值是一个连续场的采样我们可以对它做各种数学运算和重分类这就像对地面上的物理量做算术。第二种是“邻域思维”——某一个像元的特征往往由周围像元决定空间自相关就在这儿起作用坡度、粗糙度、聚集度全是邻域计算的产物。第三种是“尺度思维”——在不同分辨率下做相同分析你会得出不同结论这不是结果错了而是空间尺度变了。这三种思维贯穿整个栅格空间分析过程。比如你在做城市形态栅格数据集分析时想算“透水面率”这个指标本质上就是一个邻域窗口内的均值计算想算“建成区连片度”就得靠空间统计判断一个像元是否被同类像元包围。只要把这三个思维模式建立起来后面所有工具都只是实现这些思想的按钮而已。另外我想提一个非常容易犯的错把栅格分析当成栅格“处理”。处理只是改数、裁剪、拼接、提取分析是在处理之上加入空间关系和统计逻辑。比如把一个NDVI栅格按阈值分成植被和非植被这是处理但如果你接着算植被覆盖像元在某一半径内的连通性或者计算林地斑块之间的隔离距离这才是分析。后面的实操里我会把这条界线说得很具体。2.3 工具选型建议免费的QGIS配Python最实用工具方面现在市面主流是ArcGIS和QGIS二选一再加上基于Python的rasterio、geopandas、numpy组合。我自己的习惯是日常空间分析用QGIS做交互式操作批量处理或需要复现的工作交给Python脚本ArcGIS偶尔用来对接老项目数据。之所以主推QGIS倒不是它比ArcGIS强到哪里去而是它免费开源许可证上没有坑放在个人项目里也完全够用。而且QGIS的栅格工具集非常完整从栅格计算器、重分类到坡度坡向一应俱全处理流程和ArcGIS高度相似你在QGIS里练熟了再切ArcGIS几乎没有学习成本。至于Python我是强烈建议从第一次碰栅格分析就开始接触的哪怕只是用脚本批量裁剪、拼接、投影转换都能帮你省下大量重复劳动。后续系列如果写深了我还会专门出一篇纯Python做栅格分析的配置和套路这次先留在篇尾简单带一嘴。说到底空间分析的成败不取决于软件牌子而是取决于数据质量和你的参数选择。接下来咱们进入实操前的关键环节——把数据底子打好。3. 准备底子从拿到数据到能动手分析3.1 你这套数据是不是“真”栅格先说一个我见过太多次的乌龙有些数据源所谓的“栅格”其实是一个esri grid文件夹或者一个带有.asc后缀的文本栅格甚至有些“shp栅格”实际是打了网格的矢量面。比如“我国地下水位栅格数据shp”这种名字一听就混了——shp是矢量格式栅格通常应该是tif、img、asc或者grid。如果不信这个邪你直接把shp扔进QGIS里看属性表保准发现一堆多边形而不是一个个像素。所以拿到数据第一步要做三件事确认格式、确认波段、确认像元类型。在QGIS里加载栅格后右键图层选“属性”在“信息”里能看到数据格式、尺寸、波段数、像元大小、像元深度。像元类型如果是整型Int还是浮点型Float很关键因为它决定你能否做除法、平均或者浮点阈值判断。地下水位这类连续变量正确的栅格通常应该是浮点型像素值有小数的可能性很高。如果你发现像元类型是整型要小心可能数据被四舍五入过了这种精度损失会直接影响插值后的统计口径。城市形态栅格数据集一般都是整型编码比如1代表建设用地2代表绿地这种分类栅格不适合做平滑和插值只适合做重分类和邻域统计。3.2 坐标系不统一一切免谈坐标系是栅格空间分析里最琐碎但最致命的问题。我见过有同事拿着WGS84经纬度坐标系的高程栅格和UTM投影坐标系的行政区面叠在一起做裁剪结果裁出来的图整个歪到海里。QGIS其实会在坐标系不一致时弹默认提示但很多人习惯直接点掉后面分析结果就莫名其妙。栅格数据的坐标系检查并不复杂在数据源加载后打开图层属性查看“CRS”字段。如果是WGS84通常显示“EPSG:4326”这是经纬度坐标单位是度如果是投影坐标则会显示类似“EPSG:32650”这样的编号单位是米。做面积、距离、坡度这类分析最好统一到投影坐标系因为经纬度在不同纬度下格网代表的实际面积差异很大直接拿它算距离会严重失真。具体操作上如果数据的坐标系不一致可以写一个小命令用GDAL统一转换。我最常用的命令行如下gdalwarp -t_srs EPSG:32650 -r bilinear input.tif output_utm.tif这里的-r参数是重采样方法如果你处理的是高程或连续变量用bilinear或cubic更顺滑如果你处理的是土地利用分类千万别用bilinear很容易插出“1.7类地物”这种荒诞结果老老实实用nearest最近邻。这种因为重采样方法没选对导致分类栅格被污染的问题我后面还会再提一次。3.3 裁剪、镶嵌与NoData预处理有了不同来源的数据后你会发现各自的范围边界不一致有的覆盖全图有的只覆盖局部裁边界线还可能有锯齿因为矢量边界和栅格像元的对齐方式不同。不管之后做坡度计算还是做叠加分析必须先把范围统一。QGIS里裁剪栅格最常用的是“按掩膜图层裁剪栅格”工具Mask图层选你要用的边界shp输入栅格选待裁剪文件注意在“输出无数据值”选项里填-9999或者直接读取源文件的NoData值。这里就涉及到一个NoData的概念——NoData代表这个像元没有有效观测值它不等于零所以分析时必须区分对待。我遇到过这样的问题地下水位栅格在山区有大量NoData直接做邻域平均或者栅格计算器运算时NoData会像瘟疫一样扩散导致结果图秃一块。处理办法是有两种先判断NoData区域面积是否过大如果不大可以用邻域填充或者中值平滑补上如果过大就尽量保持原始NoData范围在做最终出图时把NoData区域透明化并且在图例上标注清楚。预处理之后还要注意数据的空间对齐。做栅格叠加分析时如果两个栅格的像元大小不同、原点不同QGIS会自动重采样对齐但默认的重采样方法和对齐策略未必合理。我的习惯是在关键分析前手动把所有参与计算的栅格重采样到同一个分辨率、同一个原点。QGIS有个工具叫“对齐栅格”在工具箱搜索“align”就能找到勾选统一的CRS和像元大小后执行能有效减少后续“神秘偏移”问题。4. 核心实操五个高频场景的栅格分析全流程4.1 从DEM到坡度坡向一套完整的邻域运算逻辑先拿坡度坡向开刀因为它是栅格邻域分析的经典入门案例也是后面做水文分析、太阳能辐射评估、土壤侵蚀分析的基础。坡度本质上是地形表面在某一点的最大斜率用数学语言说就是该点高程对水平距离的方向导数而在离散栅格里这个导数是用3乘3的窗口边上的像元差分来拟合出来的。如果你用QGIS在工具箱里找到“坡度”输入DEM栅格输出度数或百分比即可。这里的“Z因子”参数很多人忽略它的作用是调整水平单位和垂直单位之间的比例如果你的DEM是经纬度坐标单位是度而高程单位是米那么不设Z因子会导致坡度全部失真。换算方式不复杂大致用111320乘以该纬度余弦的倒数但更稳妥的建议是做坡度分析前先把DEM投影到以米为单位的投影坐标系这样Z因子直接设为1就完了。坡向分析同样是3乘3邻域的事。QGIS的“坡向”工具会输出0到360度的方位角并同时定义平坦区域为-1。拿到坡向数据后可以继续把它重分类成九个方向分类图比如北坡、东北坡、东坡等这在后续分析植被分布或者山区建筑选址时非常常见。实操过程中有个小细节值得注意坡度和坡向计算对DEM的噪声很敏感如果原始DEM看起来有“条带感”或“台阶感”结果图会很花。我一般会在坡度分析前对DEM做一个轻度的平滑滤波比如QGIS里的“焦点统计”取3乘3窗口的中值连续做一次就行不要做太多否则会把真实地形磨平。4.2 栅格计算器把多张栅格做成一张科学结论栅格计算器是整个栅格空间分析里最实用的工具没有之一。它的逻辑和Excel公式很像你可以拿两张或更多栅格做加减乘除、布尔判断、条件运算。比如你要做一个“高风险滑坡区”识别公式可以是“坡度大于30度且植被覆盖率低于0.2的区域”在QGIS的栅格计算器里就能一行公式跑出来(slope1 30) AND (ndvi1 0.2)这段公式返回的结果是布尔栅格1代表条件成立0代表不成立。你随时可以把它乘以其他数值生成“风险等级”连续值。写公式时需要注意的是QGIS里1代表第一个波段如果你数据是多波段影像比如RGB三波段就得明确到底用哪个波段的数值运算。栅格计算器还有一个高频场景是用来做多因子叠加评分比如“城市扩张适宜性分析”你把地形坡度、交通距离、生态红线限制、人口密度等不同因子分别赋予权重最后用一条长公式加总。这里的关键问题是量纲统一坡度是0到90人口密度可能是0到几千直接相加没有意义。正确做法是先对各因子做0到1的归一化我再补充一句归一化方式最好选择“最大值最小值拉伸”并且注意避开异常高值否则结果会被几个极端点绑架。做个简单示范假设你有“距道路距离”栅格dist_road想转成“交通便利度”得分分值越高越便利。公式可以写成1 - (dist_road1 - 0) / (10000 - 0)这里的10000是人为设定的最大影响距离超过这个值就不再有区分度。用栅格计算器配合归一化你就能灵活组合各种数据做出一张真正属于自己的“叠加分析成果图”。4.3 重分类给连续栅格分等级的艺术重分类工具就是按阈值把连续栅格离散化。它的用途太广泛了DEM按海拔分成平原、丘陵、山地地下水位栅格按埋深分成浅、中、深夜间灯光影像按亮度分成低、中、高城市活动区。QGIS的重分类工具在“栅格”菜单下叫“重分类栅格”它会弹出一张表让你填重新分类的阈值区间和对应新值。比如对于城市形态栅格原编码假设1是建设用地、2是绿地、3是水域你想变成“建设/非建设”的二值图那么区间设置就是1到1映射为12到2映射为03到3映射为0NoData保持NoData即可。这里我特别想提醒新手一点重分类的归类方式必须要能解释、可追溯。你定义一个比如地下水位风险阈值埋深小于5米为高风险5到15米为中风险大于15米为低风险这个阈值不是随便拍的它要有水文地质依据要么引文献要么参照当地规范。在输出成果时最好把阈值依据写进说明文档否则评审专家一问“你为什么把10米作为分界”回答不上来就很尴尬。另外一个技巧是在做重分类之前先把栅格的直方图拉开看一下在QGIS图层面板右键“属性”在“直方图”选项卡里能看到像素值分布。比如地下水位值的分布如果集中在10到20米区间你就知道分级区间要在这个区间里多划分几级别把第一级设为0到50米结果一张图上九成都是同一个颜色等于什么也没分出来。4.4 距离分析算“影响范围”最强工具空间距离分析在栅格空间分析里占据非常特殊的地位因为很多现实问题都跟“距某个东西有多远”有关居民点到避难场所的距离、城市中心到边缘的扩张梯度、污染源周围的影响半径都属于这个范畴。矢量空间里做距离分析通常很笨重因为要算大量几何距离到了栅格里每个像元到最近目标像元的距离能一次性批量算完。QGIS里对应的工具叫“欧氏距离”在“栅格分析”菜单下。输入一个源栅格比如坑塘、商业中心、地铁站输出得到一张新的栅格每个像元代表该位置到最近“源”的距离。配合上重分类你就可以快速画出“距市中心20公里/40公里/80公里”的同心圆辐射范围图这在城市形态数据集分析里几乎是必用操作。需要注意的参数是“距离单位”因为输入输出栅格必须是有投影坐标系的QGIS默认按照投影坐标系的单位输出距离所以用经纬度坐标系做欧氏距离的结果是以度为单位的基本等于废纸一张。因此我又要强调一次预处理的必要性——距离分析前必须投影转换到米制坐标系。更进阶的中心位置分析是把距离栅格和另外的栅格一起叠加比如“距水源距离”乘上“坡度因子”模拟取水便利度。你可以在栅格计算器里组合也可以是“成本距离分析”工具这个工具会考虑地表阻力差异模拟出的不是轴对称的欧氏圆而是真实受力条件下的等时圈、等成本圈。这里先埋个伏笔后续写到水文和可达性分析时再展开。4.5 邻域统计滑动窗口里的省域城市形态邻域统计是栅格分析里真正的“高手操作”因为空间自相关性、平滑、聚集度、斑块密度全靠它。它的核心思想是在每一个像元周围开一个固定大小的窗口比如3乘3、5乘5或者90米乘90米然后计算窗口内所有像元的某种统计量比如平均值、最大值、标准差、众数把这个统计量赋给中心像元。QGIS里的工具叫“焦点统计”Focal Statistics。举个例子你想看某省建设用地密度分布原始栅格只有0和1两个值0代表非建设用地1代表建设用地。直接用原始栅格绘制等密度图图面会非常破碎根本看不出空间聚集趋势。如果你用焦点统计统计类型选“平均值”窗口选一个比如15乘15像元新栅格的值就变成每个中心位置周边225个像元内建设用地的占比也就是“局部建成区比例”。这个指标连续、平滑出图后一眼就能看出哪个地区城市化连片度高哪个地区呈星星点点孤立分布。另外“标准差”型焦点统计常用于识别空间异质性比如城市地表温度的空间变异热点温度标准差高的区域意味着热环境剧烈跳跃往往是城市热岛边界或水体边缘。窗口大小没有绝对标准但必须交代清楚我通常根据数据的空间分辨率和分析尺度来定窗口半径像30米分辨率的建设用地数据想看宏观城市形态用15到25个像元窗口如果只想看周边社区尺度3到5个像元窗口就足够。窗口越大平滑力度越强但也会掩盖局地细节。实操中我会做两到三组窗口大小对比图选看起来既能抓到主要城市群形态又不至于碎成渣的窗口尺寸这个判断准则很难量化但做得多了自然有手感。5. 分析中常见的坑和排查套路5.1 NoData引发的“鬼影”区域做栅格计算器或者焦点统计时最常遇到的怪现象是结果图上有大片空白的洞而且这些洞不在数据源头存在或者越扩越大。这基本就是NoData扩散导致的。因为一旦参与计算的某个像元是NoData默认规则会把结果像元也变成NoData窗口滑过去之后原本有数据的边界也会被逐渐“啃掉”形成所谓的NoData扩张效应。解决办法有几种。如果只是小范围缺失用“焦点统计”的“忽略NoData”选项计算局部统计时跳过空值填一个估计值如果缺失面积较大可以用多时相数据或者插值先补一次底图。QGIS里比较朴素的补洞方式是把NoData区域设置成一个极端的哨兵值比如-99999再用“栅格计算器”判断哨兵位置附近的有效值取平均但这种方法后处理复杂不如从源头检查一下原始数据有没有异常。我自己的排查顺序是先在原图上做一遍“NoData的可视化”把NoData单色显示看看它们是否连片、是否出现在关键分析区。如果NoData区域的分布本身没有空间规律比如不像水体或边界导致说明数据预处理有问题比如投影范围设置错误导致重采样时部分区域没有源像素落进来。5.2 坐标系不一致导致的“乾坤大挪移”坐标系问题在实操里几乎天天见。表现形态主要有两种第一是从不同部门拿到的数据坐标系代号不同比如一个用CGCS2000一个用WGS84由于两者椭球体参数非常接近很多工具显示不出明显区别但叠加后会有几十米甚至上百米的系统性偏移。对于公里级分析可能无所谓但做地下水位或建设用地这种城市尺度的精细分析偏移一旦出现后续像元对齐就全是错的。第二种情况是投影坐标系和地理坐标系混用这种更明显图面直接错位。遇到这种问题我建议在项目启动前就定一个“项目坐标系”全流程统一到它。通常城市尺度我会选CGCS2000 / Gauss-Kruger zone对应的中央经线全国尺度则用Albers等积投影。定好之后所有输入栅格第一件事就是检查CRS不一致就执行gdalwarp重投影不要心存侥幸。5.3 数据量大、内存不足导致运行迟缓栅格文件往往非常大全国尺度的城市形态栅格动辄几个GB如果在QGIS里直接跑焦点统计或者欧氏距离可能运行几分钟甚至直接崩溃。这时候有几个务实办法。第一先裁剪到需要研究的目标区域别全图硬跑第二降低分辨率比如从30米重采样到100米或者250米有些宏观形态分析在250米分辨率下结论依然稳健第三利用金字塔结构QGIS会在加载大栅格时自动生成金字塔可以明显加速显示但计算工具往往是逐像元处理的金字塔帮不上太多忙。我用rasterio做Python批处理时习惯分块读取和计算避免一次性把全图压进内存。比如用rasterio打开栅格后用window参数读入小块数据处理完再写回。这种分块处理对超大栅格几乎是最优解。但如果只是偶尔做一两次交互式分析QGIS的裁剪和重采样就够用。5.4 常见问题速查表现象可能原因排查步骤解决方案裁剪结果有大片黑色区域NoData被当成0值查看栅格“无数据值”重设NoData值并重裁坡度值全部异常偏大DEM坐标系为经纬度检查CRS投影到米制坐标再算重采样后分类栅格出现小数使用了连续重采样方法查看像元类型用最近邻法重采样道路距离分析中心偏移输入矢量未投影检查源的CRS统一投影后再转栅格焦点统计结果出白斑窗口内NoData未处理检查NoData区域开启忽略NoData选项两张栅格叠加后有错位原点或像素大小不一致对比栅格属性信息用对齐栅格工具统一网格6. Python脚本批量处理补充如果你要处理的数据多了比如上百个城市的城市形态栅格或者多年份的地下水位数据再手动在QGIS里一个个点工具是真的浪费时间。建议把高频操作用Python封装起来核心库是rasterio、numpy、geopandas再加上一些GDAL的命令行调用。举个例子批量裁剪所有栅格到某个shp边界范围import rasterio from rasterio.mask import mask import geopandas as gpd shape gpd.read_file(region.shp) with rasterio.open(input.tif) as src: out_image, out_transform mask(src, shape.geometry, cropTrue) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_crop.tif, w, **out_meta) as dst: dst.write(out_image)注意在mask前确保shape和栅格使用相同坐标系如果坐标系不一致先调用to_crs统一。再比如批量栅格计算器公式本质上就是把两个rasterio打开的文件读成numpy数组做数组运算后写新tif。numpy数组的广播逻辑和QGIS栅格计算器几乎一样只要注意NoData掩膜的处理就行import numpy as np slope rasterio.open(slope.tif).read(1) ndvi rasterio.open(ndvi.tif).read(1) nodata_slope -9999 slope np.where(slope nodata_slope, np.nan, slope) result np.where((slope 30) (ndvi 0.2), 1, 0)这种脚本的最大优势是可复现今天跑一遍和半年后跑一遍参数一致结果就一致不会因为手动误操作产生差异。我对所有做空间分析的团队都建议至少把预处理环节脚本化它带来的效率提升在后续系列中会越来越明显。7. 最后的经验之谈在我做项目的时候经常有人问为什么同一份数据两个人跑出两张完全不同的图。答案往往不是操作不同而是他们在分析前对数据的理解不一样——一个检查过坐标系、处理过NoData、选择过合适的重采样方法另一个直接打开原始数据飞速点击工具。栅格空间分析这个领域前期的数据整理工作大约占整个项目一半时间这不是浪费时间而是给结果上保险。所以我的建议是不要急着追求工具数量和技术炫技先把每个栅格像元代表什么物理意义搞清楚再把分辨率、坐标系、NoData、波段信息这四个基础属性吃透。地基打牢了后面什么坡度坡向、距离分析、邻域统计、成本分析真的只是多了一个又一个按钮而已。最后再分享一个小技巧吧。我每次拿到一份新的栅格数据都会先在QGIS里用单波段伪彩色方式看一遍原始值同时打开直方图拖拽色带拉伸范围。这一个动作看起来不起眼但它能帮你快速发现像元值的异常分布比如有特别离谱的高值或者大量集中在某个区间这些信息会直接影响你后面怎么设阈值、怎么归类。栅格分析大多数时候不是算法不够强而是你对输入的“手感”不够准。这份手感只能靠一次次打开直方图、一次次试不同参数慢慢磨出来。 栅格空间分析实战：从像素思维到核心工具一次讲清
1. 栅格空间分析第一课先把“像素思维”扭转过来做GIS的人绕不开栅格数据。哪怕你今天只处理矢量数据只要你做坡度分析、做流域提取、做城市扩张监测、做地下水位模拟最终都会被拉回栅格的世界。原因很简单——很多自然现象是连续分布的而矢量结构本质上描述的是离散边界你没法用一个个多边形把一个光滑的水位面、温度场或者人口密度面管理好这时候必须靠规则格网。我在过去的项目里经常遇到这种场景甲方给了一个全国地下水位观测点shp要求做区域水位插值插完还要分区分级出图或者给你一个城市形态栅格数据集让你算不同城市建成区的破碎度、连片度指标。这些任务如果直接把矢量点和面拿来算你会发现操作链非常冗长而一旦转入栅格语言很多指标几行工具就能跑通。所以我打算把栅格数据空间分析做成一个系列从最常用、最容易踩坑的基础操作开始讲第一篇先说清楚三个问题栅格空间分析的底层逻辑是什么、拿到数据后怎么把底子打牢、以及最常用到的几个核心分析工具应该怎么用。这一篇适合谁适合刚入门的GISer也适合那些已经把ArcGIS或者QGIS工具栏翻过一遍但对栅格一直“知其一不知其二”的朋友。你放心我不会直接丢出一堆功能菜单让你无脑点我会把每个操作背后的原理讲清楚这样以后无论换什么软件、什么版本的Python库你都能很快迁移过去。2. 栅格数据的本质和空间分析思维2.1 一个栅格背后藏着哪些信息栅格数据说白了就是一块规则铺开的棋盘每个格子叫像元Cell或者像素Pixel每个像素里填一个数值这个数值可以代表高程、气温、降雨、地下水位、土地利用类型编码等任何空间属性。你可以把它理解成一个带地理坐标的Excel表格每一行每一列都有实际的位置信息而单元格里的数字就是我们要分析的对象。从这个角度看栅格有一个矢量无法替代的优势空间连续。比如地下水位今天你去井里测只能得到离散点但真实的水位面是一个连续曲面你必须借助栅格插值把它表达成一个个像元上的估计值。同一个逻辑也适用于城市形态分析——建筑物密度、透水面率、夜间灯光亮度这些指标本身就适合按格网统计因为城市边界不是一条干净利落的线而是慢慢过渡的场。栅格的“分辨率”是一个极其关键的概念它指每个像元在地面上代表的实际尺寸。比如一个30米分辨率的栅格意味着每个格子代表地面30乘30米的区域。分辨率越高数据量越大分析精度也越高但噪声和存储压力同时上升。做空间分析时我每次拿到栅格数据都会先看一眼像元大小因为它直接决定后续邻域分析、坡度计算时的尺度感。分辨率也决定了你“敢不敢”做某些分析。比如用30米DEM算坡度坡面逻辑是成立的如果你用一公里分辨率的DEM去算坡度那算出来的东西只能叫“宏观地势起伏趋势”不能叫工程意义上的坡度。很多新手拿到全国栅格就一通百通地分析最后结果图看起来挺漂亮实际已经失真了。这一点要刻在脑门上。2.2 栅格空间分析的三大思维模式做栅格分析前我建议你先建立三种思维模式。第一种是“字段思维”——像元上的数值是一个连续场的采样我们可以对它做各种数学运算和重分类这就像对地面上的物理量做算术。第二种是“邻域思维”——某一个像元的特征往往由周围像元决定空间自相关就在这儿起作用坡度、粗糙度、聚集度全是邻域计算的产物。第三种是“尺度思维”——在不同分辨率下做相同分析你会得出不同结论这不是结果错了而是空间尺度变了。这三种思维贯穿整个栅格空间分析过程。比如你在做城市形态栅格数据集分析时想算“透水面率”这个指标本质上就是一个邻域窗口内的均值计算想算“建成区连片度”就得靠空间统计判断一个像元是否被同类像元包围。只要把这三个思维模式建立起来后面所有工具都只是实现这些思想的按钮而已。另外我想提一个非常容易犯的错把栅格分析当成栅格“处理”。处理只是改数、裁剪、拼接、提取分析是在处理之上加入空间关系和统计逻辑。比如把一个NDVI栅格按阈值分成植被和非植被这是处理但如果你接着算植被覆盖像元在某一半径内的连通性或者计算林地斑块之间的隔离距离这才是分析。后面的实操里我会把这条界线说得很具体。2.3 工具选型建议免费的QGIS配Python最实用工具方面现在市面主流是ArcGIS和QGIS二选一再加上基于Python的rasterio、geopandas、numpy组合。我自己的习惯是日常空间分析用QGIS做交互式操作批量处理或需要复现的工作交给Python脚本ArcGIS偶尔用来对接老项目数据。之所以主推QGIS倒不是它比ArcGIS强到哪里去而是它免费开源许可证上没有坑放在个人项目里也完全够用。而且QGIS的栅格工具集非常完整从栅格计算器、重分类到坡度坡向一应俱全处理流程和ArcGIS高度相似你在QGIS里练熟了再切ArcGIS几乎没有学习成本。至于Python我是强烈建议从第一次碰栅格分析就开始接触的哪怕只是用脚本批量裁剪、拼接、投影转换都能帮你省下大量重复劳动。后续系列如果写深了我还会专门出一篇纯Python做栅格分析的配置和套路这次先留在篇尾简单带一嘴。说到底空间分析的成败不取决于软件牌子而是取决于数据质量和你的参数选择。接下来咱们进入实操前的关键环节——把数据底子打好。3. 准备底子从拿到数据到能动手分析3.1 你这套数据是不是“真”栅格先说一个我见过太多次的乌龙有些数据源所谓的“栅格”其实是一个esri grid文件夹或者一个带有.asc后缀的文本栅格甚至有些“shp栅格”实际是打了网格的矢量面。比如“我国地下水位栅格数据shp”这种名字一听就混了——shp是矢量格式栅格通常应该是tif、img、asc或者grid。如果不信这个邪你直接把shp扔进QGIS里看属性表保准发现一堆多边形而不是一个个像素。所以拿到数据第一步要做三件事确认格式、确认波段、确认像元类型。在QGIS里加载栅格后右键图层选“属性”在“信息”里能看到数据格式、尺寸、波段数、像元大小、像元深度。像元类型如果是整型Int还是浮点型Float很关键因为它决定你能否做除法、平均或者浮点阈值判断。地下水位这类连续变量正确的栅格通常应该是浮点型像素值有小数的可能性很高。如果你发现像元类型是整型要小心可能数据被四舍五入过了这种精度损失会直接影响插值后的统计口径。城市形态栅格数据集一般都是整型编码比如1代表建设用地2代表绿地这种分类栅格不适合做平滑和插值只适合做重分类和邻域统计。3.2 坐标系不统一一切免谈坐标系是栅格空间分析里最琐碎但最致命的问题。我见过有同事拿着WGS84经纬度坐标系的高程栅格和UTM投影坐标系的行政区面叠在一起做裁剪结果裁出来的图整个歪到海里。QGIS其实会在坐标系不一致时弹默认提示但很多人习惯直接点掉后面分析结果就莫名其妙。栅格数据的坐标系检查并不复杂在数据源加载后打开图层属性查看“CRS”字段。如果是WGS84通常显示“EPSG:4326”这是经纬度坐标单位是度如果是投影坐标则会显示类似“EPSG:32650”这样的编号单位是米。做面积、距离、坡度这类分析最好统一到投影坐标系因为经纬度在不同纬度下格网代表的实际面积差异很大直接拿它算距离会严重失真。具体操作上如果数据的坐标系不一致可以写一个小命令用GDAL统一转换。我最常用的命令行如下gdalwarp -t_srs EPSG:32650 -r bilinear input.tif output_utm.tif这里的-r参数是重采样方法如果你处理的是高程或连续变量用bilinear或cubic更顺滑如果你处理的是土地利用分类千万别用bilinear很容易插出“1.7类地物”这种荒诞结果老老实实用nearest最近邻。这种因为重采样方法没选对导致分类栅格被污染的问题我后面还会再提一次。3.3 裁剪、镶嵌与NoData预处理有了不同来源的数据后你会发现各自的范围边界不一致有的覆盖全图有的只覆盖局部裁边界线还可能有锯齿因为矢量边界和栅格像元的对齐方式不同。不管之后做坡度计算还是做叠加分析必须先把范围统一。QGIS里裁剪栅格最常用的是“按掩膜图层裁剪栅格”工具Mask图层选你要用的边界shp输入栅格选待裁剪文件注意在“输出无数据值”选项里填-9999或者直接读取源文件的NoData值。这里就涉及到一个NoData的概念——NoData代表这个像元没有有效观测值它不等于零所以分析时必须区分对待。我遇到过这样的问题地下水位栅格在山区有大量NoData直接做邻域平均或者栅格计算器运算时NoData会像瘟疫一样扩散导致结果图秃一块。处理办法是有两种先判断NoData区域面积是否过大如果不大可以用邻域填充或者中值平滑补上如果过大就尽量保持原始NoData范围在做最终出图时把NoData区域透明化并且在图例上标注清楚。预处理之后还要注意数据的空间对齐。做栅格叠加分析时如果两个栅格的像元大小不同、原点不同QGIS会自动重采样对齐但默认的重采样方法和对齐策略未必合理。我的习惯是在关键分析前手动把所有参与计算的栅格重采样到同一个分辨率、同一个原点。QGIS有个工具叫“对齐栅格”在工具箱搜索“align”就能找到勾选统一的CRS和像元大小后执行能有效减少后续“神秘偏移”问题。4. 核心实操五个高频场景的栅格分析全流程4.1 从DEM到坡度坡向一套完整的邻域运算逻辑先拿坡度坡向开刀因为它是栅格邻域分析的经典入门案例也是后面做水文分析、太阳能辐射评估、土壤侵蚀分析的基础。坡度本质上是地形表面在某一点的最大斜率用数学语言说就是该点高程对水平距离的方向导数而在离散栅格里这个导数是用3乘3的窗口边上的像元差分来拟合出来的。如果你用QGIS在工具箱里找到“坡度”输入DEM栅格输出度数或百分比即可。这里的“Z因子”参数很多人忽略它的作用是调整水平单位和垂直单位之间的比例如果你的DEM是经纬度坐标单位是度而高程单位是米那么不设Z因子会导致坡度全部失真。换算方式不复杂大致用111320乘以该纬度余弦的倒数但更稳妥的建议是做坡度分析前先把DEM投影到以米为单位的投影坐标系这样Z因子直接设为1就完了。坡向分析同样是3乘3邻域的事。QGIS的“坡向”工具会输出0到360度的方位角并同时定义平坦区域为-1。拿到坡向数据后可以继续把它重分类成九个方向分类图比如北坡、东北坡、东坡等这在后续分析植被分布或者山区建筑选址时非常常见。实操过程中有个小细节值得注意坡度和坡向计算对DEM的噪声很敏感如果原始DEM看起来有“条带感”或“台阶感”结果图会很花。我一般会在坡度分析前对DEM做一个轻度的平滑滤波比如QGIS里的“焦点统计”取3乘3窗口的中值连续做一次就行不要做太多否则会把真实地形磨平。4.2 栅格计算器把多张栅格做成一张科学结论栅格计算器是整个栅格空间分析里最实用的工具没有之一。它的逻辑和Excel公式很像你可以拿两张或更多栅格做加减乘除、布尔判断、条件运算。比如你要做一个“高风险滑坡区”识别公式可以是“坡度大于30度且植被覆盖率低于0.2的区域”在QGIS的栅格计算器里就能一行公式跑出来(slope1 30) AND (ndvi1 0.2)这段公式返回的结果是布尔栅格1代表条件成立0代表不成立。你随时可以把它乘以其他数值生成“风险等级”连续值。写公式时需要注意的是QGIS里1代表第一个波段如果你数据是多波段影像比如RGB三波段就得明确到底用哪个波段的数值运算。栅格计算器还有一个高频场景是用来做多因子叠加评分比如“城市扩张适宜性分析”你把地形坡度、交通距离、生态红线限制、人口密度等不同因子分别赋予权重最后用一条长公式加总。这里的关键问题是量纲统一坡度是0到90人口密度可能是0到几千直接相加没有意义。正确做法是先对各因子做0到1的归一化我再补充一句归一化方式最好选择“最大值最小值拉伸”并且注意避开异常高值否则结果会被几个极端点绑架。做个简单示范假设你有“距道路距离”栅格dist_road想转成“交通便利度”得分分值越高越便利。公式可以写成1 - (dist_road1 - 0) / (10000 - 0)这里的10000是人为设定的最大影响距离超过这个值就不再有区分度。用栅格计算器配合归一化你就能灵活组合各种数据做出一张真正属于自己的“叠加分析成果图”。4.3 重分类给连续栅格分等级的艺术重分类工具就是按阈值把连续栅格离散化。它的用途太广泛了DEM按海拔分成平原、丘陵、山地地下水位栅格按埋深分成浅、中、深夜间灯光影像按亮度分成低、中、高城市活动区。QGIS的重分类工具在“栅格”菜单下叫“重分类栅格”它会弹出一张表让你填重新分类的阈值区间和对应新值。比如对于城市形态栅格原编码假设1是建设用地、2是绿地、3是水域你想变成“建设/非建设”的二值图那么区间设置就是1到1映射为12到2映射为03到3映射为0NoData保持NoData即可。这里我特别想提醒新手一点重分类的归类方式必须要能解释、可追溯。你定义一个比如地下水位风险阈值埋深小于5米为高风险5到15米为中风险大于15米为低风险这个阈值不是随便拍的它要有水文地质依据要么引文献要么参照当地规范。在输出成果时最好把阈值依据写进说明文档否则评审专家一问“你为什么把10米作为分界”回答不上来就很尴尬。另外一个技巧是在做重分类之前先把栅格的直方图拉开看一下在QGIS图层面板右键“属性”在“直方图”选项卡里能看到像素值分布。比如地下水位值的分布如果集中在10到20米区间你就知道分级区间要在这个区间里多划分几级别把第一级设为0到50米结果一张图上九成都是同一个颜色等于什么也没分出来。4.4 距离分析算“影响范围”最强工具空间距离分析在栅格空间分析里占据非常特殊的地位因为很多现实问题都跟“距某个东西有多远”有关居民点到避难场所的距离、城市中心到边缘的扩张梯度、污染源周围的影响半径都属于这个范畴。矢量空间里做距离分析通常很笨重因为要算大量几何距离到了栅格里每个像元到最近目标像元的距离能一次性批量算完。QGIS里对应的工具叫“欧氏距离”在“栅格分析”菜单下。输入一个源栅格比如坑塘、商业中心、地铁站输出得到一张新的栅格每个像元代表该位置到最近“源”的距离。配合上重分类你就可以快速画出“距市中心20公里/40公里/80公里”的同心圆辐射范围图这在城市形态数据集分析里几乎是必用操作。需要注意的参数是“距离单位”因为输入输出栅格必须是有投影坐标系的QGIS默认按照投影坐标系的单位输出距离所以用经纬度坐标系做欧氏距离的结果是以度为单位的基本等于废纸一张。因此我又要强调一次预处理的必要性——距离分析前必须投影转换到米制坐标系。更进阶的中心位置分析是把距离栅格和另外的栅格一起叠加比如“距水源距离”乘上“坡度因子”模拟取水便利度。你可以在栅格计算器里组合也可以是“成本距离分析”工具这个工具会考虑地表阻力差异模拟出的不是轴对称的欧氏圆而是真实受力条件下的等时圈、等成本圈。这里先埋个伏笔后续写到水文和可达性分析时再展开。4.5 邻域统计滑动窗口里的省域城市形态邻域统计是栅格分析里真正的“高手操作”因为空间自相关性、平滑、聚集度、斑块密度全靠它。它的核心思想是在每一个像元周围开一个固定大小的窗口比如3乘3、5乘5或者90米乘90米然后计算窗口内所有像元的某种统计量比如平均值、最大值、标准差、众数把这个统计量赋给中心像元。QGIS里的工具叫“焦点统计”Focal Statistics。举个例子你想看某省建设用地密度分布原始栅格只有0和1两个值0代表非建设用地1代表建设用地。直接用原始栅格绘制等密度图图面会非常破碎根本看不出空间聚集趋势。如果你用焦点统计统计类型选“平均值”窗口选一个比如15乘15像元新栅格的值就变成每个中心位置周边225个像元内建设用地的占比也就是“局部建成区比例”。这个指标连续、平滑出图后一眼就能看出哪个地区城市化连片度高哪个地区呈星星点点孤立分布。另外“标准差”型焦点统计常用于识别空间异质性比如城市地表温度的空间变异热点温度标准差高的区域意味着热环境剧烈跳跃往往是城市热岛边界或水体边缘。窗口大小没有绝对标准但必须交代清楚我通常根据数据的空间分辨率和分析尺度来定窗口半径像30米分辨率的建设用地数据想看宏观城市形态用15到25个像元窗口如果只想看周边社区尺度3到5个像元窗口就足够。窗口越大平滑力度越强但也会掩盖局地细节。实操中我会做两到三组窗口大小对比图选看起来既能抓到主要城市群形态又不至于碎成渣的窗口尺寸这个判断准则很难量化但做得多了自然有手感。5. 分析中常见的坑和排查套路5.1 NoData引发的“鬼影”区域做栅格计算器或者焦点统计时最常遇到的怪现象是结果图上有大片空白的洞而且这些洞不在数据源头存在或者越扩越大。这基本就是NoData扩散导致的。因为一旦参与计算的某个像元是NoData默认规则会把结果像元也变成NoData窗口滑过去之后原本有数据的边界也会被逐渐“啃掉”形成所谓的NoData扩张效应。解决办法有几种。如果只是小范围缺失用“焦点统计”的“忽略NoData”选项计算局部统计时跳过空值填一个估计值如果缺失面积较大可以用多时相数据或者插值先补一次底图。QGIS里比较朴素的补洞方式是把NoData区域设置成一个极端的哨兵值比如-99999再用“栅格计算器”判断哨兵位置附近的有效值取平均但这种方法后处理复杂不如从源头检查一下原始数据有没有异常。我自己的排查顺序是先在原图上做一遍“NoData的可视化”把NoData单色显示看看它们是否连片、是否出现在关键分析区。如果NoData区域的分布本身没有空间规律比如不像水体或边界导致说明数据预处理有问题比如投影范围设置错误导致重采样时部分区域没有源像素落进来。5.2 坐标系不一致导致的“乾坤大挪移”坐标系问题在实操里几乎天天见。表现形态主要有两种第一是从不同部门拿到的数据坐标系代号不同比如一个用CGCS2000一个用WGS84由于两者椭球体参数非常接近很多工具显示不出明显区别但叠加后会有几十米甚至上百米的系统性偏移。对于公里级分析可能无所谓但做地下水位或建设用地这种城市尺度的精细分析偏移一旦出现后续像元对齐就全是错的。第二种情况是投影坐标系和地理坐标系混用这种更明显图面直接错位。遇到这种问题我建议在项目启动前就定一个“项目坐标系”全流程统一到它。通常城市尺度我会选CGCS2000 / Gauss-Kruger zone对应的中央经线全国尺度则用Albers等积投影。定好之后所有输入栅格第一件事就是检查CRS不一致就执行gdalwarp重投影不要心存侥幸。5.3 数据量大、内存不足导致运行迟缓栅格文件往往非常大全国尺度的城市形态栅格动辄几个GB如果在QGIS里直接跑焦点统计或者欧氏距离可能运行几分钟甚至直接崩溃。这时候有几个务实办法。第一先裁剪到需要研究的目标区域别全图硬跑第二降低分辨率比如从30米重采样到100米或者250米有些宏观形态分析在250米分辨率下结论依然稳健第三利用金字塔结构QGIS会在加载大栅格时自动生成金字塔可以明显加速显示但计算工具往往是逐像元处理的金字塔帮不上太多忙。我用rasterio做Python批处理时习惯分块读取和计算避免一次性把全图压进内存。比如用rasterio打开栅格后用window参数读入小块数据处理完再写回。这种分块处理对超大栅格几乎是最优解。但如果只是偶尔做一两次交互式分析QGIS的裁剪和重采样就够用。5.4 常见问题速查表现象可能原因排查步骤解决方案裁剪结果有大片黑色区域NoData被当成0值查看栅格“无数据值”重设NoData值并重裁坡度值全部异常偏大DEM坐标系为经纬度检查CRS投影到米制坐标再算重采样后分类栅格出现小数使用了连续重采样方法查看像元类型用最近邻法重采样道路距离分析中心偏移输入矢量未投影检查源的CRS统一投影后再转栅格焦点统计结果出白斑窗口内NoData未处理检查NoData区域开启忽略NoData选项两张栅格叠加后有错位原点或像素大小不一致对比栅格属性信息用对齐栅格工具统一网格6. Python脚本批量处理补充如果你要处理的数据多了比如上百个城市的城市形态栅格或者多年份的地下水位数据再手动在QGIS里一个个点工具是真的浪费时间。建议把高频操作用Python封装起来核心库是rasterio、numpy、geopandas再加上一些GDAL的命令行调用。举个例子批量裁剪所有栅格到某个shp边界范围import rasterio from rasterio.mask import mask import geopandas as gpd shape gpd.read_file(region.shp) with rasterio.open(input.tif) as src: out_image, out_transform mask(src, shape.geometry, cropTrue) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(output_crop.tif, w, **out_meta) as dst: dst.write(out_image)注意在mask前确保shape和栅格使用相同坐标系如果坐标系不一致先调用to_crs统一。再比如批量栅格计算器公式本质上就是把两个rasterio打开的文件读成numpy数组做数组运算后写新tif。numpy数组的广播逻辑和QGIS栅格计算器几乎一样只要注意NoData掩膜的处理就行import numpy as np slope rasterio.open(slope.tif).read(1) ndvi rasterio.open(ndvi.tif).read(1) nodata_slope -9999 slope np.where(slope nodata_slope, np.nan, slope) result np.where((slope 30) (ndvi 0.2), 1, 0)这种脚本的最大优势是可复现今天跑一遍和半年后跑一遍参数一致结果就一致不会因为手动误操作产生差异。我对所有做空间分析的团队都建议至少把预处理环节脚本化它带来的效率提升在后续系列中会越来越明显。7. 最后的经验之谈在我做项目的时候经常有人问为什么同一份数据两个人跑出两张完全不同的图。答案往往不是操作不同而是他们在分析前对数据的理解不一样——一个检查过坐标系、处理过NoData、选择过合适的重采样方法另一个直接打开原始数据飞速点击工具。栅格空间分析这个领域前期的数据整理工作大约占整个项目一半时间这不是浪费时间而是给结果上保险。所以我的建议是不要急着追求工具数量和技术炫技先把每个栅格像元代表什么物理意义搞清楚再把分辨率、坐标系、NoData、波段信息这四个基础属性吃透。地基打牢了后面什么坡度坡向、距离分析、邻域统计、成本分析真的只是多了一个又一个按钮而已。最后再分享一个小技巧吧。我每次拿到一份新的栅格数据都会先在QGIS里用单波段伪彩色方式看一遍原始值同时打开直方图拖拽色带拉伸范围。这一个动作看起来不起眼但它能帮你快速发现像元值的异常分布比如有特别离谱的高值或者大量集中在某个区间这些信息会直接影响你后面怎么设阈值、怎么归类。栅格分析大多数时候不是算法不够强而是你对输入的“手感”不够准。这份手感只能靠一次次打开直方图、一次次试不同参数慢慢磨出来。