ArcGIS夜间灯光数据提取与强度统计实战:从掩膜提取到分区统计全流程解析
2026/9/23 5:56:33 网站建设 项目流程

简介:在地理信息系统与遥感技术中,栅格数据的区域统计分析是空间信息提取的关键环节。不同来源的夜间灯光数据(如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 = r"D:\night_light_data" arcpy.env.overwriteOutput = True # 行政边界和灯光栅格列表 boundary = r"D:\gis_data\taiyuan_boundary.shp" field = "name" light_rasters = arcpy.ListRasters("*.tif") for ras in light_rasters: # 输出掩膜栅格 out_mask = r"D:\night_light_data\mask_" + ras out_table = r"D:\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(r"D:\data\viirs.tif").spatialReference shp_sr = arcpy.Describe(r"D:\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 字段获得),得到“单位面积亮度”,如果这个值在相邻区域之间存在突变,尤其是没有山脉河流阻隔的地方,往往说明边界提取或掩膜过程中出了问题。最后,把行政区面要素转成栅格,再和灯光掩膜相减,可以直观看出哪些位置的灯光被多扣或少扣。这套查错顺序做下来,至少能过滤掉九成的低级错误。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询