简介这份资源面向从事城市扩张、社会经济估算、生态环境监测等研究的遥感与GIS人员提供2012—2020年NPP-VIIRS夜间灯光数据分辨率500米已完成基础预处理以WGS84中国矢量边界裁剪出全国范围VNL2数据并转换为Asia_Lambert_Conformal_Conic平面投影坐标系同时参照Kang Wu等2019以472.86为全国像元辐射阈值剔除极端异常值并去除负值像元可直接用于时序分析与建模。压缩包共36个文件约176.73MB以9个tif栅格为主体配套tfw坐标信息、xml元数据与ovr金字塔文件便于在ArcGIS等软件中快速加载与配准。目前已有3508人学习下载适合需要现成、规范夜光数据支撑论文或项目分析的研究者使用。1. 拿到一份“预处理后”的夜光遥感数据先别急着画图夜光遥感数据在从业者手里最常见的用法是拿来做区域经济活动强度、城镇化扩张、人口分布代理指标的估算。但真正下过手的人都知道原始 NPP-VIIRS 月度产品里塞满了杂散光、极光、火点、油气燃烧、背景噪声直接拿来算均值出来的曲线能把你带沟里去。这份 2012–2020 年预处理后的夜光遥感数据集价值就在“预处理”三个字上——它把跨传感器校正、异常值剔除、负值与背景噪声处理这些脏活提前做完了你拿到的是可以直接进入统计和建模环节的年际序列。适合谁做区域经济、城市研究、资源环境遥感的研究生和一线分析人员尤其是那些不想在数据清洗上再耗两周、只想快速验证假设的人。下面按“这份数据是什么结构 → 怎么读进来 → 怎么算指标 → 坑在哪 → 怎么验证”的顺序拆。2. 数据组织与读取先搞清楚 DN 值和辐射定标的关系2.1 预处理到底动了哪些手脚NPP-VIIRS 的夜间灯光产品VCMSLCFG 系列原始 DN 值并不是物理辐射量而是经过缩放处理的整数。常见的预处理链条包括去除受杂散光影响的像元、剔除闪电与火点造成的瞬时高值、对负值做归零或掩膜、跨年辐射定标一致性调整。这份 2012–2020 数据既然标注“预处理后”意味着上述步骤已经完成你拿到的 DN 值在年际间具备可比性——这一点非常关键因为原始月度产品在不同年份的定标系数不同直接跨年比较会引入系统性偏差。需要提醒的是“预处理后”不等于“已投影到统一网格”。不同来源的预处理产品有的保留地理经纬度坐标WGS84有的重采样到等面积投影。读取前先看元数据里的坐标系和像元大小别默认它是某个固定投影。常见做法是用rasterio或gdalinfo先探一下。# 先看栅格元信息坐标系、像元大小、NoData 值、波段数 gdalinfo NPP_VIIRS_2015_preprocessed.tif这条命令输出里重点盯三样Coordinate System是否为地理坐标、Pixel Size的绝对值约 0.00417 度对应约 500m 赤道分辨率、NoData Value。如果 NoData 是 -9999 之类后续统计必须显式排除否则均值会被拉垮。2.2 用 Python 把年度序列读成数组单年读取简单但做趋势分析要的是 2012–2020 的堆叠。我一般用rasterio逐文件读统一到同一网格后再堆成三维数组。下面这段是可直接抄的骨架import rasterio import numpy as np import glob # 按年份排序确保时间序列顺序正确 files sorted(glob.glob(data/NPP_VIIRS_*_preprocessed.tif)) years list(range(2012, 2021)) stack [] for f in files: with rasterio.open(f) as src: arr src.read(1).astype(float32) nodata src.nodata # 关键把 NoData 和负值统一置为 NaN避免污染统计 if nodata is not None: arr[arr nodata] np.nan arr[arr 0] np.nan stack.append(arr) cube np.stack(stack, axis0) # 形状 (9, H, W) print(cube.shape, np.nanmean(cube[0]))逻辑说明sorted保证文件按年份排列否则时间轴错位会让趋势完全反掉astype(float32)是为了后续能存 NaN整型数组塞不进 NaN负值置 NaN 是因为部分预处理产品保留了负 DN直接参与均值会低估亮度。参数上如果你的文件命名不是按年份可排序的务必手动建一个year - filepath的字典别赌文件名。2.3 裁剪到研究区再算别对全球求均值全球均值对区域研究毫无意义。裁剪有两种路子按矢量边界掩膜或按经纬度范围切片。前者更严谨后者更快。用rasterio.mask做矢量裁剪时注意cropTrue会改变数组尺寸后续堆叠要保证所有年份裁出来尺寸一致——如果边界跨了不同分幅先做镶嵌再裁。from rasterio.mask import mask import geopandas as gpd boundary gpd.read_file(study_area.shp) with rasterio.open(files[0]) as src: geom [boundary.to_crs(src.crs).geometry.unary_union] out_image, out_transform mask(src, geom, cropTrue, nodatanp.nan)这里to_crs(src.crs)不能省矢量与栅格坐标系不一致是裁剪失败的头号原因报错往往还很隐晦。nodatanp.nan要求输出为浮点若源数据是整型会报错所以前面统一转 float32 是有必要的。3. 从 DN 值到可用指标均值、总和与灯光面积怎么选3.1 三种常用统计量的适用边界拿到裁剪后的数组下一步是把它压成一个能进回归或时序分析的数字。常见三种区域内 DN 均值、DN 总和、以及阈值以上的灯光像元面积。它们不是随便挑的各有适用场景。指标计算方式适用场景主要风险DN 均值有效像元 DN 求和 / 有效像元数比较亮度强度变化对研究区面积变化敏感DN 总和全部有效像元 DN 求和总量型经济代理受区域面积和像元数影响灯光面积DN 阈值的像元数 × 像元面积城镇化扩张、建成区提取阈值选择主观性强均值适合研究区边界固定的场景比如同一个行政区多年对比总和适合做经济总量代理但行政区划调整的年份会跳变灯光面积对阈值极其敏感阈值差几个 DN面积能差出百分之十几。我的习惯是三个都算互相印证哪个出现异常跳变就回头查那年的数据。3.2 阈值怎么定别拍脑袋用分位数或文献锚点灯光面积的阈值没有统一标准文献里从 1 到 30 都有人用。实操中我一般先看研究区 DN 的分布取某个分位数作为相对阈值或者用已知建成区样本反推。下面这段给出分位数法和固定阈值法两种def lit_area(cube, pixel_area_km2, methodquantile, q0.9, fixed8): areas [] for t in range(cube.shape[0]): band cube[t] valid band[~np.isnan(band)] if method quantile: thr np.quantile(valid, q) # 相对阈值随年份自适应 else: thr fixed # 固定阈值跨年可比 areas.append(np.sum(band thr) * pixel_area_km2) return np.array(areas)逻辑说明分位数法每年阈值随分布漂移适合看“相对亮区”的扩张固定阈值法跨年可比适合做时间序列回归。pixel_area_km2必须按纬度换算500m 分辨率在赤道约 0.25 km²高纬度会缩小粗略估算可用(0.00417 * 111.32 * cos(lat))²。参数q取 0.9 意味着最亮 10% 像元算作灯光区这个值要结合你的研究区亮度分布调城市密集区可以取 0.95。3.3 把年度指标落成一张可回归的表算完指标输出成 CSV 才能进 Stata、R 或 Python 的统计模型。别小看这一步列名和单位不统一后面合并社会经济数据时全是麻烦。import pandas as pd df pd.DataFrame({ year: years, mean_dn: [np.nanmean(cube[i]) for i in range(cube.shape[0])], sum_dn: [np.nansum(cube[i]) for i in range(cube.shape[0])], lit_area_km2: lit_area(cube, pixel_area_km20.25, methodfixed, fixed8) }) df.to_csv(nightlight_indicators_2012_2020.csv, indexFalse, encodingutf-8-sig)encodingutf-8-sig是为了 Excel 打开中文列名不乱码这个细节踩过坑的人都懂。np.nansum与np.nanmean会自动跳过 NaN但前提是你前面确实把 NoData 置成了 NaN否则它们会把 -9999 当有效值加进去。4. 避坑与排查预处理数据也会翻车的地方4.1 现象某年指标突然断崖式下跌 → 原因该年文件波段或定标版本不同 → 解决逐年核对元数据跨年序列最怕某一年的产品版本和其他年份不一致。表现是均值曲线在某个年份突然掉一截看着像经济衰退其实是数据问题。排查方法是逐年打印gdalinfo的波段统计和元数据里的定标参数对比 DN 的最大值量级。如果某年最大值只有其他年份的十分之一基本可以判定定标系数不同需要按元数据里的缩放因子还原或者直接剔除该年。4.2 现象裁剪后数组尺寸不一致堆叠报错 → 原因矢量跨分幅或cropTrue行为差异 → 解决先镶嵌到统一网格再裁mask的cropTrue会按几何外接矩形裁不同年份若源数据分幅边界不同裁出来尺寸就差几个像素。解决办法是先对所有年份做rasterio.merge镶嵌到同一网格或者用reproject统一到目标网格后再裁。别试图用np.resize硬凑那会把地理对应关系彻底搞乱。4.3 现象均值算出来是负数或异常小 → 原因NoData 未识别-9999 参与了运算 → 解决读入即掩膜这是最高频的翻车点。很多预处理产品的 NoData 是 -9999 或 -3.4e38如果只判断 0置 NaN-9999 确实被处理了但有些产品 NoData 是 0 而有效暗区也是 0一刀切会误伤。正确做法是优先用src.nodata读元数据里的 NoData 值精确匹配匹配不到再考虑负值掩膜。算完均值后随手print一下最小值和有效像元占比占比低于 50% 就要警惕。4.4 现象灯光面积年际波动剧烈和常识不符 → 原因阈值用了分位数法阈值随分布漂移 → 解决时间序列分析改用固定阈值分位数法每年重新定阈值等于每年换了一把尺子面积序列反映的是相对变化而非绝对扩张。做趋势回归时必须用固定阈值且阈值要在所有年份上统一。如果研究区跨越不同亮度等级比如既有大城市又有乡村单一固定阈值可能两头不讨好这时可以分区域分别定阈值但要在论文里说清楚。4.5 现象高纬度研究区面积估算偏大 → 原因像元面积按赤道值算 → 解决按纬度逐像元换算面积地理坐标下的像元面积随纬度收缩用固定 0.25 km² 在高纬度会高估。严谨做法是用pyproj或rasterio的变换参数逐行计算像元实际面积再求和。粗略修正可以用研究区中心纬度的余弦因子但跨纬度大的区域误差仍不可忽略。5. 验证与进阶用不变区域和辅助数据给结果上保险数据算完不能直接交差得验证。我常用的两个手段一是找一块理论上灯光变化极小的区域比如大片自然水体或无人区做“不变区域”检验如果它的 DN 均值年际波动超过 5%说明预处理或读取环节有问题二是拿已知的建成区矢量或统计年鉴里的城镇化率做相关性对照夜光指标和城镇化率的相关系数低于 0.7 就要回头查阈值和裁剪。进阶用法上这份 2012–2020 的连续序列很适合做趋势分解。简单线性回归只能看单调趋势但夜光受经济周期影响明显用 Theil-Sen 斜率比最小二乘稳健能抗异常年干扰。下面这段给出逐像元趋势的写法适合做空间上的“哪里在变亮”from scipy.stats import theilslopes h, w cube.shape[1], cube.shape[2] slope np.full((h, w), np.nan, dtypefloat32) for i in range(h): for j in range(w): ts cube[:, i, j] if np.isnan(ts).any(): continue # Theil-Sen 对异常值稳健适合夜光这种有火点残留的序列 s, _, _, _ theilslopes(ts, years) slope[i, j] s逐像元循环在 Python 里慢大区域建议用xarray的polyfit或分块并行。theilslopes返回的斜率单位是 DN/年正值表示变亮。算完把 slope 栅格导出叠加行政区边界就能直观看到哪些区域在扩张、哪些在收缩。这一步做完再回头对照前面算的区域均值两者方向一致才算过关。血泪经验是夜光数据的“预处理后”只是起点不是终点。我现在的习惯是每次拿到新序列先跑一遍不变区域检验和 NoData 占比检查再动统计模型——这两步花不了十分钟但能省掉后面返工两周的后悔药。希望帮到你。本文还有配套的精品资源点击获取