简介本资源为2012—2020年NPP-VIIRS夜间灯光遥感数据VNL2的预处理成果面向从事城市扩张、经济活动评估、碳排放与人口空间化等研究的高校师生及科研人员帮助其跳过繁琐的裁剪与异常值处理环节直接获得可直接用于建模分析的全国尺度数据。数据分辨率500m已用WGS84中国矢量地图完成全国范围裁剪并转换为Asia_Lambert_Conformal_Conic平面投影坐标系同时参考Kang Wu等2019方法以472.86为全国像元辐射阈值剔除极端异常值并去除负值像元。压缩包共36个文件以9个tif栅格影像为主配套tfw坐标文件、xml元数据与ovr金字塔文件整体约176.73MB按年份命名便于逐年调用。目前已有3505人学习下载适合需要长时序、已清洗夜光数据支撑空间分析的研究者直接使用。1. 预处理后的2012-2020年NPP-VIIRS夜光遥感数据从原始栅格到能直接跑模型的干净面板如果你拿到的夜光遥感数据是2012–2020年NPP-VIIRS的“预处理后”版本那它大概率已经替你干掉了三件最磨人的事去除了火山、油气田燃烧等短暂异常光源滤掉了背景噪声做了跨年辐射定标。但别急着直接丢进模型——我见过太多人拿到所谓“预处理后”的数据第一反应是打开看数值范围发现最小值是负数最大值能到几千然后就开始怀疑人生。这个标题真正指向的是一套能让你跳过清洗环节、直接进入时空分析或社会经济指标回归的年度合成夜光面板。它适合做城市扩张监测、GDP或人口空间化、灾后恢复评估的人也适合想拿夜光当协变量但不想在数据清洗上耗三个月的新手。核心价值就一句话把“能用的夜光”从原始月度栅格变成年度、去噪、可比的面板你只需要关心怎么用而不是怎么洗。2. 先搞懂预处理到底动了哪些刀年度合成、异常剔除与跨年定标2.1 从月度到年度为什么不是简单求平均NPP-VIIRS的原始产品是月度无云合成但月度数据里混着两类东西一类是稳定的人类活动光源另一类是渔船灯光、天然气火炬、火山岩浆这些“闪一下就走”的短暂光源。如果你直接对12个月求平均那些只在某个月出现的强光源会被摊薄但不会消失反而会在年均值里留下一个不高不低的“鬼影”。常见做法是先做月度层面的异常像元标记再在年度合成时用中位数或修剪均值。我一般会看年度合成用的是哪种统计量如果是中位数对极端值更稳如果是均值那就要确认异常剔除是否彻底。预处理后的数据通常会附带一个质量波段告诉你每个像元参与合成的有效月份数这个波段比夜光值本身还重要——有效月份少于6的像元年均值基本不可信。2.2 跨年定标2012和2020的数值为什么不能直接比NPP-VIIRS的传感器在2012年刚上线时辐射定标参数和后续年份不一样导致2012年的夜光值整体偏低2013年之后才逐渐稳定。如果你拿2012年的均值和2020年的均值直接做差值得到的“增长”里有一部分是传感器差异不是真实的光源变化。预处理后的数据一般会做两种处理之一要么把2012年数据通过线性模型调整到后续年份的基准要么直接建议你从2013年开始用。我见过一个做城市扩张的案例没注意这点结论是“2012到2013年全国夜光暴涨”后来发现是定标跳变。所以拿到数据第一件事是画一条全国或区域总夜光的时间序列看2012到2013之间有没有一个不自然的台阶。如果有要么剔除2012要么找数据说明里的定标系数自己调。2.3 负值和背景噪声为什么你的最小值不是零夜光遥感反映的是辐射亮度经过大气校正和背景噪声去除后理论上无光区域应该接近零但实际数据里经常出现负值。这不是错误而是噪声去除算法比如阈值法或掩膜法的副产品。预处理后的数据通常会把负值截断为零或者保留负值但告诉你“小于零的像元视为无光”。如果你要做对数变换或比值计算负值会直接让程序报错。我一般会先统计负值像元的比例如果低于1%直接截断为零如果高于5%就要怀疑背景噪声去除是否过度可能把微弱但真实的光源也压成了负值。这个比例在数据说明里通常不会写得自己算。3. 用Python把年度夜光面板读进来栅格拼接、重采样与区域统计3.1 读取与拼接把分幅数据拼成一张全国图预处理后的NPP-VIIRS数据通常按经纬度分幅一年可能有多个瓦片。如果你只研究一个城市读对应瓦片就行但如果做全国或区域分析就得先拼接。我一般用rasterio读用rioxarray做拼接因为它能自动处理坐标参考系和nodata。下面这段代码假设你有一年的多个GeoTIFF瓦片放在同一个文件夹里。import rioxarray import xarray as xr from pathlib import Path # 假设瓦片按年份放在 data/2015/ 下文件名类似 NPP_2015_tile_*.tif tile_dir Path(data/2015) tiles sorted(tile_dir.glob(*.tif)) # 用rioxarray打开每个瓦片沿x和y方向拼接 # 注意如果瓦片有重叠merge会取第一个非nodata值建议先检查重叠区 arrays [rioxarray.open_rasterio(t, maskedTrue) for t in tiles] merged xr.concat(arrays, dimband) # 先沿band堆叠 # 实际拼接用combine_by_coords更稳但要求瓦片坐标不重叠 merged rioxarray.merge.merge_arrays(arrays) # 写出拼接后的单年栅格 merged.rio.to_raster(data/2015_merged.tif)逻辑说明rioxarray.open_rasterio带maskedTrue会把nodata转成NaN避免后续统计被污染。merge_arrays会自动对齐坐标但如果瓦片之间有重叠且值不一致它会取第一个所以拼接前最好用gdalinfo或rioxarray检查一下重叠区。参数上maskedTrue是关键不加的话nodata会参与均值计算结果偏小。3.2 重采样到统一网格为什么不能直接算区域均值不同年份的瓦片可能因为传感器扫描角度或处理流程网格原点有微小偏移。如果你直接按行政边界做区域统计边界像元可能落在不同网格上导致年际比较出现伪变化。常见做法是先重采样到一个统一的参考网格比如1km或500m的等经纬度网格。我一般用rioxarray的reproject_match以某一年的栅格为基准把其他年份对齐过去。import rioxarray # 以2015年为基准网格 base rioxarray.open_rasterio(data/2015_merged.tif, maskedTrue) target rioxarray.open_rasterio(data/2020_merged.tif, maskedTrue) # 把2020重采样到2015的网格用双线性插值 # 夜光值连续双线性比最近邻更合适 target_aligned target.rio.reproject_match(base, resampling2) # 2对应双线性 target_aligned.rio.to_raster(data/2020_aligned.tif)逻辑说明reproject_match的resampling参数用数字2代表双线性适合连续型夜光值如果用最近邻0会在边界产生锯齿。注意重采样会引入平滑如果你研究的是点状光源如单个工厂双线性可能把光晕扩大这时候用最近邻更保真。参数选择取决于你的分析尺度区域统计用双线性点源识别用最近邻。3.3 按行政区做区域统计用zonal stats提取面板拼接和重采样之后下一步是按行政区或自定义网格提取夜光总和、均值、最大值。我一般用rasterstats的zonal_stats它支持GeoJSON或Shapefile边界。下面以某虚拟区域的GeoJSON为例。from rasterstats import zonal_stats import geopandas as gpd # 读取行政边界假设是WGS84经纬度 gdf gpd.read_file(data/region_boundary.geojson) # 对每个多边形计算夜光总和、均值、有效像元数 stats zonal_stats( data/2015_merged.tif, gdf.geometry, stats[sum, mean, count], nodata-9999, # 根据数据说明改常见是-9999或0 geojson_outFalse ) # 把结果挂回GeoDataFrame gdf[nightlight_sum_2015] [s[sum] for s in stats] gdf[nightlight_mean_2015] [s[mean] for s in stats] gdf[valid_pixels_2015] [s[count] for s in stats] gdf.to_file(data/region_nightlight_2015.geojson, driverGeoJSON)逻辑说明nodata参数必须和数据里的nodata值一致否则会把nodata当成真实夜光算进去。count返回的是非nodata像元数如果某个区域的count远小于总面积对应的像元数说明该区域有大量nodata统计值不可信。我一般会设一个阈值比如有效像元占比低于70%的区域就标记为“数据不足”不参与后续年际比较。4. 避坑与排查预处理后数据最常见的5个翻车点4.1 现象2012年总夜光比2013年高出一截原因定标跳变未完全消除解决画时间序列确认后剔除2012或自行调平这个坑我踩过两次。第一次是做全国夜光总量趋势发现2012到2013是下降的但所有经济数据都说那几年在增长。后来把各省分开画发现有的省2012高、有的省2013高说明不是真实变化而是定标残留。预处理后的数据虽然号称做了跨年定标但不同区域的调整系数可能不一样尤其是高纬度地区。解决办法很简单把2012和2013的年度栅格相减看差值空间分布。如果差值在空间上随机说明定标没问题如果差值呈现明显的南北或东西梯度说明定标不彻底。这时候要么从2013年开始用要么用2013–2015的稳定期做基准把2012线性映射过去。4.2 现象区域夜光均值逐年下降但总和在上升原因有效像元数变化解决用总和而非均值做趋势或固定像元集这个坑更隐蔽。假设你研究一个城市2015年有效像元1000个2020年有效像元1200个因为城市边缘原本无光的区域被新光源点亮纳入了有效像元。这时候均值可能被新纳入的低值像元拉低但总和是上升的。如果你只看均值会得出“夜光变暗”的错误结论。我一般会同时输出总和、均值和有效像元数三个指标做趋势时优先看总和。如果非要用均值就固定一个像元集比如只取2015年有效像元在2020年对应的值忽略新增像元。这个操作在zonal_stats里不好直接做需要先用rioxarray做掩膜。4.3 现象重采样后点状光源消失原因双线性插值把峰值摊平解决点源分析用最近邻区域分析用双线性有一次帮人做单个港口的夜光变化重采样到500m后港口的夜光峰值从3000降到了800差点得出“港口衰落”的结论。后来发现是双线性插值把3×3像元的峰值平均了。夜光数据里点状光源港口、机场、单个工业区的峰值很高但空间范围小双线性会把它摊到周围像元导致峰值降低、光晕扩大。如果你研究的是点源重采样一定用最近邻如果研究的是区域总量双线性更平滑不影响总和。这个选择没有绝对对错取决于你的分析单元。4.4 现象负值截断为零后对数变换仍然报错原因零值取对数无定义解决加一个小常数或改用平方根变换夜光数据做回归时常有人用对数变换来压缩动态范围。但预处理后的数据即使截断了负值仍然有大量零值。np.log(0)是负无穷直接报错。常见做法是加1再取对数即log(x1)这样零值变成零不影响单调性。但要注意如果数据里有很多小于1的值加1会改变相对关系。另一种做法是平方根变换对零值友好且比对数更保序。我一般先看零值比例如果低于20%用log(x1)如果高于20%用平方根避免加1带来的偏差。4.5 现象不同年份的栅格范围不一致原因年度合成时裁剪范围不同解决统一用基准年份的extent做裁剪预处理后的数据可能因为每年处理的瓦片数量不同导致最终栅格的范围有差异。比如2012年只覆盖了东部2020年覆盖了全国。如果你直接做年际差值程序会因为形状不一致报错。解决办法是选一个覆盖最全的年份作为基准用rioxarray的rio.clip_box把其他年份裁到相同范围。注意裁剪后要检查nodata比例如果某年裁剪后nodata超过50%说明该年数据本身就不全不适合参与比较。5. 进阶用法把夜光面板变成可解释指标的两个技巧5.1 用夜光总和做GDP空间化的滑动窗口校准夜光总和和GDP的关系不是线性的而且不同区域的斜率差异很大。我一般不用全国统一模型而是用滑动窗口以每个城市为中心取周围一定半径内的城市用这些城市的夜光总和和GDP做线性回归得到局部斜率再把这个斜率应用到中心城市的夜光上。这样做的好处是能捕捉区域经济结构差异比如东部沿海的夜光效率高西部低。具体操作上先用geopandas计算城市间距离矩阵对每个城市取最近的20个邻居用scipy.stats.linregress做回归把斜率和截距存下来。注意窗口大小要试20个邻居在东部可能半径很小在西部可能半径很大我一般会设一个距离上限比如500km超过就不纳入。5.2 用年度夜光变化率识别恢复期以灾后重建为例灾后恢复评估是夜光数据的经典应用。假设某区域在2015年发生灾害夜光骤降之后逐年恢复。你可以计算每个像元的年度变化率然后看恢复曲线的形状。我一般会做三件事第一用2013–2014年的均值作为基准计算2015年的下降幅度第二计算2016–2020年相对于基准的恢复比例第三用逻辑斯蒂曲线拟合恢复过程提取恢复速率。逻辑斯蒂拟合用scipy.optimize.curve_fit参数包括恢复上限、速率和拐点。注意夜光恢复不等于经济恢复有些区域夜光恢复快是因为临时照明不是永久重建。所以我会叠加土地利用数据看恢复的夜光是否落在建成区。import numpy as np from scipy.optimize import curve_fit # 假设某像元2013-2020年的夜光值 years np.arange(2013, 2021) nightlight np.array([120, 125, 40, 60, 90, 110, 118, 122]) # 逻辑斯蒂恢复模型L / (1 exp(-k*(t - t0))) baseline def logistic(t, L, k, t0, baseline): return L / (1 np.exp(-k * (t - t0))) baseline # 初始猜测上限120速率0.8拐点2017基准40 p0 [80, 0.8, 2017, 40] popt, pcov curve_fit(logistic, years, nightlight, p0p0, maxfev5000) # popt里依次是L, k, t0, baseline print(f恢复上限: {popt[0]:.1f}, 速率: {popt[1]:.2f}, 拐点: {popt[2]:.1f}, 基准: {popt[3]:.1f})逻辑说明curve_fit对初始值敏感p0里的拐点最好设在夜光开始明显回升的年份。maxfev调大是为了防止迭代次数不够。拟合出来的k越大恢复越快t0是恢复过半的年份。如果pcov对角线很大说明拟合不稳定可能是数据点太少或恢复曲线不明显。这个模型适合单像元或区域均值不适合直接套到整张栅格上因为每个像元的恢复形态不同。5.3 一个我常犯的错误忽略数据说明里的“建议使用范围”最后说个血泪经验。预处理后的数据通常会在说明里写“建议用于1km及以上尺度的分析”但很多人拿它做500m甚至更细的网格。夜光数据的原始分辨率是750m左右重采样到更细网格不会增加信息量只会引入插值误差。我早期做过一个500m网格的城市内部夜光差异分析结论是“市中心和郊区的夜光梯度很陡”后来发现那个梯度有一部分是重采样造成的。所以拿到数据先看说明里的建议尺度别跟分辨率较劲。希望帮到你。本文还有配套的精品资源点击获取