500米分辨率夜间灯光数据实战用Python处理2000-2023年全国NPP-VIIRS数据附完整代码夜间灯光数据已成为研究城市化进程、经济活动分布和能源消耗的重要指标。对于GIS开发者和数据分析师而言掌握从原始数据到可视化分析的全流程处理技术至关重要。本文将带你用Python完整实现NPP-VIIRS夜间灯光数据的处理涵盖数据读取、跨年校准、水域掩膜处理等核心环节。1. 环境准备与数据加载处理夜间灯光数据需要配置专门的Python环境。推荐使用conda创建独立环境conda create -n nightlight python3.9 conda activate nightlight pip install rasterio geopandas numpy matplotlib关键库说明rasterio专业的栅格数据处理库geopandas处理地理空间矢量数据numpy数值计算核心库matplotlib数据可视化加载.tif格式的夜间灯光数据import rasterio def load_nightlight_data(filepath): with rasterio.open(filepath) as src: data src.read(1) profile src.profile bounds src.bounds return data, profile, bounds # 示例加载2015年数据 data_2015, profile_2015, bounds_2015 load_nightlight_data(NPP_2015.tif)注意原始数据通常采用WGS84坐标系(GCS_WGS_1984)处理时需保持坐标系一致2. 数据预处理与质量控制2.1 异常值处理夜间灯光数据常见问题包括负值传感器异常异常高值火光等临时光源缺失值云层覆盖等处理代码示例import numpy as np def clean_nightlight_data(data): # 处理负值 data[data 0] 0 # 去除极端高值99.9百分位以上 threshold np.percentile(data[data 0], 99.9) data[data threshold] threshold # 填充缺失值通常用邻近像元均值 from scipy.ndimage import generic_filter def fill_nan(window): window window[~np.isnan(window)] return window.mean() if len(window) 0 else np.nan if np.isnan(data).any(): data generic_filter(data, fill_nan, size3) return data2.2 跨传感器校准DMSP-OLS与NPP-VIIRS数据校准是关键挑战。基于Chen等(2021)的方法实现校准系数计算def cross_sensor_calibration(dmsp_data, viirs_data, maskNone): dmsp_data: DMSP-OLS数据(2000-2012) viirs_data: NPP-VIIRS数据(2013-2018) mask: 可选的水域掩膜 if mask is not None: dmsp_masked dmsp_data[mask] viirs_masked viirs_data[mask] else: dmsp_masked dmsp_data viirs_masked viirs_data # 计算校准系数 valid_pixels (dmsp_masked 0) (viirs_masked 0) x dmsp_masked[valid_pixels] y viirs_masked[valid_pixels] # 稳健回归避免异常值影响 from sklearn.linear_model import TheilSenRegressor model TheilSenRegressor().fit(x.reshape(-1,1), y) slope model.coef_[0] intercept model.intercept_ return slope, intercept3. 水域掩膜处理技术水域区域灯光数据通常需要特殊处理。以下是创建水域掩膜的完整流程3.1 获取水域数据推荐使用全球地表水数据集GLWD或HydroSHEDSimport geopandas as gpd def load_water_mask(water_shapefile, target_crs, target_bounds): water_shapefile: 水域矢量数据路径 target_crs: 目标坐标系 target_bounds: 目标范围 water gpd.read_file(water_shapefile).to_crs(target_crs) # 创建栅格掩膜 from rasterio.features import geometry_mask mask geometry_mask( water.geometry, out_shape(int((target_bounds[3]-target_bounds[1])/0.00416667), # 500m分辨率 int((target_bounds[2]-target_bounds[0])/0.00416667)), transformrasterio.transform.from_bounds(*target_bounds, widthint((target_bounds[2]-target_bounds[0])/0.00416667), heightint((target_bounds[3]-target_bounds[1])/0.00416667)), invertTrue ) return mask3.2 应用水域掩膜def apply_water_mask(data, water_mask, fill_value0): data: 夜间灯光数据 water_mask: 水域掩膜(True表示水域) fill_value: 水域填充值 masked_data data.copy() masked_data[water_mask] fill_value return masked_data4. 时间序列分析与可视化4.1 多年数据整合构建2000-2023年时间序列数据集import os import pandas as pd def build_time_series(data_dir, water_maskNone): years range(2000, 2024) time_series [] for year in years: file os.path.join(data_dir, fNPP_{year}.tif) data, _, _ load_nightlight_data(file) if water_mask is not None: data apply_water_mask(data, water_mask) # 计算年度总灯光强度 total_light data[data 0].sum() # 计算灯光区域平均强度 mean_light data[data 0].mean() # 计算灯光区域面积(km²) light_area (data 0).sum() * 0.25 # 500m分辨率每像素0.25km² time_series.append({ year: year, total_light: total_light, mean_light: mean_light, light_area: light_area }) return pd.DataFrame(time_series)4.2 可视化分析使用matplotlib创建专业可视化import matplotlib.pyplot as plt def plot_light_trends(df): fig, (ax1, ax2, ax3) plt.subplots(3, 1, figsize(12, 15)) # 总灯光强度趋势 ax1.plot(df[year], df[total_light], o-, color#FF7F0E) ax1.set_title(Total Nighttime Light Intensity (2000-2023)) ax1.set_ylabel(Total Light (nW/cm²/sr)) # 平均灯光强度趋势 ax2.plot(df[year], df[mean_light], s-, color#1F77B4) ax2.set_title(Average Nighttime Light Intensity (2000-2023)) ax2.set_ylabel(Mean Light (nW/cm²/sr)) # 灯光区域面积趋势 ax3.plot(df[year], df[light_area], d-, color#2CA02C) ax3.set_title(Lighted Area (2000-2023)) ax3.set_ylabel(Area (km²)) ax3.set_xlabel(Year) for ax in [ax1, ax2, ax3]: ax.grid(True, linestyle--, alpha0.6) plt.tight_layout() plt.savefig(light_trends.png, dpi300, bbox_inchestight) plt.show()5. 空间分析与热点检测5.1 灯光变化热点分析使用空间统计方法检测显著变化区域from pysal.explore import esda from pysal.lib import weights def detect_hotspots(data_early, data_late, threshold0.05): 检测灯光强度变化显著的热点区域 data_early: 早期数据(如2000年) data_late: 后期数据(如2023年) threshold: 显著性水平 # 计算变化率 change (data_late - data_early) / (data_early 1e-10) # 避免除零 # 创建空间权重矩阵(皇后邻接) w weights.Queen.from_array(np.ones_like(change)) # 计算Getis-Ord Gi*统计量 gi esda.G_Local(change, w, transformB) # 创建热点掩膜(显著性水平0.05) hotspots (gi.z_sim 1.96) (change 0) coldspots (gi.z_sim -1.96) (change 0) return hotspots, coldspots5.2 空间可视化def plot_spatial_changes(data_early, data_late, hotspots, coldspots, save_path): fig, (ax1, ax2, ax3) plt.subplots(1, 3, figsize(18, 6)) # 早期灯光 im1 ax1.imshow(data_early, cmaphot, vmaxnp.percentile(data_early, 99)) ax1.set_title(Early Period (2000)) fig.colorbar(im1, axax1, labelLight Intensity) # 后期灯光 im2 ax2.imshow(data_late, cmaphot, vmaxnp.percentile(data_late, 99)) ax2.set_title(Late Period (2023)) fig.colorbar(im2, axax2, labelLight Intensity) # 变化热点 from matplotlib.colors import ListedColormap cmap ListedColormap([blue, grey, red]) bounds [-1, -0.5, 0.5, 1] norm plt.Normalize(-1, 1) change_map np.zeros_like(hotspots, dtypeint) change_map[coldspots] -1 change_map[hotspots] 1 im3 ax3.imshow(change_map, cmapcmap, normnorm) ax3.set_title(Change Hotspots (p0.05)) cbar fig.colorbar(im3, axax3, ticks[-0.66, 0, 0.66]) cbar.ax.set_yticklabels([Coldspots, No change, Hotspots]) plt.tight_layout() plt.savefig(save_path, dpi300, bbox_inchestight) plt.show()6. 实际应用案例城市化进程评估结合夜间灯光数据与统计数据进行城市发展评估import statsmodels.api as sm def urban_growth_analysis(light_df, stats_df): light_df: 灯光数据DataFrame(含year,total_light等列) stats_df: 统计数据DataFrame(含year,population,gdp等列) merged pd.merge(light_df, stats_df, onyear) # 标准化变量 for col in [total_light, population, gdp]: merged[col_norm] (merged[col] - merged[col].min()) / (merged[col].max() - merged[col].min()) # 建立回归模型 X merged[[total_light_norm]] y_pop merged[population_norm] y_gdp merged[gdp_norm] model_pop sm.OLS(y_pop, sm.add_constant(X)).fit() model_gdp sm.OLS(y_gdp, sm.add_constant(X)).fit() # 绘制关系图 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) ax1.scatter(X, y_pop, colorblue) ax1.plot(X, model_pop.predict(sm.add_constant(X)), colorred) ax1.set_title(Nighttime Light vs Population) ax1.set_xlabel(Normalized Light Intensity) ax1.set_ylabel(Normalized Population) ax2.scatter(X, y_gdp, colorgreen) ax2.plot(X, model_gdp.predict(sm.add_constant(X)), colorred) ax2.set_title(Nighttime Light vs GDP) ax2.set_xlabel(Normalized Light Intensity) ax2.set_ylabel(Normalized GDP) plt.tight_layout() plt.show() return model_pop, model_gdp夜间灯光数据的处理流程中最耗时的环节通常是大数据量的栅格运算。在实际项目中我通常会使用Dask进行并行计算特别是处理全国范围的多年度数据时性能提升非常明显。另一个实用技巧是预先计算并存储中间结果比如水域掩膜和校准系数避免重复计算。
500米分辨率夜间灯光数据实战：用Python处理2000-2023年全国NPP-VIIRS数据（附完整代码）
发布时间：2026/5/19 18:15:44
500米分辨率夜间灯光数据实战用Python处理2000-2023年全国NPP-VIIRS数据附完整代码夜间灯光数据已成为研究城市化进程、经济活动分布和能源消耗的重要指标。对于GIS开发者和数据分析师而言掌握从原始数据到可视化分析的全流程处理技术至关重要。本文将带你用Python完整实现NPP-VIIRS夜间灯光数据的处理涵盖数据读取、跨年校准、水域掩膜处理等核心环节。1. 环境准备与数据加载处理夜间灯光数据需要配置专门的Python环境。推荐使用conda创建独立环境conda create -n nightlight python3.9 conda activate nightlight pip install rasterio geopandas numpy matplotlib关键库说明rasterio专业的栅格数据处理库geopandas处理地理空间矢量数据numpy数值计算核心库matplotlib数据可视化加载.tif格式的夜间灯光数据import rasterio def load_nightlight_data(filepath): with rasterio.open(filepath) as src: data src.read(1) profile src.profile bounds src.bounds return data, profile, bounds # 示例加载2015年数据 data_2015, profile_2015, bounds_2015 load_nightlight_data(NPP_2015.tif)注意原始数据通常采用WGS84坐标系(GCS_WGS_1984)处理时需保持坐标系一致2. 数据预处理与质量控制2.1 异常值处理夜间灯光数据常见问题包括负值传感器异常异常高值火光等临时光源缺失值云层覆盖等处理代码示例import numpy as np def clean_nightlight_data(data): # 处理负值 data[data 0] 0 # 去除极端高值99.9百分位以上 threshold np.percentile(data[data 0], 99.9) data[data threshold] threshold # 填充缺失值通常用邻近像元均值 from scipy.ndimage import generic_filter def fill_nan(window): window window[~np.isnan(window)] return window.mean() if len(window) 0 else np.nan if np.isnan(data).any(): data generic_filter(data, fill_nan, size3) return data2.2 跨传感器校准DMSP-OLS与NPP-VIIRS数据校准是关键挑战。基于Chen等(2021)的方法实现校准系数计算def cross_sensor_calibration(dmsp_data, viirs_data, maskNone): dmsp_data: DMSP-OLS数据(2000-2012) viirs_data: NPP-VIIRS数据(2013-2018) mask: 可选的水域掩膜 if mask is not None: dmsp_masked dmsp_data[mask] viirs_masked viirs_data[mask] else: dmsp_masked dmsp_data viirs_masked viirs_data # 计算校准系数 valid_pixels (dmsp_masked 0) (viirs_masked 0) x dmsp_masked[valid_pixels] y viirs_masked[valid_pixels] # 稳健回归避免异常值影响 from sklearn.linear_model import TheilSenRegressor model TheilSenRegressor().fit(x.reshape(-1,1), y) slope model.coef_[0] intercept model.intercept_ return slope, intercept3. 水域掩膜处理技术水域区域灯光数据通常需要特殊处理。以下是创建水域掩膜的完整流程3.1 获取水域数据推荐使用全球地表水数据集GLWD或HydroSHEDSimport geopandas as gpd def load_water_mask(water_shapefile, target_crs, target_bounds): water_shapefile: 水域矢量数据路径 target_crs: 目标坐标系 target_bounds: 目标范围 water gpd.read_file(water_shapefile).to_crs(target_crs) # 创建栅格掩膜 from rasterio.features import geometry_mask mask geometry_mask( water.geometry, out_shape(int((target_bounds[3]-target_bounds[1])/0.00416667), # 500m分辨率 int((target_bounds[2]-target_bounds[0])/0.00416667)), transformrasterio.transform.from_bounds(*target_bounds, widthint((target_bounds[2]-target_bounds[0])/0.00416667), heightint((target_bounds[3]-target_bounds[1])/0.00416667)), invertTrue ) return mask3.2 应用水域掩膜def apply_water_mask(data, water_mask, fill_value0): data: 夜间灯光数据 water_mask: 水域掩膜(True表示水域) fill_value: 水域填充值 masked_data data.copy() masked_data[water_mask] fill_value return masked_data4. 时间序列分析与可视化4.1 多年数据整合构建2000-2023年时间序列数据集import os import pandas as pd def build_time_series(data_dir, water_maskNone): years range(2000, 2024) time_series [] for year in years: file os.path.join(data_dir, fNPP_{year}.tif) data, _, _ load_nightlight_data(file) if water_mask is not None: data apply_water_mask(data, water_mask) # 计算年度总灯光强度 total_light data[data 0].sum() # 计算灯光区域平均强度 mean_light data[data 0].mean() # 计算灯光区域面积(km²) light_area (data 0).sum() * 0.25 # 500m分辨率每像素0.25km² time_series.append({ year: year, total_light: total_light, mean_light: mean_light, light_area: light_area }) return pd.DataFrame(time_series)4.2 可视化分析使用matplotlib创建专业可视化import matplotlib.pyplot as plt def plot_light_trends(df): fig, (ax1, ax2, ax3) plt.subplots(3, 1, figsize(12, 15)) # 总灯光强度趋势 ax1.plot(df[year], df[total_light], o-, color#FF7F0E) ax1.set_title(Total Nighttime Light Intensity (2000-2023)) ax1.set_ylabel(Total Light (nW/cm²/sr)) # 平均灯光强度趋势 ax2.plot(df[year], df[mean_light], s-, color#1F77B4) ax2.set_title(Average Nighttime Light Intensity (2000-2023)) ax2.set_ylabel(Mean Light (nW/cm²/sr)) # 灯光区域面积趋势 ax3.plot(df[year], df[light_area], d-, color#2CA02C) ax3.set_title(Lighted Area (2000-2023)) ax3.set_ylabel(Area (km²)) ax3.set_xlabel(Year) for ax in [ax1, ax2, ax3]: ax.grid(True, linestyle--, alpha0.6) plt.tight_layout() plt.savefig(light_trends.png, dpi300, bbox_inchestight) plt.show()5. 空间分析与热点检测5.1 灯光变化热点分析使用空间统计方法检测显著变化区域from pysal.explore import esda from pysal.lib import weights def detect_hotspots(data_early, data_late, threshold0.05): 检测灯光强度变化显著的热点区域 data_early: 早期数据(如2000年) data_late: 后期数据(如2023年) threshold: 显著性水平 # 计算变化率 change (data_late - data_early) / (data_early 1e-10) # 避免除零 # 创建空间权重矩阵(皇后邻接) w weights.Queen.from_array(np.ones_like(change)) # 计算Getis-Ord Gi*统计量 gi esda.G_Local(change, w, transformB) # 创建热点掩膜(显著性水平0.05) hotspots (gi.z_sim 1.96) (change 0) coldspots (gi.z_sim -1.96) (change 0) return hotspots, coldspots5.2 空间可视化def plot_spatial_changes(data_early, data_late, hotspots, coldspots, save_path): fig, (ax1, ax2, ax3) plt.subplots(1, 3, figsize(18, 6)) # 早期灯光 im1 ax1.imshow(data_early, cmaphot, vmaxnp.percentile(data_early, 99)) ax1.set_title(Early Period (2000)) fig.colorbar(im1, axax1, labelLight Intensity) # 后期灯光 im2 ax2.imshow(data_late, cmaphot, vmaxnp.percentile(data_late, 99)) ax2.set_title(Late Period (2023)) fig.colorbar(im2, axax2, labelLight Intensity) # 变化热点 from matplotlib.colors import ListedColormap cmap ListedColormap([blue, grey, red]) bounds [-1, -0.5, 0.5, 1] norm plt.Normalize(-1, 1) change_map np.zeros_like(hotspots, dtypeint) change_map[coldspots] -1 change_map[hotspots] 1 im3 ax3.imshow(change_map, cmapcmap, normnorm) ax3.set_title(Change Hotspots (p0.05)) cbar fig.colorbar(im3, axax3, ticks[-0.66, 0, 0.66]) cbar.ax.set_yticklabels([Coldspots, No change, Hotspots]) plt.tight_layout() plt.savefig(save_path, dpi300, bbox_inchestight) plt.show()6. 实际应用案例城市化进程评估结合夜间灯光数据与统计数据进行城市发展评估import statsmodels.api as sm def urban_growth_analysis(light_df, stats_df): light_df: 灯光数据DataFrame(含year,total_light等列) stats_df: 统计数据DataFrame(含year,population,gdp等列) merged pd.merge(light_df, stats_df, onyear) # 标准化变量 for col in [total_light, population, gdp]: merged[col_norm] (merged[col] - merged[col].min()) / (merged[col].max() - merged[col].min()) # 建立回归模型 X merged[[total_light_norm]] y_pop merged[population_norm] y_gdp merged[gdp_norm] model_pop sm.OLS(y_pop, sm.add_constant(X)).fit() model_gdp sm.OLS(y_gdp, sm.add_constant(X)).fit() # 绘制关系图 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) ax1.scatter(X, y_pop, colorblue) ax1.plot(X, model_pop.predict(sm.add_constant(X)), colorred) ax1.set_title(Nighttime Light vs Population) ax1.set_xlabel(Normalized Light Intensity) ax1.set_ylabel(Normalized Population) ax2.scatter(X, y_gdp, colorgreen) ax2.plot(X, model_gdp.predict(sm.add_constant(X)), colorred) ax2.set_title(Nighttime Light vs GDP) ax2.set_xlabel(Normalized Light Intensity) ax2.set_ylabel(Normalized GDP) plt.tight_layout() plt.show() return model_pop, model_gdp夜间灯光数据的处理流程中最耗时的环节通常是大数据量的栅格运算。在实际项目中我通常会使用Dask进行并行计算特别是处理全国范围的多年度数据时性能提升非常明显。另一个实用技巧是预先计算并存储中间结果比如水域掩膜和校准系数避免重复计算。
