从下载到出图：一份给GIS新手的VIIRS夜光数据保姆级处理指南（附Python代码）

从零开始处理VIIRS夜光数据：Python实战指南

夜光遥感数据正在成为城市发展、经济活动和环境监测的重要指标。作为NASA/NOAA联合推出的新一代夜光数据产品，VIIRS（Visible Infrared Imaging Radiometer Suite）凭借其高分辨率和稳定的数据质量，在学术研究和商业分析中展现出巨大价值。但对于刚接触GIS和遥感数据的开发者来说，从原始数据到可用成果的完整流程往往充满挑战——复杂的文件命名规则、庞大的TIFF文件体积、专业术语的解读障碍，都可能让初学者望而却步。

本文将采用"问题导向"的实战路线，带您完整走通VIIRS夜光数据的处理全流程。不同于传统教程的理论讲解，我们聚焦于实际项目中遇到的典型问题：如何快速定位需要的数据产品？如何高效处理动辄数GB的GeoTIFF文件？如何针对特定区域进行精确裁剪？我们将使用Python生态中的主流工具链（rasterio、GDAL、NumPy等），通过可复用的代码示例解决这些实际问题。无论您是城市规划专业的学生，还是需要夜光数据支撑业务分析的数据工程师，都能从中获得可直接落地的技术方案。

1. 数据准备与环境配置

1.1 理解VIIRS数据产品体系

VIIRS夜光数据主要分为三个核心产品线，各自适用于不同分析场景：

提示：初学者建议从VCM月度数据入手，其稳定的质量控制和相对较小的数据量更适合练手。中国区域的T3级别数据（文件名含75N060E）覆盖全国范围。

1.2 Python环境准备

推荐使用conda创建专属的GIS处理环境，避免库版本冲突：

conda create -n viirs python=3.9 conda activate viirs conda install -c conda-forge rasterio geopandas matplotlib numpy jupyterlab

关键库的作用说明：

• rasterio：高效读写GeoTIFF文件，支持内存映射处理大文件
• geopandas：处理矢量边界数据（如行政区划shp文件）
• matplotlib：数据可视化与出图
• numpy：底层数组运算与数据掩膜处理

验证安装是否成功：

import rasterio print(rasterio.__version__) # 应输出1.3.0以上版本

2. 数据下载与文件解析

2.1 从NOAA官网获取数据

VIIRS夜光数据的官方下载入口是NOAA的Earth Observation Group页面。实际下载时需要注意：

1. 选择年度/月度复合产品（避免原始每日数据量过大）
2. 确认产品类型（vcm、vcm-orm-ntl等）
3. 检查TILE编号（中国全境通常选择75N060E）
4. 注意版本号（v21c为当前稳定版本）

典型文件名示例：VNL_v2_npp_2020_global_vcmslcfg_c202102211500.avg_rade9h.tif.gz

各字段含义解析：

• VNL_v2：VIIRS夜光数据版本2
• npp：卫星平台（Suomi NPP）
• 2020：数据年份
• global：全球覆盖
• vcmslcfg：数据产品类型（VCM月度合成）
• c202102211500：产品生成时间戳
• avg_rade9h：辐射校正后的平均辐射值

2.2 文件解压与预处理

下载的.gz压缩文件需要解压后才能使用。在Linux/macOS下可直接使用gunzip命令，Windows用户推荐7-Zip工具。Python中也可直接处理：

import gzip import shutil with gzip.open('input.tif.gz', 'rb') as f_in: with open('output.tif', 'wb') as f_out: shutil.copyfileobj(f_in, f_out)

解压后的TIFF文件通常大小在1-3GB之间，建议检查文件完整性：

with rasterio.open('output.tif') as src: print(f"宽度: {src.width} 像素") print(f"高度: {src.height} 像素") print(f"波段数: {src.count}") print(f"CRS: {src.crs}") # 应显示EPSG:4326

3. 数据处理核心流程

3.1 数据裁剪与区域提取

实际分析中，我们通常只需要特定区域的数据。以下示例展示如何使用省级行政区划裁剪数据：

import geopandas as gpd from rasterio.mask import mask # 加载行政区划矢量数据 shp_path = "province_boundary.shp" gdf = gpd.read_file(shp_path) geometry = gdf[gdf['NAME']=='广东省'].geometry # 执行裁剪 with rasterio.open('output.tif') as src: out_image, out_transform = mask( src, geometry, crop=True, all_touched=True ) out_meta = src.meta.copy() # 更新元数据并保存 out_meta.update({ "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform }) with rasterio.open('guangdong.tif', 'w', **out_meta) as dest: dest.write(out_image)

注意：all_touched=True确保边界像素不被错误剔除，但会增加计算量。对精确分析建议开启此选项。

3.2 数据可视化技巧

夜光数据的有效可视化需要特殊的配色方案。以下代码创建专业级的夜间灯光图：

import matplotlib.pyplot as plt import numpy as np def plot_nightlight(data, title, save_path=None): # 创建自定义colormap colors = ['black', 'blue', 'cyan', 'yellow', 'white'] cmap = plt.colors.LinearSegmentedColormap.from_list('nightlight', colors) # 对数变换增强低值区可视性 norm = plt.colors.LogNorm( vmin=np.nanmin(data)+1e-9, vmax=np.nanmax(data) ) # 绘制图像 fig, ax = plt.subplots(figsize=(12, 10)) img = ax.imshow(data, cmap=cmap, norm=norm) plt.colorbar(img, ax=ax, label='辐射值 (nW/cm²/sr)') ax.set_title(title, fontsize=14) ax.axis('off') if save_path: plt.savefig(save_path, dpi=300, bbox_inches='tight') plt.close() # 使用示例 with rasterio.open('guangdong.tif') as src: data = src.read(1) plot_nightlight(data, '广东省2020年夜间灯光分布', 'gd_2020.png')

关键参数说明：

• LogNorm：对数归一化处理，使弱光区域更明显
• 自定义colormap：模拟专业遥感软件的夜间灯光配色
• dpi=300：保证出版级图像质量

4. 进阶分析与应用

4.1 灯光指数计算

通过简单的NumPy运算，我们可以从原始辐射值派生出有意义的指标：

def calculate_light_indices(data): # 替换无效值 data = np.where(data == -999, np.nan, data) # 计算各类指标 total_light = np.nansum(data) mean_light = np.nanmean(data) light_area = np.sum(data > 5) # 假设5为有效灯光阈值 # 灯光聚集度指数 top10 = np.nanpercentile(data, 90) concentration = np.nansum(data[data > top10]) / total_light return { '总辐射量': total_light, '平均辐射': mean_light, '亮区面积(像素)': light_area, '聚集度指数': concentration }

4.2 时间序列分析

对比不同年份数据可以分析城市发展动态。以下是核心比对方法：

def compare_years(file_2020, file_2015, region_mask): with rasterio.open(file_2020) as src: data_2020 = src.read(1) with rasterio.open(file_2015) as src: data_2015 = src.read(1) # 应用统一掩膜 mask = ~np.isnan(region_mask) data_2020 = np.where(mask, data_2020, np.nan) data_2015 = np.where(mask, data_2015, np.nan) # 计算变化指标 growth = np.nansum(data_2020) / np.nansum(data_2015) - 1 new_hotspots = np.sum((data_2020 > 5) & (data_2015 <= 5)) # 生成变化矩阵 change_matrix = np.zeros_like(data_2020) change_matrix[(data_2020 > 5) & (data_2015 <= 5)] = 1 # 新增亮区 change_matrix[(data_2020 <= 5) & (data_2015 > 5)] = -1 # 消失亮区 return { '增长率': growth, '新增亮区': new_hotspots, '变化矩阵': change_matrix }

4.3 性能优化技巧

处理省级以上规模数据时，内存可能成为瓶颈。以下是两个关键优化策略：

分块处理大文件

# 创建内存映射避免全量加载 with rasterio.open('large_file.tif') as src: block_shapes = src.block_shapes for ji, window in src.block_windows(): data = src.read(window=window) # 在此处理每个分块...

使用Zstandard压缩

# 在保存结果时启用高效压缩 profile = { 'driver': 'GTiff', 'compress': 'zstd', 'zstd_level': 3, 'predictor': 2, 'tiled': True, 'blockxsize': 256, 'blockysize': 256 } with rasterio.open('compressed.tif', 'w', **profile) as dst: dst.write(data)

5. 常见问题解决方案

5.1 坐标系转换问题

当需要将结果与其他数据集叠加时，可能需要坐标系转换：

from rasterio.warp import calculate_default_transform, reproject def reproject_raster(input_path, output_path, target_crs='EPSG:3857'): with rasterio.open(input_path) as src: transform, width, height = calculate_default_transform( src.crs, target_crs, src.width, src.height, *src.bounds ) metadata = src.meta.copy() metadata.update({ 'crs': target_crs, 'transform': transform, 'width': width, 'height': height }) with rasterio.open(output_path, 'w', **metadata) as dst: for i in range(1, src.count + 1): reproject( source=rasterio.band(src, i), destination=rasterio.band(dst, i), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=target_crs, resampling=rasterio.enums.Resampling.bilinear )

5.2 异常值处理

VIIRS数据中常见的异常值包括：

• 负值：通常表示填充值（如-999）
• 极高值：可能是闪电等临时光源
• NaN值：无效数据区域

标准化处理流程：

def clean_data(data, min_val=0, max_val=1000): data = np.where(data < min_val, np.nan, data) data = np.where(data > max_val, np.nan, data) return data

5.3 批量处理脚本

对于需要处理多年份/多区域的情况，建议构建自动化流水线：

import os from pathlib import Path def batch_process(input_dir, output_dir, shp_file): Path(output_dir).mkdir(exist_ok=True) for file in Path(input_dir).glob('*.tif'): year = file.stem.split('_')[3] # 从文件名提取年份 out_path = Path(output_dir) / f"processed_{year}.tif" try: # 执行裁剪和重投影 clip_and_reproject( str(file), str(out_path), shp_file ) print(f"成功处理: {file.name}") except Exception as e: print(f"处理失败 {file.name}: {str(e)}")

实际项目中，我们会将上述代码封装成命令行工具，配合日志系统和断点续处理功能，构建健壮的生产级数据处理流程。