简介本资源为2012—2020年NPP/VIIRS夜间灯光数据集面向城市遥感、区域经济与地理信息分析方向的研究人员及高年级学生用于解决长时间序列夜间灯光影像难以直接获取与使用的问题。原始影像已完成年度合成、去噪与连续性校正可直接投入城市建成区提取、GDP空间化、人口分布模拟等社会经济指标的空间化分析。压缩包共38个文件约130.74MB以9个tif栅格影像为主体配套9个tfw坐标文件、9个ovr金字塔文件及11个xml元数据兼顾影像读取、坐标定位与快速显示覆盖2012至2020年逐年数据。目前已有1999人学习下载说明该数据在相关研究中具备一定认可度。对于需要构建长时序灯光面板、开展建成区扩张或经济空间格局分析的读者这份经过预处理的成品数据可省去大量清洗与校正环节直接进入建模与制图阶段。1. 2012-2020年NPP/VIIRS夜间灯光数据集从下载到能跑模型中间隔着多少坑如果你拿到的是一份 2012-2020 年 NPP/VIIRS 夜间灯光数据集压缩包解压后大概率会看到一堆按年份或月份命名的 TIF 文件单个文件几百 MB 到几个 GB 不等。很多人第一反应是直接丢进 GIS 软件出图结果发现不同年份的像元值对不上、边缘有异常负值、月度合成里混着杂散光。NPP/VIIRS 夜间灯光数据集的核心价值在于它比 DMSP/OLS 有更高的辐射分辨率和更细的空间尺度但代价是原始产品分了好几个版本每个版本的处理链不一样。这篇内容面向需要把这份数据真正用起来的人——做城市扩张分析、GDP 空间化、碳排放估算或者电力消费建模的从业者。我会按“数据是什么→怎么预处理→怎么验证→坑在哪”的顺序把 2012-2020 年这个时间跨度里最容易翻车的环节讲清楚。2. 先搞清楚你手里的是哪个版本NPP/VIIRS 夜间灯光数据集的三种常见形态2.1 月度合成、年度合成和掩膜版用错一个后面全白做NPP/VIIRS 夜间灯光数据从 2012 年开始由 Suomi NPP 卫星搭载的 VIIRS 传感器采集常见的分发形态有三种。第一种是月度无云合成产品通常命名为VNL_v2_npp_YYYYMM_global_vcmcfg_c2022xxxxx.tif这类格式它做了云掩膜和杂散光校正但保留了月度内的光照变化。第二种是年度合成产品把 12 个月的月度数据做平均或中值合成文件数量少适合做长时间序列趋势分析。第三种是经过进一步掩膜的版本去掉了极光、火点、渔船灯光等非城市光源像元值更“干净”但会损失一部分真实信号。我一般会先看文件名里的vcmcfg还是vcmslcfg。前者是严格云掩膜配置后者是允许更多观测的配置后者在冬季高纬度地区覆盖更好但杂散光残留风险更高。如果你要做 2012-2020 年的连续分析必须确认所有年份用的是同一种配置否则 2012 年和 2020 年的像元值差异里会混进处理链变化带来的偏差。2.2 像元值单位不是辐射亮度直接当 DN 用会出大问题很多人拿到 TIF 后直接读数组发现值域大概在 0 到几百之间就当成 DMSP/OLS 那种 DN 值用了。NPP/VIIRS 月度产品的像元值单位是nW/cm²/sr是辐射亮度不是无量纲的 DN。这意味着两件事第一不同月份的绝对值可以直接比较但要做跨传感器融合时得先做辐射定标第二负值在理论上不应该出现但实际数据里边缘区域会有 -0.5 到 -1.5 左右的负值这是杂散光校正的残留。处理负值的常见做法有两种直接截断为 0或者用邻域中值替换。我倾向于先统计负值像元占比如果小于 0.1%直接截断对整体分析影响可以忽略如果超过 1%说明这个月份的杂散光校正有问题建议换用年度合成产品或者做局部掩膜。2.3 2012 年和 2013 年的数据为什么总被单独讨论2012 年 4 月到 2013 年期间VIIRS 的杂散光校正算法还在迭代导致部分月份在高纬度夏季出现明显的条带噪声。如果你做的是全球尺度分析这两年数据可以用但要做局部城市提取时建议把 2012-2013 年单独做一次阈值标定。我通常会把 2014 年作为基准年因为从这一年开始校正算法稳定后续年份的像元值分布一致性更好。3. 用 Python 把月度 TIF 拼成年度序列从读取到重投影的完整链路3.1 读取单个 TIF 并检查元数据别急着做统计拿到数据后第一步不是算均值而是把元数据看清楚。用rasterio打开一个文件看 CRS、transform、nodata 值和像元尺寸。NPP/VIIRS 全球产品通常是地理坐标系 WGS84像元大小 15 弧秒大约 500 米分辨率。如果你要做区域分析重投影到投影坐标系是必须的但重投影会引入重采样误差所以要在统计之前想清楚。import rasterio import numpy as np # 打开单个月度文件 with rasterio.open(VNL_v2_npp_202001_global_vcmcfg.tif) as src: print(CRS:, src.crs) print(Transform:, src.transform) print(Shape:, src.shape) print(Nodata:, src.nodata) # 读取第一波段 data src.read(1) # 统计负值和零值占比 neg_ratio np.sum(data 0) / data.size zero_ratio np.sum(data 0) / data.size print(f负值占比: {neg_ratio:.4f}, 零值占比: {zero_ratio:.4f}) # 有效值范围 valid data[data 0] print(f有效值范围: {valid.min():.2f} ~ {valid.max():.2f})这段代码的作用是快速判断一个文件是否“干净”。负值占比超过 0.5% 就要警惕零值占比过高可能是云掩膜过度。src.nodata有时候是 None这时候零值既可能是无数据也可能是真实无灯光需要结合掩膜文件判断。3.2 月度合成年度用中值还是均值取决于你的应用场景把 12 个月的月度数据合成年度序列时均值和中值的选择会影响城市边缘区域的提取结果。均值对夏季高值敏感会让城市核心区偏亮中值对异常月份更稳健但会低估季节性灯光变化明显的区域。我一般做城市扩张分析时用中值做电力消费估算时用均值因为后者更接近全年平均辐射水平。import glob import rasterio import numpy as np def monthly_to_annual(year, input_dir, output_path): # 找到该年份所有月度文件 files sorted(glob.glob(f{input_dir}/VNL_v2_npp_{year}*_global_vcmcfg.tif)) if len(files) ! 12: print(f警告: {year}年只有{len(files)}个月度文件) # 读取所有月份数据 stack [] for f in files: with rasterio.open(f) as src: data src.read(1) # 负值截断为0 data[data 0] 0 stack.append(data) profile src.profile # 按像元计算中值 stack np.array(stack) annual_median np.median(stack, axis0) # 写出结果 profile.update(dtyperasterio.float32, count1, compresslzw) with rasterio.open(output_path, w, **profile) as dst: dst.write(annual_median.astype(np.float32), 1) print(f{year}年合成完成输出: {output_path}) # 批量处理2012-2020年 for y in range(2012, 2021): monthly_to_annual(y, ./monthly_data, f./annual/VNL_annual_{y}.tif)这里有几个参数需要根据实际情况调整。data[data 0] 0是简单截断如果你希望保留负值信息用于后续校正可以跳过这一步但在合成前要记录负值比例。np.median在内存里会生成一个三维数组如果处理全球数据12 个月 × 全球范围大约需要 20-30 GB 内存建议分块处理或者用dask做延迟计算。3.3 重投影到区域坐标系别用默认的最近邻插值做区域分析时把全球地理坐标系重投影到 Albers 等面积投影或者 UTM 是常见操作。rasterio的warp默认用最近邻插值这对分类数据没问题但对连续辐射值会引入块状伪影。我一般用双线性插值重采样后的像元值更平滑但会轻微改变极值。from rasterio.warp import calculate_default_transform, reproject, Resampling def reproject_to_albers(src_path, dst_path, dst_crsEPSG:5070): with rasterio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height, dtype: rasterio.float32 }) with rasterio.open(dst_path, w, **kwargs) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear # 关键参数 ) print(f重投影完成: {dst_path})EPSG:5070是北美 Albers 等面积投影如果你做中国区域可以用EPSG:4490或者自定义 Albers。Resampling.bilinear适合连续数据Resampling.cubic更平滑但计算慢Resampling.average适合降尺度。重投影后像元大小会变统计面积时要按新的 transform 重新计算。4. 数据质量验证怎么判断你处理完的 NPP/VIIRS 夜间灯光数据集能不能用4.1 用城市核心区灯光总量做年际一致性检查处理完 2012-2020 年的年度序列后不要直接跑模型。先选一个已知城市核心区比如北京五环内或者上海外环内统计每年的灯光总量和均值。如果某一年突然下降 20% 以上要么是那年的月度数据缺失严重要么是杂散光校正出了问题。我一般会画一条时间序列折线图肉眼扫一遍异常年份再回去查原始月度文件。import rasterio import numpy as np import matplotlib.pyplot as plt def check_city_trend(city_bbox, annual_dir): # city_bbox (left, bottom, right, top) 在对应CRS下 years range(2012, 2021) sums [] for y in years: with rasterio.open(f{annual_dir}/VNL_annual_{y}.tif) as src: # 按窗口读取城市区域 window rasterio.windows.from_bounds(*city_bbox, src.transform) data src.read(1, windowwindow) data[data 0] 0 sums.append(np.sum(data)) plt.plot(years, sums, markero) plt.xlabel(年份) plt.ylabel(灯光总量 (nW/cm²/sr)) plt.title(城市核心区灯光总量年际变化) plt.grid(True) plt.savefig(city_trend.png, dpi150) # 计算年际变化率 for i in range(1, len(sums)): change (sums[i] - sums[i-1]) / sums[i-1] * 100 print(f{years[i]}年变化率: {change:.2f}%) check_city_trend((115.0, 39.5, 117.0, 40.5), ./annual)如果某年变化率超过 ±15%就要标记为可疑年份。注意城市扩张本身会带来灯光增长所以轻微上升是正常的突然下降才是问题。4.2 和 DMSP/OLS 重叠年份做交叉验证2012 和 2013 年有 DMSP/OLS 和 NPP/VIIRS 的重叠观测可以用这两年的数据做交叉验证。把 DMSP/OLS 的 DN 值和 NPP/VIIRS 的辐射值做回归如果 R² 低于 0.6说明两者的空间分布差异较大后续做长时间序列融合时要谨慎。我一般会选 10 个以上城市样本分别提取核心区均值做散点图看趋势。4.3 检查月度文件数量是否完整2012-2020 年一共 108 个月如果你的月度文件夹里少于 100 个文件年度合成的可靠性会下降。缺失月份超过 2 个的年份建议用相邻年份插值或者直接标记为低质量年份。我见过有人直接用 10 个月的数据合成年度结果城市边缘区域出现明显条带这就是月度覆盖不均导致的。5. 避坑与排查NPP/VIIRS 夜间灯光数据预处理里最容易翻车的 5 个地方5.1 现象年度合成后城市核心区出现大面积零值原因月度文件里的零值被当成有效值参与了中值计算而零值在月度数据里既可能是无灯光也可能是云掩膜残留。如果某个月份云掩膜过度该月城市区域大量像元为 0中值合成后就会把真实灯光抹掉。解决在合成前把零值替换为 NaN然后用np.nanmedian计算。同时检查每个月的零值占比超过 30% 的月份直接剔除。data src.read(1).astype(np.float32) data[data 0] np.nan # 零值和负值都设为NaN # 合成时用nanmedian annual np.nanmedian(stack, axis0)5.2 现象重投影后区域面积对不上统计结果偏大或偏小原因重投影时像元大小变了但统计时还在用原始像元面积乘像元数量。地理坐标系下 15 弧秒的像元面积随纬度变化在高纬度地区实际面积比赤道小很多。解决重投影后按新的 transform 计算每个像元的实际面积或者直接用等面积投影。统计总量时用像元值 × 像元面积再求和不要用像元值求和 × 固定面积。5.3 现象2012 年和 2013 年数据和其他年份拼接后趋势线断裂原因这两年部分月份杂散光校正不完善高纬度夏季出现异常高值导致年度均值偏高。如果直接和 2014 年以后的数据拼接趋势线会在 2013-2014 之间出现台阶。解决把 2012-2013 年单独做一次阈值标定或者用 2014 年的城市灯光分布做掩膜只保留稳定灯光区域做趋势分析。我一般会在论文或报告里明确标注这两年数据的处理方式。5.4 现象用rasterio读取大文件时内存溢出原因全球范围的浮点型 TIF 单个文件解压后可能超过 10 GB直接src.read(1)会把整个数组加载到内存。解决用窗口分块读取或者用rasterio的out_shape参数做降采样读取。如果要做全量计算用dask.array配合rioxarray做延迟计算。import rioxarray as rxr import dask.array as da # 用dask分块读取 data rxr.open_rasterio(VNL_annual_2020.tif, chunks{x: 2000, y: 2000}) # 后续计算自动分块 mean_value data.mean().compute()5.5 现象不同年份的 TIF 文件 CRS 或 transform 不一致原因NPP/VIIRS 产品在 2012-2020 年间经历过几次重处理不同批次的文件可能用了略微不同的投影参数或像元对齐方式。解决在批量处理前先写一个脚本检查所有文件的 CRS 和 transform把不一致的文件列出来。如果只是微小偏移可以用rasterio的align功能对齐如果 CRS 不同必须先重投影到统一坐标系再合成。import glob import rasterio files glob.glob(./annual/*.tif) ref_crs None ref_transform None for f in sorted(files): with rasterio.open(f) as src: if ref_crs is None: ref_crs src.crs ref_transform src.transform else: if src.crs ! ref_crs: print(fCRS不一致: {f}) if src.transform ! ref_transform: print(fTransform不一致: {f})6. 进阶技巧用 NPP/VIIRS 夜间灯光数据集做城市扩张分析的阈值选择方法6.1 为什么固定阈值法在 2012-2020 年序列里会失效很多人用 DMSP/OLS 时代的经验设一个固定 DN 阈值比如 10 来提取城市建成区。但 NPP/VIIRS 的辐射值单位不同而且 2012-2020 年间传感器衰变和校正算法变化会导致同一城市的灯光值整体漂移。固定阈值在 2012 年可能提取出 1000 平方公里到 2020 年可能变成 1500 平方公里其中一部分是真实扩张一部分是数据漂移。我一般用两种方法做交叉验证。第一种是相对阈值法先统计整个研究区的灯光值分布取第 95 百分位作为城市核心区阈值第 80 百分位作为城市边缘阈值。第二种是突变检测法对每个像元做时间序列分析找到灯光值突然上升的年份作为城市扩张的起始年。import numpy as np import rasterio def extract_urban_area(tif_path, percentile95): with rasterio.open(tif_path) as src: data src.read(1) data[data 0] 0 # 只统计有灯光的像元 valid data[data 0] if len(valid) 0: return None threshold np.percentile(valid, percentile) urban_mask data threshold # 计算城市面积假设像元面积已知 pixel_area_km2 (500 * 500) / 1e6 # 500米分辨率 area np.sum(urban_mask) * pixel_area_km2 print(f阈值: {threshold:.2f}, 城市面积: {area:.2f} km²) return urban_mask, threshold # 对2012和2020年分别提取 mask_2012, th_2012 extract_urban_area(./annual/VNL_annual_2012.tif, 95) mask_2020, th_2020 extract_urban_area(./annual/VNL_annual_2020.tif, 95)6.2 用灯光值突变点做城市扩张时间定位相对阈值法能告诉你城市范围变了多少但不知道具体哪一年开始扩张。我通常会对每个像元做 Mann-Kendall 趋势检验或者简单的滑动窗口突变检测。具体做法是对每个像元的 2012-2020 年时间序列计算相邻年份的差值如果连续两年差值超过该像元历史差值的 2 倍标准差就标记为突变年。def detect_change_year(city_bbox, annual_dir): years list(range(2012, 2021)) stack [] for y in years: with rasterio.open(f{annual_dir}/VNL_annual_{y}.tif) as src: window rasterio.windows.from_bounds(*city_bbox, src.transform) data src.read(1, windowwindow) data[data 0] 0 stack.append(data) stack np.array(stack) # shape: (9, H, W) # 计算每个像元的年际差值 diff np.diff(stack, axis0) # 计算每个像元差值的标准差 std_diff np.std(diff, axis0) # 找到突变年 change_year np.full(stack.shape[1:], -1, dtypenp.int8) for i in range(diff.shape[0]): mask diff[i] 2 * std_diff change_year[mask (change_year -1)] years[i1] # 统计突变年份分布 unique, counts np.unique(change_year[change_year 0], return_countsTrue) for u, c in zip(unique, counts): print(f{u}年突变像元数: {c}) return change_year detect_change_year((115.0, 39.5, 117.0, 40.5), ./annual)这个方法的参数2 * std_diff可以根据研究区调整城市边缘区域建议用 1.5 倍核心区用 2.5 倍。突变检测的结果可以和 Landsat 影像做交叉验证看看突变年份是否对应实际的新区建设或道路开通。6.3 一个我踩过的坑忽略传感器衰变导致趋势被高估VIIRS 传感器从 2012 年到 2020 年有轻微衰变虽然官方产品做了辐射定标但残余衰变仍然存在。我早期做城市扩张分析时发现所有城市的灯光总量都在上升后来用稳定沙漠区域的灯光值做参考发现 2020 年比 2012 年整体高了约 3-5%。这个偏差在城市扩张速率计算里会放大尤其是扩张缓慢的城市。修正方法很简单选一个 2012-2020 年没有明显人类活动的区域比如沙漠或深海统计其灯光均值年际变化把这个变化率作为传感器衰变参考从所有像元值里扣除。我一般会选撒哈拉沙漠中心或者太平洋偏远海域这些区域在 NPP/VIIRS 数据里灯光值接近零但仍有微小波动。做这个方向的研究数据预处理花的时间往往比建模还多。我自己的习惯是拿到任何年份的 NPP/VIIRS 数据先跑一遍负值统计、零值占比、月度完整性检查三个指标都过了再进入合成流程。2012-2020 年这个时间跨度里2012 和 2013 年要单独对待2014 年以后相对稳定。如果你要做长时间序列分析建议把每年的处理日志留下来包括负值比例、零值比例、异常月份标记后面写论文或报告时这些记录能帮你省很多回头查的时间。希望帮到你。本文还有配套的精品资源点击获取