简介本资源为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(fGeoTIFF 总数: {len(tif_files)}) # 解压其中一个 tif 到临时目录做元数据检查 sample tif_files[0] with zipfile.ZipFile(zip_path, r) as z: z.extract(sample, pathtmp_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默认cropTrue会裁掉外围但如果你要保留完整边界形状设cropFalse。我一般会先做一次小范围测试确认裁出来的像元数和边界面积大致匹配再跑全量。3.2 逐区域裁剪与均值计算的代码模板下面这个函数是我反复用过的模板核心是用rasterio.mask按几何裁剪然后对有效像元求均值和总和。注意nodata的处理先转成浮点再把nodata设为np.nan最后用np.nanmean。def extract_light_stats(raster_path, geometry, nodataNone): 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], cropTrue, nodatand) 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返回的是裁剪后的数组和新的仿射变换cropTrue会紧贴几何边界减少内存。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, indexFalse, encodingutf-8-sig) print(panel.shape)这段代码跑起来可能比较慢因为每个县每景都要做一次mask。优化办法是先把所有栅格用rioxarray读成一个DataArray然后用rasterstats的zonal_stats批量算速度能快几倍。但rasterstats对nodata的处理不如手写灵活如果你数据里杂散光很多还是建议手写掩膜。参数上encodingutf-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(maskedTrue)让 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, k1.5): valid data[~np.isnan(data)] if valid.size 0: return np.zeros_like(data, dtypebool) thresh np.nanmean(valid) k * np.nanstd(valid) binary data thresh # 闭运算连接邻近像元 binary ndimage.binary_closing(binary, structurenp.ones((3,3))) return binary参数k控制建成区范围k越大越保守。我试过k1.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 月的数据整体偏低排查后发现是压缩包里那个月的文件其实是半成品只有半景。希望帮到你。本文还有配套的精品资源点击获取