简介在地理信息系统与遥感技术中栅格数据的区域统计分析是空间信息提取的关键环节。不同来源的夜间灯光数据如DMSP-OLS与VIIRS具有不同的分辨率与像元值特征而行政区划边界与栅格影像之间的坐标系一致性直接影响提取结果的准确性。通过掩膜提取实现研究区裁剪再借助分区统计工具计算均值、总和等指标是城市灯光强度研究的标准流程。针对多期数据处理利用ArcPy脚本可显著提升批量处理效率。本文以ArcGIS平台为例系统梳理夜间灯光数据预处理、像元大小控制、统计口径选择及专题图制作中的常见问题帮助工程与科研人员快速获得可靠的城市灯光特征参数。1. 夜间灯光提取不只是“裁剪”这么简单夜间灯光数据现在已经是城市规划、经济估算和人口分布研究的常客但真正在 ArcGIS 里把一整幅全球栅格变成某个县级市的分区亮度统计表中间隔着三道容易翻车的工序区域选取、掩膜提取、分区统计。很多人以为用“按掩膜提取”把行政边界以外的像元砍掉就算完事结果拿去和统计年鉴对数字时均值偏了好几倍。问题往往出在坐标系没有统一、像元类型不对、或者连接字段选错。这篇内容基于我实际拆过的一套作业流程把夜间灯光的提取和强度计算完整串起来包含 ArcGIS 界面操作、参数选择逻辑以及用 ArcPy 批量复现的脚本。新手能照步骤跑通老手可以对比一下自己在分区统计时有没有漏掉重采样设置。2. 数据准备从 VIIRS/DMSP-OLS 到统一坐标系2.1 夜间灯光数据的两个主流来源与选择夜间灯光栅格目前最常用的是两个系列DMSP-OLS 和 VIIRS。DMSP-OLS 数据时间跨度是 1992—2013 年像元亮度值范围是 0—63空间分辨率约 1km容易饱和城市中心亮度值会被“顶满”。VIIRS 是新一代传感器2012 年之后的数据更细腻亮度值没有 63 的上限而且去掉火光等杂散光源的月合成产品做得更干净。如果你研究的是近十年的事优先找 VIIRS VNL 月度或年度合成数据如果做历史长序列对比必须用 DMSP-OLS 且要做连续性校正。二者在 ArcGIS 里的处理逻辑一样但 VIIRS 的浮点型像元值在“以表格显示分区统计”时建议选用“均值”而不是“总和”因为不同月份的像元值基数不同总和容易受面积影响。下载时要注意文件名里的分辨率标识。有的平台提供“viirs_2022_vcmcfg”这种版本VCMCFG 是排除杂散光版本更适合城市灯光研究。拿到数据后先别急打开属性表看像元类型——大部分是整型少数是浮点型。这会影响后面的统计口径我一般会在原始 tif 上右键看“源”里的像素类型记录下来。2.2 ArcGIS 中行政矢量与栅格坐标系的统一这是一个看起来基础但最容易埋雷的环节。行政边界 shp 常用的是 CGCS2000 或 WGS84 地理坐标系而灯光栅格很多是 WGS84 或者 Sinusoidal 投影。直接拖进 ArcMap如果软件启用了“后台地理处理”的自动投影界面可能会给你“假装”对齐但实际掩膜提取时输出的范围可能跑偏。我的习惯是先把矢量数据和栅格数据都检查一遍在内容列表里右键每个图层打开“属性→源”看“空间参考”那一栏。如果两者的坐标系名称不一致哪怕都是地理坐标只是基准面不同比如 WGS84 和 CGCS2000 之间差几十厘米到一米对整个城市尺度的灯光提取来说误差不大但为了严谨我会用“投影”或“定义投影”把行政矢量统一到灯光栅格的坐标系。操作路径是ArcToolbox → 数据管理工具 → 投影与变换 → 要素 → 投影输入 shp选择输出坐标系时直接从下拉列表里选“与图层相同”或者手动导入灯光 tif 的坐标系。如果栅格是投影坐标系行政矢量是地理坐标系直接用“按掩膜提取”也能跑但输出栅格的像元大小可能被默认重新采样成不适合的值。因此建议统一坐标系后再记录一下灯光栅格的像元大小比如 15 弧秒或 0.004166 度后续掩膜提取时手动指定。2.3 用按属性选择提取目标行政区含多条件写法在 ArcMap 里打开行政矢量 shp点击“选择”→“按属性选择”。双击字段名“name”点“获取唯一值”能看到所有区域名。单选一个区域时条件写成name 太原多个区域时用 Or 拼接name 太原 Or name 晋中 Or name 吕梁点击“应用”后选中的要素会高亮。这时右键图层选择“数据”→“导出数据”导出范围选“所选要素”输出要素类命名成 taiyuan_boundary.shp。注意保存路径不要包含中文否则后续有些工具会报“数据集不存在”的莫名其妙的错。导出后最好再在内容列表里把原始 shp 的勾选去掉只保留导出图层这样后面做掩膜时不会选错输入。这段操作的逻辑是直接从全国或全省的行政区划里切出研究区边界而不是手动裁剪。好处是边界完全来自权威数据不会因为你手工描边画歪了导致灯光像元被多算或少算。如果你研究的是乡镇尺度同样用这个办法只是字段名可能是“XZQMC”这类拼音缩写。建议在导出前先打开属性表确认字段名称避免把“PAC”这种代码字段当成名称字段。3. 掩膜提取与强度统计核心参数决定结果差异3.1 按掩膜提取的工具参数与输出设置拿到研究区边界后开始执行夜间灯光的裁剪。打开 ArcToolbox → Spatial Analyst 工具 → 提取分析 → 按掩膜提取。这里输入的“输入栅格”是原始夜间灯光 tif“输入栅格数据或要素掩膜数据”是刚才导出的行政边界。关键参数是“输出栅格”的位置和名称。我会把输出名写成 mask_light.tif 并放在专门的工作目录里。很多人忽略的一个设置是环境变量。在“按掩膜提取”对话框左下角有“环境”按钮点开后看“处理范围”和“栅格分析”里的“像元大小”。如果这里不设置结果可能比原始栅格的范围小一圈或者像元大小被自动改成和默认值不同的值。我的固定做法是处理范围选“与掩膜相同”像元大小选“与输入栅格相同”。这样能保证提取出来的像元行列数正确不会因为重采样造成亮度值被平滑。还有一个坑如果掩膜要素是多边形且内部有孔洞输出栅格在孔洞位置会是 NoData。这在实际计算区域均值时会被忽略但如果后续需要把灯光栅格转成面或点参与其他分析NoData 会造成边界不连续。这时候可以先用“提取分析→按属性提取”或者对掩膜要素做“融合”后再执行。大多数行政区是单面要素不会遇到孔洞问题但含飞地的区划就要注意。执行完后内容列表里会出现两个同名图层一个是原始 tif 的裁剪结果一个是 ArcMap 自动添加的颜色带渲染。建议右键输出栅格打开属性查看“源”里的像元大小和空间参考确认和预期一致。如果发现像元大小变成 0.01 度之类的值说明环境设置没生效需要重新执行。3.2 以表格显示分区统计字段与统计类型怎么选灯光强度计算的核心工具是“Spatial Analyst 工具 → 区域分析 → 以表格显示分区统计”。这里有四个必填项输入栅格数据或要素区域数据、区域字段、输入赋值栅格、输出表。区域数据就是行政区边界 shp比如 taiyuan_boundary.shp。区域字段要看属性表里的唯一标识字段常见的是“name”或“OBJECTID”。我建议用“name”因为后续连接时不容易因为 OBJECTID 顺序变化而错位。输入赋值栅格是掩膜提取后的 mask_light.tif注意要选裁剪后的不要选原始全球数据否则统计结果会把所有区域都算一遍而且属性表会巨大。统计类型默认是“SUM”。下拉框里有 MEAN、MINIMUM、MAXIMUM、RANGE、STD 等。对于夜间灯光强度MEAN 是普遍需要的指标因为它反映区域平均灯光亮度总和 SUM 也常用但受区域面积影响大适合做总量比较。如果你的研究是“灯光强度与 GDP 关系”建议同时输出 SUM 和 MEAN 两张表或者直接输出一个包含多种统计类型的表。工具里只能选一种统计类型要么分两次执行。我通常先跑 MEAN再跑 SUM然后通过连接合并。输出表是 dbf 格式。字段结构大概是VALUE区域字段、AREA、COUNT以及你选择的统计字段。COUNT 表示参与统计的像元个数可以用来计算有效灯光面积占比。这里有个细节如果掩膜提取后有些像元是 NoData它们不会计入 COUNT。所以当区域边缘有大量 NoData 时MEAN 可能会略微偏大因为它在计算时只用了有效像元。3.3 一个 ArcPy 批处理脚本示例如果只需要做一次界面操作足够。但当你面对几十个年份的灯光数据或者需要把每个县单独提取并计算就要用脚本批量处理。下面这段我常用的 ArcPy 脚本核心是循环遍历 tif 文件用同一个行政边界做掩膜和分区统计。import arcpy from arcpy.sa import * # 设置工作空间和覆盖选项 arcpy.env.workspace rD:\night_light_data arcpy.env.overwriteOutput True # 行政边界和灯光栅格列表 boundary rD:\gis_data\taiyuan_boundary.shp field name light_rasters arcpy.ListRasters(*.tif) for ras in light_rasters: # 输出掩膜栅格 out_mask rD:\night_light_data\mask_ ras out_table rD:\night_light_data\stats_ ras.replace(.tif, .dbf) # 按掩膜提取 arcpy.gp.ExtractByMask_sa(ras, boundary, out_mask) # 分区统计统计类型为 MEAN arcpy.gp.ZonalStatisticsAsTable_sa(boundary, field, out_mask, out_table, DATA, MEAN) print(完成: ras)这段脚本里arcpy.gp.ExtractByMask_sa是“按掩膜提取”的命令行接口arcpy.gp.ZonalStatisticsAsTable_sa是分区统计接口。注意DATA参数表示忽略 NoData如果改成NODATA则会把 NoData 当作 0 参与统计这在夜间灯光里会导致 MEAN 被严重拉低尤其是郊区大量像元本身是 0 值时。分区统计里的 MEAN 只统计有效像元所以DATA是正确的选择。脚本方案的优势是跨年份的处理口径完全一致。但你必须先确认所有 tif 文件的坐标系、像元大小一致否则不同年份的统计结果不具备可比性。如果发现不一致需要在循环里先调用arcpy.Resample_management统一像元大小再用arcpy.ProjectRaster_management统一投影。4. 从统计表到专题图连接、符号化与出图4.1 连接和关联为什么分区统计表会“对不上”分区统计生成的是 dbf 表需要把它连接到行政区矢量图上。右键行政边界图层选择“连接和关联”→“连接”。第一个下拉框选“某个表的字段”第二个选刚才生成的 dbf 表然后指定两个表关联的字段。行政边界这边选“name”dbf 表那边选“VALUE”——默认生成的分区统计表第一列就叫 VALUE存放的是区域字段值。但这里有一个常见的坑如果行政边界属性表里有重名区域比如两个乡镇都叫“城关镇”连接后会变成一对多结果只有第一条记录被连接其他记录显示为 NULL。解决的办法是在分区统计之前先给行政边界添加一个唯一 ID 字段比如用OBJECTID然后用OBJECTID作为区域字段执行统计连接时也用它。虽然图面上不直观但为了数据准确我一般都会新建一个整型字段FID2用字段计算器赋值为OBJECTID再拿FID2做区域字段。连接后建议把连接结果导出成新要素类避免会话关闭后连接丢失。操作方法右键图层 → 数据 → 导出数据 → 导出为 shp。导出后再打开属性表能看到所有字段包括 MEAN、SUM、COUNT。此时可以再用“字段计算器”根据 COUNT 和像元大小计算有效面积有效面积_km2 COUNT * 像元大小_km2如果像元大小是 0.004166 度需要先转换成米。更稳妥的方式是在投影坐标系下统计让像元面积直接是平方米。4.2 分级色彩符号化手动分段与自然间断点连接好之后右键图层 → 属性 → 符号系统 → 数量 → 分级色彩。字段选 MEAN配色建议用“黄-橙-红”或“深蓝-浅蓝”因为夜间灯光本质上是强度由暗到亮。关键在“分类”按钮。ArcGIS 默认用“自然间断点Jenks”分级它会让类内方差最小适合展示空间差异。但如果你想做多年对比图不能用每次自动算出的间断点否则颜色深浅含义不同跨年份看图会误导。我的做法是固定分级阈值比如 0、5、10、20、40、60。在“分类”对话框里选择“手动”然后输入断点。这样每一年的地图颜色深浅含义一致直接放一起对比。对于 DMSP-OLS 数据0—63 的亮度区间本身不大5 级左右就够对于 VIIRS 数据亮度可能到几百建议先用直方图观察数据分布再决定断点。如果大多数像元集中在 0—5而少数城市中心是 100 以上自然间断点会把 0—5 切成好几段反而弱化了城乡差异。我会优先用分位数或者手动按对数间隔分确保低亮度区有梯度。4.3 布局视图中的图例、比例尺与指北针专题图最终要放到布局视图里。点击左下角的“布局视图”按钮可以看到纸张模型。插入图名、图例、比例尺、指北针在“插入”菜单里都有。图例默认会把所有子图层的名称都列出来很啰嗦。双击图例在“项目”选项卡里删掉不需要的只保留“MEAN”这一项。比例尺类型建议选“交替单位”比如公里和英里都显示。指北针选一个简洁的正北方向即可。这里要提到的细节是出图精度。灯光强度专题图通常还要叠加行政边界线右键行政边界图层 → 属性 → 显示 → 勾选“符号级别”把边界线放在填充色上面设置白色细线或深灰色线区分度更高。导出图片时使用“文件 → 导出地图”分辨率设到 300 dpi格式选 TIFF 或 PNG避免期刊投稿时图太小。5. 实战中的几个坑坐标系、像元大小与统计口径5.1 坐标系不一致导致“空白掩膜”有一次我直接拿 WGS84 的灯光栅格去对 CGCS2000 的县界做掩膜工具没报错但输出栅格全黑。后来检查发现两个数据的坐标系虽然都叫“地理坐标”但基准面不同边界在几百公里外偏移掩膜范围完全落在灯光数据之外。解决方法是先在 ArcToolbox 里用“投影”把矢量转成 WGS84或者用“投影栅格”把灯光转到 CGCS2000。注意不要用“定义投影”因为那只是修改元数据不会真的改变像元位置。判断坐标系是否真正的统一一个简单的办法是把灯光和矢量都拖进 ArcMap开启“视图中 → 数据框属性 → 坐标系”查看数据框的投影然后缩放至两个图层的交集。如果边界线和灯光亮区有肉眼可见的错位说明坐标系不匹配。另一种是直接用 ArcPy 读取空间参考比较import arcpy rast_sr arcpy.Describe(rD:\data\viirs.tif).spatialReference shp_sr arcpy.Describe(rD:\data\boundary.shp).spatialReference print(rast_sr.name, shp_sr.name)如果两者名称不同就执行投影转换。坐标统一不是可选项是前提条件。5.2 像元大小与重采样对强度均值的影响夜间灯光栅格的分辨率从 0.004 度到 0.01 度不等。当你用行政边界做掩膜后如果 ArcGIS 环境变量里“像元大小”不是“与输入相同”它可能会自动用边界要素的分辨率去重采样。比如边界数据没有像元大小概念系统可能会给一个默认值比如 0.001 度导致输出栅格比原始栅格格子更细产生大量重复插值。插值后的 MEAN 可能变化不大但 SUM 会因为你把每个像元面积算错而失真。正确做法是分区统计前在“环境”里强制把“栅格分析”的“像元大小”设置为原始灯光栅格的值。重采样方法也值得注意。如果需要重新投影灯光数据默认使用“双线性”或“三次卷积”会平滑亮度值如果使用“最邻近”则保持原始亮度值但边缘锯齿明显。对于灯光这种连续型亮度变量双线性更合理可以避免高空值被低估。但如果你后续要做“灯光面积”统计关注的是亮像元的个数那么用最邻近法会更保守不会因为平滑把 0 值周围的低值像素变成有值像素。5.3 灯光强度口径均值/总和适合什么场景分区统计的 MEAN 和 SUM 各有适用场景。研究城市化强度或夜间活动水平MEAN 更合适因为它不受区域面积影响能直接比较不同大小的行政区。比如太原市的 MEAN 和晋中市的 MEAN 是可以直接比的而 SUM 会明显偏向面积更大的区域。但 SUM 在城市总体经济规模估算里很有用因为灯光总量可以近似区域活动总量。实际操作里我会同时保留两个字段论文里如果需要“单位面积灯光强度”就用 MEAN需要“总灯光辐射量”就用 SUM。此外还有一个容易被忽略的指标COUNT 字段对应的像元个数结合像元大小可以反推“有效灯光覆盖面积”。当你需要计算“建成区范围内灯光占比”时这个值比单纯看 MEAN 更稳定。具体做法是用属性选择找出 COUNT 大于 0 的区域然后求和。这个指标在长时间序列分析中能规避传感器增益变化带来的亮度值漂移问题。5.4 统计结果的验证方法拿到分区统计表后不要直接信输出。我习惯随机抽两三个行政区用“多值提取至点”工具把灯光栅格的像元值提取到随机点上然后对比点均值和分区表里的 MEAN。如果相差不大10% 以内说明统计路径没问题。另一个验证是看总和 SUM 除以区域面积从矢量属性表的 Shape_Area 字段获得得到“单位面积亮度”如果这个值在相邻区域之间存在突变尤其是没有山脉河流阻隔的地方往往说明边界提取或掩膜过程中出了问题。最后把行政区面要素转成栅格再和灯光掩膜相减可以直观看出哪些位置的灯光被多扣或少扣。这套查错顺序做下来至少能过滤掉九成的低级错误。本文还有配套的精品资源点击获取
