简介：本资源为2012—2020年NPP/VIIRS夜间灯光数据集，面向城市研究、遥感与地理信息分析人员，以及从事社会经济空间化建模的科研与教学用户。原始灯光影像经年度合成、去噪与连续性校正处理，形成可直接使用的长时间序列数据，适用于城市建成区提取、GDP空间化、人口分布及各类社会经济指标的空间化分析。压缩包共38个文件，约130.74MB，以9个tif栅格影像为主体，配套9个tfw坐标文件、9个ovr金字塔文件及11个xml元数据，便于在GIS软件中快速定位、显示与批量处理。目前已有1999人学习下载。数据覆盖2012至2020年逐年影像，时间连续、格式规范，可省去繁琐的预处理环节，帮助读者直接开展建成区扩张监测、经济指标空间化与区域发展对比等研究，适合作为中高级空间分析项目的底图数据。
1. 从一份 2012-2020 年 NPP/VIIRS 夜间灯光数据集压缩包说起
如果你手里正好有一个名为2012-2020年NPP/VIIRS夜间灯光数据集.zip的压缩包，或者正准备去下载这样一份数据，那你大概率已经知道它能干什么：用夜间灯光强度去代理经济活动、城市扩张、人口分布、灾后恢复，甚至做碳排放的代理变量。但真正让一线做空间分析的人头疼的，从来不是「有没有数据」，而是「这份数据拿到手之后，怎么从一堆 GeoTIFF 变成能进模型、能画图、能对比的干净面板」。NPP/VIIRS 夜间灯光数据集本身是 NASA 的 Suomi NPP 卫星上 VIIRS 传感器做的日/月合成产品，2012 年之后逐步替代了 DMSP/OLS，最大的好处是辐射定标更稳、没有 DMSP 那种饱和溢出，但代价是原始月度产品里混着杂散光、极光、火灾和油气燃烧的瞬时亮斑，直接拿来算城市灯光均值，结果会非常玄学。
这份 2012-2020 年的压缩包，通常覆盖的是月度或年度合成影像，空间分辨率在 15 弧秒（约 500 米）左右，格式以 GeoTIFF 为主，可能还带一些辅助的 CSV 或说明文档。它适合做长时间序列的省/市/县灯光均值提取、建成区范围识别、灯光基尼系数计算，也适合跟 GDP、人口栅格做回归。但前提是你得先过预处理这一关——热搜里「npp夜间灯光数据预处理」被反复搜，说明翻车的人不在少数。下面我就按自己实际处理这类数据的顺序，把选型、步骤、参数和踩过的坑一次讲清楚。
2. 先搞懂 NPP/VIIRS 月度与年度产品的差异，再决定用哪一层
2.1 月度产品、年度产品与「杂散光校正版」到底怎么选
打开压缩包，你可能会看到类似VNL_npp_2012-2020_global或者按年月分文件夹的结构。这里第一个要做的决策是：用月度还是年度？用原始版还是杂散光校正版？NPP/VIIRS 的月度产品（VCMSLCFG 之类）保留了最多的时相信息，适合做季节波动、短期事件（比如疫情初期灯光骤降）分析；但月度数据里每个像元的质量参差不齐，尤其是高纬度夏季的杂散光会污染整景。年度产品（VNL V2 年度合成）已经做了时序中值合成和部分杂散光剔除，拿来算年度均值更稳，但会抹掉年内变化。
我一般会这样选：如果研究目标是「2012-2020 年某区域灯光总量的趋势」，直接用年度杂散光校正版，省掉大量清洗工作；如果要做「某次灾害前后 3 个月的灯光变化」，那必须用月度，并且要自己写掩膜把杂散光区域抠掉。压缩包里如果同时有monthly和yearly两个目录，别偷懒只解压一个，先各拿一景在 QGIS 里对比一下直方图，你会看到月度数据的最大值能飙到几千甚至上万，而年度数据通常在几百以内——这不是数据错了，是月度没做离群值压制。
2.2 用 Python 快速检查压缩包内文件结构与坐标系
在动手写提取脚本之前，先花五分钟把压缩包里的文件清单和坐标系摸清楚。很多人直接gdalinfo都不跑，结果后面用错投影，算出来的面积差出几个数量级。下面这段代码用zipfile和rasterio做一次快速体检：
import zipfile import rasterio import os zip_path = "2012-2020年NPP/VIIRS夜间灯光数据集.zip" # 只列出压缩包内前 20 个文件，避免刷屏 with zipfile.ZipFile(zip_path, 'r') as z: names = z.namelist() for n in names[:20]: print(n) # 统计 tif 数量 tif_files = [n for n in names if n.lower().endswith(('.tif', '.tiff'))] print(f"GeoTIFF 总数: {len(tif_files)}") # 解压其中一个 tif 到临时目录做元数据检查 sample = tif_files[0] with zipfile.ZipFile(zip_path, 'r') as z: z.extract(sample, path="tmp_check") with rasterio.open(os.path.join("tmp_check", sample)) as src: print("CRS:", src.crs) print("分辨率:", src.res) print("范围:", src.bounds) print("波段数:", src.count) print("数据类型:", src.dtypes) print("NoData:", src.nodata)这段代码的逻辑很直白：先看压缩包内有没有按年份分目录，再确认 GeoTIFF 的坐标系是不是地理坐标（EPSG:4326）还是投影坐标。NPP/VIIRS 官方产品大多是 EPSG:4326，分辨率 15 弧秒，但有些二次分发的压缩包会被重投影成 Albers 或 Web Mercator，如果你不检查就直接按经纬度算面积，结果会偏。参数上重点看nodata，很多月度产品的背景值是 -999 或 255，如果你不设掩膜，算均值时会把背景算进去，灯光均值直接掉一个量级。这一步做完，你才能决定后面是用rasterio还是xarray来批量处理。
3. 用 Python 批量提取省/市/县灯光均值的完整流程
3.1 准备行政边界与统一投影
提取灯光均值的第一步不是写循环，而是把行政边界和栅格统一到同一个投影下。我通常用geopandas读 Shapefile，然后检查它的 CRS 是否和栅格一致。如果不一致，用to_crs转过去。注意：NPP/VIIRS 是地理坐标，直接算面积会得到平方度，所以如果你后面要算「灯光总量/面积」，最好先投影到等面积投影，比如 Albers（中国区域常用EPSG:6373或自定义 Albers）。但如果你只是算区域内灯光均值，地理坐标也能用，只是要注意像元面积随纬度变化。
import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np import pandas as pd # 读行政边界 gdf = gpd.read_file("boundaries/counties.shp") print("原始 CRS:", gdf.crs) # 读一景灯光栅格 with rasterio.open("tmp_check/" + sample) as src: raster_crs = src.crs print("栅格 CRS:", raster_crs) # 统一到栅格 CRS if gdf.crs != raster_crs: gdf = gdf.to_crs(raster_crs) # 只保留需要的区域，加速后续裁剪 gdf = gdf[gdf['province'].isin(['广东省', '广西壮族自治区'])]这里的关键参数是to_crs的目标 CRS 必须和栅格完全一致，否则mask会报错或者裁出空数组。另外，行政边界如果有飞地或岛屿，mask默认crop=True会裁掉外围，但如果你要保留完整边界形状，设crop=False。我一般会先做一次小范围测试，确认裁出来的像元数和边界面积大致匹配，再跑全量。
3.2 逐区域裁剪与均值计算的代码模板
下面这个函数是我反复用过的模板，核心是用rasterio.mask按几何裁剪，然后对有效像元求均值和总和。注意nodata的处理：先转成浮点，再把nodata设为np.nan，最后用np.nanmean。
def extract_light_stats(raster_path, geometry, nodata=None): with rasterio.open(raster_path) as src: # 如果栅格 nodata 未定义，用传入值 nd = nodata if nodata is not None else src.nodata try: out_image, out_transform = mask(src, [geometry], crop=True, nodata=nd) except ValueError: return None # 几何与栅格无交集 data = out_image[0].astype('float32') if nd is not None: data[data == nd] = np.nan # 有些产品用 0 表示背景，按需过滤 data[data <= 0] = np.nan valid = data[~np.isnan(data)] if valid.size == 0: return {'mean': np.nan, 'sum': np.nan, 'count': 0} return { 'mean': float(np.nanmean(valid)), 'sum': float(np.nansum(valid)), 'count': int(valid.size) } # 批量跑 results = [] for idx, row in gdf.iterrows(): geom = row.geometry stats = extract_light_stats("tmp_check/" + sample, geom) if stats: stats['county'] = row.get('name', idx) results.append(stats) df = pd.DataFrame(results) print(df.head())逻辑说明：mask返回的是裁剪后的数组和新的仿射变换，crop=True会紧贴几何边界，减少内存。data[data <= 0] = np.nan这一行是为了过滤掉背景 0 值，但要注意有些真实灯光极弱的区域可能真的是 0，如果你研究的是偏远地区，这行要慎用，可以改成只过滤nodata。参数上，nodata如果栅格自带就自动用，没有就手动传 -999。跑完一个区域后，最好把count和该区域的像元总数对比一下，如果count远小于预期，说明几何可能没落在栅格范围内，或者投影没对齐。
3.3 把 2012-2020 年所有月份拼成面板数据
单景提取只是热身，真正的活是循环 2012 到 2020 所有月份，再把结果拼成county-year-month的面板。这里有两个坑：一是文件命名不统一，有的用201201，有的用2012_01；二是不同年份的栅格范围可能略有偏移，导致同一个县在边缘年份被裁掉。我的做法是先按文件名解析出年月，排序后逐个提取，最后用pandas.concat合并，并对缺失值做前后向填充标记。
import re import glob all_files = sorted(glob.glob("tmp_check/*.tif")) records = [] for f in all_files: # 从文件名提取年月，兼容 201201 和 2012_01 m = re.search(r'(20\d{2})[_\-]?(0[1-9]|1[0-2])', f) if not m: continue year, month = int(m.group(1)), int(m.group(2)) for idx, row in gdf.iterrows(): stats = extract_light_stats(f, row.geometry) if stats: stats.update({'county': row.get('name', idx), 'year': year, 'month': month}) records.append(stats) panel = pd.DataFrame(records) panel.to_csv("nightlight_panel_2012_2020.csv", index=False, encoding='utf-8-sig') print(panel.shape)这段代码跑起来可能比较慢，因为每个县每景都要做一次mask。优化办法是先把所有栅格用rioxarray读成一个DataArray，然后用rasterstats的zonal_stats批量算，速度能快几倍。但rasterstats对nodata的处理不如手写灵活，如果你数据里杂散光很多，还是建议手写掩膜。参数上，encoding='utf-8-sig'是为了 Excel 打开不乱码，这个细节很多人忽略，结果给合作者发 CSV 被吐槽。
4. 避坑与排查：NPP/VIIRS 预处理里最容易翻车的 5 件事
4.1 现象：灯光均值逐年下降，但经济数据在涨
原因：月度产品里的杂散光在早期年份（2012-2014）污染更严重，尤其是夏季高纬度区域，导致早期均值被抬高，后期校正后反而显得下降。另外，如果你没做离群值压制，个别油气燃烧的亮斑（比如中东、西伯利亚）会把区域均值拉高。解决：用年度杂散光校正版做趋势分析，或者在月度数据里先做 99 分位数截断，把超过阈值的像元用邻域中值替换。我一般会先画一个全国灯光总值的时序图，如果 2012 到 2014 有个明显的台阶，基本就是杂散光没清干净。
4.2 现象：裁剪出来的区域全是 NoData
原因：行政边界的坐标系和栅格不一致，或者边界文件是经纬度但栅格被重投影过。另一个常见原因是mask的invert参数用反了，把内部裁掉了。解决：先打印gdf.crs和src.crs，确认一致；再用gpd.clip做一次可视化，看边界是否落在栅格范围内。如果边界是跨 180 度经线的，还要处理经度环绕问题，这个在 NPP/VIIRS 全球产品里偶尔遇到。
4.3 现象：同一区域不同月份的像元数量差异巨大
原因：不同月份的栅格范围或分辨率被二次分发者改过，或者你用的 Shapefile 有简化版本，在不同月份裁剪时边界略有出入。解决：统一用同一份边界文件，并且在提取前先用rasterio.warp.reproject把所有栅格重采样到同一网格。重采样方法选bilinear或cubic，不要用nearest，否则灯光值会跳变。重采样后像元数量就一致了，面板数据也更好对齐。
4.4 现象：算出来的灯光总和是负数
原因：nodata值没设对，比如原始数据用 -999 表示背景，但你用 0 去过滤，结果 -999 被当成有效值参与求和。解决：在rasterio.open后立刻读src.nodata，如果为None，用src.read(masked=True)让 rasterio 自动掩膜。或者手动data[data < 0] = np.nan，因为灯光辐射值不可能为负。这个坑我踩过两次，第一次是算某省灯光总量得到负值，排查半天才发现是背景值没清。
4.5 现象：面板数据里某些县某些年份整段缺失
原因：压缩包里某些年份的文件缺失，或者文件名解析正则没匹配上。解决：先ls一遍所有文件名，用pandas生成一个完整的county-year-month笛卡尔积，再和提取结果左连接，缺失的标记出来。如果缺失是数据源本身没有，那就只能接受，并在论文里说明。如果是解析问题，调整正则，比如兼容VNL_npp_2015_01.tif和201501.tif两种命名。
5. 进阶技巧：用灯光数据做建成区提取与跨传感器校正
5.1 用阈值法快速提取建成区范围
拿到干净的年度灯光栅格后，一个常见需求是提取建成区。最简单的方法是阈值法：先算区域灯光均值和标准差，取mean + 1.5 * std作为阈值，高于阈值的像元判为建成区。这个方法在 NPP/VIIRS 上比 DMSP 好用，因为 VIIRS 没有饱和，城市核心和郊区的梯度更明显。但阈值不是固定的，不同区域差异很大，我一般会先画几个典型城市的灯光剖面图，手动定一个初始阈值，再用scipy.ndimage做形态学闭运算，把零散像元连成片。
from scipy import ndimage def extract_urban(data, k=1.5): valid = data[~np.isnan(data)] if valid.size == 0: return np.zeros_like(data, dtype=bool) thresh = np.nanmean(valid) + k * np.nanstd(valid) binary = data > thresh # 闭运算连接邻近像元 binary = ndimage.binary_closing(binary, structure=np.ones((3,3))) return binary参数k控制建成区范围，k越大越保守。我试过k=1.0到2.0，一般1.5比较接近目视解译。闭运算的结构元素用 3x3 就够了，太大容易把农村灯光也吞进去。提取完可以叠加行政边界算建成区面积占比，和统计年鉴里的建成区面积对比，如果差太多，回头调k。
5.2 跨传感器校正：让 NPP/VIIRS 和 DMSP/OLS 能接上
如果你的研究要往前延伸到 2012 年之前，就必须把 NPP/VIIRS 和 DMSP/OLS 做校正。常见做法是找重叠年份（2012-2013），在同一个区域分别算两种数据的灯光总和，拟合一个线性或幂函数关系，然后把 DMSP 数据映射到 VIIRS 尺度。注意 DMSP 有饱和，城市核心的 DN 值封顶在 63，所以校正前要先把 DMSP 的饱和像元剔除或做去饱和处理。我一般用numpy.polyfit做一次多项式拟合，但只对非饱和区域拟合，饱和区域单独用经验公式。这个步骤没有标准答案，不同区域的拟合参数不一样，所以一定要在论文里报告你的拟合方程和 R²。
5.3 一个我坚持了很久的习惯
每次拿到新的 NPP/VIIRS 压缩包，我不会直接跑全量提取，而是先挑一个自己熟悉的城市（比如我常选广州），把 2012 到 2020 每一年的灯光均值算出来，画一条折线。如果这条线平滑上升，说明数据质量可控；如果中间有断崖或尖峰，那就得回去查那一年那一个月的原始影像。这个习惯帮我省了很多后悔药——有一次发现某年 7 月的数据整体偏低，排查后发现是压缩包里那个月的文件其实是半成品，只有半景。希望帮到你。
本文还有配套的精品资源，点击获取
