# 夜间灯光数据实战：用Python快速处理VIIRS月度压缩包（含代码示例） 深夜的城市，从卫星的视角望去，是一片璀璨的光之海洋。这些由传感器捕捉到的点点光芒，早已超越了简单的夜景范畴，成为了经济学家、城市规划师、社会学家乃至环境研究者手中极具价值的量化工具。这就是VIIRS夜间灯光数据，它像一双永不疲倦的眼睛，持续记录着人类活动的强度与变迁。然而，对于许多数据分析师和GIS开发者而言，获取并处理这些原始数据却像是一场与时间和存储空间的赛跑——动辄数GB的月度压缩包、复杂的文件结构、专业的辐射值解析，每一步都可能成为阻碍深入分析的绊脚石。 如果你也曾面对下载下来的`SVDNB_npp_20231101-20231130_75N060E_vcmcfg_v10_20231201.tgz`这类文件感到无从下手，或者厌倦了手动解压、查看、转换的繁琐流程，那么这篇文章正是为你准备的。我们将彻底抛开对数据背景的冗长介绍，直接切入核心：**如何用Python构建一套自动化、高效率的VIIRS月度灯光数据处理流水线**。从压缩包的解压、GeoTIFF文件的读取与信息提取，到辐射值的批量计算与可视化初探，我将分享一套经过实战检验的代码方案，并穿插处理过程中那些容易被忽略的“坑”与技巧。我们的目标很明确：让你拿到数据后，能在最短时间内将其转化为可分析、可应用的干净数据矩阵。 ## 1. 环境准备与数据获取 在开始编写任何代码之前，一个稳定且功能齐全的Python环境是基石。我强烈建议使用`conda`来管理环境，因为它能很好地处理地理空间数据分析中常见的、依赖关系复杂的C库。我们将主要依赖`rasterio`、`gdal`、`numpy`和`pandas`这几个核心库。 > 注意：`rasterio`和`gdal`的安装有时会因为系统依赖而失败。在Ubuntu/Debian上，你可以先运行`sudo apt-get install gdal-bin libgdal-dev`；在macOS上，使用`brew install gdal`。然后再通过`pip`安装Python绑定，这能大幅提高成功率。 创建一个专属的conda环境并安装依赖： ```bash conda create -n viirs_processing python=3.9 conda activate viirs_processing conda install -c conda-forge rasterio geopandas numpy pandas matplotlib pip install py7zr # 用于处理.tgz或.7z压缩包 ``` 接下来是数据获取。虽然本文聚焦于处理，但一个完整的自动化流程理应包含下载环节。NOAA的官方FTP服务器是数据源，但其网络稳定性对国内用户并不友好。这里提供一个**备选思路**：利用云服务商（如AWS、GCP）的公开数据集镜像，或者通过学术机构的缓存服务器进行下载，速度会快很多。以下是一个使用`requests`库进行断点续传下载的稳健函数示例，它包含了重试机制和进度显示： ```python import requests import os from tqdm import tqdm def download_viirs_monthly(url, save_path, max_retries=5): """ 下载VIIRS月度数据压缩包，支持断点续传。 :param url: 数据文件直链 :param save_path: 本地保存路径 :param max_retries: 最大重试次数 """ headers = {} if os.path.exists(save_path): # 获取已下载文件大小，用于断点续传 downloaded_size = os.path.getsize(save_path) headers['Range'] = f'bytes={downloaded_size}-' else: downloaded_size = 0 for attempt in range(max_retries): try: with requests.get(url, headers=headers, stream=True, timeout=30) as r: r.raise_for_status() total_size = int(r.headers.get('content-length', 0)) + downloaded_size mode = 'ab' if downloaded_size else 'wb' # 续传用追加模式 with open(save_path, mode) as f, tqdm( desc=os.path.basename(save_path), total=total_size, unit='iB', unit_scale=True, initial=downloaded_size ) as pbar: for chunk in r.iter_content(chunk_size=8192): size = f.write(chunk) pbar.update(size) print(f"下载完成: {save_path}") break except requests.exceptions.RequestException as e: print(f"下载尝试 {attempt + 1} 失败: {e}") if attempt == max_retries - 1: raise ``` 这个函数不仅实用，还体现了生产级代码应有的健壮性。将下载链接列表化，配合此函数，即可轻松实现月度数据的批量自动获取。 ## 2. 自动化解压与文件结构解析 下载得到的通常是一个`.tgz`（tar.gz）压缩包。解压它看似简单，但里面藏着VIIRS数据组织的逻辑。一个压缩包解压后，并非直接得到GeoTIFF文件，而是可能先得到一个`.tar`文件，需要二次解压。更关键的是，我们需要准确识别出解压出的三个核心GeoTIFF文件： * `*avg_rade9h.tif`: 包含平均辐射值，单位是纳瓦/平方厘米/球面度（nW/cm²/sr），这是我们分析的核心。 * `*cf_cvg.tif`: 无云覆盖观测次数，用于评估平均辐射值的数据质量。 * `*cvg.tif`: 总观测次数（无论天气）。 手动操作既容易出错，也无法规模化。下面这个`ViirsMonthlyArchive`类，封装了整个解压与文件识别的过程： ```python import tarfile import os import re from pathlib import Path class ViirsMonthlyArchive: def __init__(self, archive_path, extract_to='./extracted'): self.archive_path = Path(archive_path) self.extract_to = Path(extract_to) self.extract_to.mkdir(parents=True, exist_ok=True) self.tif_files = {} # 存储识别出的文件路径 def extract(self): """解压.tgz文件，并处理内部可能的.tar文件""" print(f"正在解压: {self.archive_path.name}") with tarfile.open(self.archive_path, 'r:gz') as tar: # 安全解压，避免路径遍历漏洞 for member in tar.getmembers(): member_path = self.extract_to / Path(member.name).name if member.isfile(): tar.extract(member, path=self.extract_to) # 如果解压出的是.tar文件，进行二次解压 if member_path.suffix == '.tar': self._extract_inner_tar(member_path) self._identify_tif_files() return self def _extract_inner_tar(self, inner_tar_path): with tarfile.open(inner_tar_path, 'r') as inner_tar: inner_tar.extractall(path=self.extract_to) inner_tar_path.unlink() # 删除内部的.tar文件 def _identify_tif_files(self): """识别并分类解压后的GeoTIFF文件""" pattern = re.compile(r'.*_(avg_rade9h|cf_cvg|cvg)\.tif$') for tif_file in self.extract_to.rglob('*.tif'): match = pattern.match(tif_file.name) if match: file_type = match.group(1) # avg_rade9h, cf_cvg, cvg self.tif_files[file_type] = tif_file print(f"识别到 {file_type} 文件: {tif_file.name}") if len(self.tif_files) != 3: print("警告：未识别出全部三个核心GeoTIFF文件。") else: print("文件识别完成。") def get_file_path(self, file_type): """获取指定类型文件的路径""" return self.tif_files.get(file_type) # 使用示例 archive = ViirsMonthlyArchive('SVDNB_npp_20231101-20231130_75N060E_vcmcfg_v10_20231201.tgz') archive.extract() radiance_path = archive.get_file_path('avg_rade9h') ``` 通过这个类，我们不仅自动化了解压步骤，还精准地定位了后续处理所需的关键文件，为流水线作业打下了坚实基础。 ## 3. 核心辐射值数据的读取与处理 拿到`avg_rade9h.tif`文件后，真正的数据分析才开始。使用`rasterio`读取GeoTIFF既高效又能保留丰富的地理元数据。但直接读取整个图像到内存，对于覆盖大区域（如75N060E）的高分辨率数据可能是危险的。我们需要更精细的控制。 首先，让我们看看如何安全地读取数据并理解其基本属性： ```python import rasterio import numpy as np def inspect_viirs_tif(tif_path): """ 安全打开并查看VIIRS GeoTIFF文件的基本信息。 避免一次性加载大型数组导致内存溢出。 """ with rasterio.open(tif_path) as src: print("=== 文件元数据 ===") print(f"驱动: {src.driver}") print(f"图像尺寸 (宽 x 高): {src.width} x {src.height}") print(f"波段数量: {src.count}") print(f"数据类型: {src.dtypes[0]}") print(f"NoData值: {src.nodata}") print(f"仿射变换矩阵:\n{src.transform}") print(f"坐标系 (CRS): {src.crs}") # 读取一个小的窗口样本，了解数据范围 sample_window = rasterio.windows.Window(0, 0, 100, 100) sample_data = src.read(1, window=sample_window) print(f"\n=== 数据样本 (前100x100像素) ===") print(f"形状: {sample_data.shape}") print(f"有效像素数: {np.count_nonzero(~np.isnan(sample_data))}") print(f"最小值: {np.nanmin(sample_data):.4f}") print(f"最大值: {np.nanmax(sample_data):.4f}") print(f"平均值: {np.nanmean(sample_data):.4f}") return src # 使用函数 src = inspect_viirs_tif(radiance_path) ``` 对于全幅数据的处理，我推荐使用**分块读取（Block Processing）**。`rasterio`本身支持基于块的迭代，这对于后续的统计计算或逐块应用算法（如去噪、插值）至关重要。 ```python def process_by_blocks(tif_path, chunk_size=1024): """ 按块处理大型GeoTIFF，适用于统计或逐像素运算。 :param chunk_size: 处理块的大小（像素） """ stats = {'total_pixels': 0, 'sum_radiance': 0.0, 'valid_pixels': 0} with rasterio.open(tif_path) as src: # 获取图像的分块结构 block_shapes = src.block_shapes print(f"原始数据块形状: {block_shapes}") # 自定义更大的块进行迭代，提高I/O效率 for ji, window in src.block_windows(1): data = src.read(1, window=window) valid_mask = ~np.isnan(data) & (data != src.nodata) valid_data = data[valid_mask] stats['total_pixels'] += data.size stats['valid_pixels'] += np.sum(valid_mask) stats['sum_radiance'] += np.sum(valid_data) # 可以在此处添加更复杂的处理逻辑，例如： # - 阈值过滤（剔除极低值，可能是噪声） # - 简单分类（城市、乡村、无光区） # - 实时计算并写入新的输出文件 if stats['valid_pixels'] > 0: stats['mean_radiance'] = stats['sum_radiance'] / stats['valid_pixels'] else: stats['mean_radiance'] = np.nan print(f"处理完成。有效像素比例: {stats['valid_pixels']/stats['total_pixels']*100:.2f}%") print(f"全局平均辐射值: {stats['mean_radiance']:.4f} nW/cm²/sr") return stats ``` 此外，我们经常需要将辐射值数据与其他矢量数据（如行政区划）进行关联分析。这就需要用到**栅格统计（Zonal Statistics）**。虽然`rasterstats`库很流行，但理解其原理并手动实现关键部分能给你更多控制权。下面的示例展示了如何为一个GeoDataFrame中的每个多边形计算平均灯光强度： ```python import geopandas as gpd from shapely.geometry import box import rasterio.mask def zonal_stats_viirs(tif_path, gdf, stats_field='mean_rad'): """ 为GeoDataFrame中的每个几何图形计算VIIRS辐射值的区域统计。 :param tif_path: VIIRS辐射值GeoTIFF路径 :param gdf: 包含多边形几何体的GeoDataFrame :param stats_field: 存储结果的字段名 """ results = [] with rasterio.open(tif_path) as src: for idx, row in gdf.iterrows(): geom = row.geometry try: # 1. 裁剪出多边形区域的栅格数据 out_image, out_transform = rasterio.mask.mask(src, [geom], crop=True, nodata=np.nan) # out_image形状为 (1, H, W) region_data = out_image[0] # 2. 计算统计值（忽略NaN） valid_data = region_data[~np.isnan(region_data)] if valid_data.size > 0: mean_val = np.mean(valid_data) # 还可以计算总和、标准差、最大值等 else: mean_val = np.nan results.append(mean_val) except Exception as e: print(f"处理要素 {idx} 时出错: {e}") results.append(np.nan) gdf[stats_field] = results return gdf # 假设有一个包含中国各省边界的GeoDataFrame `china_provinces_gdf` # china_provinces_with_light = zonal_stats_viirs(radiance_path, china_provinces_gdf.copy()) ``` ## 4. 质量评估与数据清洗策略 直接使用`avg_rade9h`的辐射值进行分析可能会引入误差，因为原始数据中存在**无效值、背景噪声和异常高值**。`cf_cvg`（无云覆盖次数）文件是评估数据质量的关键。通常，观测次数越多，平均辐射值的可靠性越高。 一个实用的策略是，基于`cf_cvg`创建一个**质量掩膜（Quality Mask）**。例如，我们可以设定一个阈值，只保留那些无云观测次数足够多的像素。 ```python def create_quality_mask(radiance_path, cf_cvg_path, min_clear_obs=5): """ 根据无云观测次数创建质量掩膜。 :param min_clear_obs: 认为数据可靠所需的最小无云观测次数 :return: 布尔数组掩膜，True表示高质量像素 """ with rasterio.open(radiance_path) as rad_src, rasterio.open(cf_cvg_path) as cvg_src: # 确保两个文件地理对齐（通常是一致的） radiance_data = rad_src.read(1) cf_cvg_data = cvg_src.read(1) # 创建掩膜：辐射值有效 且 无云观测次数达标 valid_radiance = ~np.isnan(radiance_data) & (radiance_data != rad_src.nodata) sufficient_coverage = cf_cvg_data >= min_clear_obs quality_mask = valid_radiance & sufficient_coverage print(f"高质量像素比例: {np.sum(quality_mask) / quality_mask.size * 100:.2f}%") return quality_mask # 获取cf_cvg文件路径 cf_cvg_path = archive.get_file_path('cf_cvg') quality_mask = create_quality_mask(radiance_path, cf_cvg_path, min_clear_obs=3) ``` 除了云覆盖，夜间灯光数据中常见的“**过饱和**”现象（城市中心像素值达到传感器上限）和“**背景噪声**”（偏远地区极低的非零值）也需要处理。这里没有放之四海而皆准的方法，需要根据研究区域和目的进行调整。一个常见的去噪方法是使用**百分比阈值法**。 ```python def denoise_radiance(radiance_data, lower_percentile=0.5, upper_percentile=99.95): """ 使用百分比阈值法进行简单的数据清洗。 :param radiance_data: 辐射值数组 :param lower_percentile: 低于此百分位的值被视为噪声，设为NaN :param upper_percentile: 高于此百分位的值可能过饱和，可裁剪或单独处理 :return: 清洗后的数组 """ valid_data = radiance_data[~np.isnan(radiance_data)] if len(valid_data) == 0: return radiance_data lower_thresh = np.percentile(valid_data, lower_percentile) upper_thresh = np.percentile(valid_data, upper_percentile) cleaned_data = radiance_data.copy() # 将极低值设为NaN cleaned_data[cleaned_data < lower_thresh] = np.nan # 对过高值进行裁剪（或记录日志供后续分析） cleaned_data[cleaned_data > upper_thresh] = upper_thresh print(f"噪声阈值: {lower_thresh:.4f}, 过饱和阈值: {upper_thresh:.4f}") print(f"被清理的噪声像素比例: {np.sum(radiance_data < lower_thresh) / radiance_data.size * 100:.4f}%") return cleaned_data ``` 将质量掩膜与去噪步骤结合，我们就能得到一份相对干净、可靠的辐射值数据集，为后续的经济指标估算、城市扩张监测等应用提供高质量输入。 ## 5. 从数据到洞察：实用案例与输出 处理干净的辐射值数据本身只是一个数字矩阵，它的价值在于与应用场景的结合。这里分享两个快速将数据转化为洞察的实用案例。 **案例一：快速生成区域灯光强度变化时间序列** 假设你已经按上述流程处理了连续多个月的数据，并计算出了每个区域（如城市）的平均辐射值。接下来，用`pandas`和`matplotlib`可以轻松实现可视化。 ```python import pandas as pd import matplotlib.pyplot as plt import seaborn as sns # 假设我们有一个DataFrame `city_light_series`，索引是时间，列是不同的城市 # 数据形如： # date | City_A | City_B | City_C # 2023-01-01 | 12.5 | 8.2 | 5.1 # 2023-02-01 | 13.1 | 8.0 | 5.3 # ... def plot_light_trends(light_df): """ 绘制多城市夜间灯光强度时间序列图。 """ plt.figure(figsize=(14, 7)) for city in light_df.columns: plt.plot(light_df.index, light_df[city], marker='o', linewidth=2, label=city) plt.title('VIIRS月度夜间灯光平均辐射强度趋势', fontsize=16, pad=20) plt.xlabel('日期', fontsize=12) plt.ylabel('平均辐射强度 (nW/cm²/sr)', fontsize=12) plt.grid(True, linestyle='--', alpha=0.7) plt.legend(title='城市', bbox_to_anchor=(1.05, 1), loc='upper left') plt.xticks(rotation=45) plt.tight_layout() plt.show() # 计算环比增长率 monthly_growth = light_df.pct_change() * 100 print("近三个月部分城市灯光环比增长率 (%):") print(monthly_growth.tail(3).round(2)) ``` **案例二：生成可供GIS软件使用的标准GeoTIFF** 经过清洗和计算的数据，我们通常需要导回为GeoTIFF，以便在QGIS或ArcGIS中进行空间分析与制图。使用`rasterio`写入数据时，务必保留原始的地理参照信息。 ```python def export_cleaned_tif(input_tif_path, cleaned_data_array, output_path, quality_mask=None): """ 将处理后的数据数组导出为新的GeoTIFF文件。 :param input_tif_path: 原始TIFF路径，用于复制元数据 :param cleaned_data_array: 处理后的numpy数组 :param output_path: 输出路径 :param quality_mask: 可选，质量掩膜。掩膜为False的位置将被设为指定的NoData值。 """ with rasterio.open(input_tif_path) as src: profile = src.profile.copy() # 更新元数据，例如数据类型、压缩方式等 profile.update( dtype=rasterio.float32, nodata=np.nan, compress='lzw' # 使用LZW压缩以减小文件体积 ) output_data = cleaned_data_array.astype(np.float32) if quality_mask is not None: # 将低质量像素设为NaN output_data[~quality_mask] = np.nan with rasterio.open(output_path, 'w', **profile) as dst: dst.write(output_data, 1) print(f"已导出清洗后的数据至: {output_path}") print(f"输出文件CRS: {profile['crs']}, 变换矩阵: {profile['transform']}") # 使用示例 # cleaned_radiance = denoise_radiance(original_radiance_data) # export_cleaned_tif(radiance_path, cleaned_radiance, './output/cleaned_radiance_202311.tif', quality_mask) ``` 在整个流程的最后，别忘了**数据与代码的整理**。为每个月的处理创建一个独立的输出目录，里面包含原始的辐射值、质量掩膜、清洗后的数据以及关键的统计日志。这不仅能让你在几个月后还能清晰地复现结果，也是团队协作的基础。我自己的项目目录通常是这样组织的： ``` viirs_processing_project/ ├── data/ │ ├── raw/ # 存放原始.tgz压缩包 │ └── extracted/ # 存放解压后的原始.tif文件 ├── scripts/ │ ├── 01_download.py │ ├── 02_extract.py │ ├── 03_process.py # 包含本文大部分函数 │ └── 04_analyze.py ├── output/ │ ├── 202311/ │ │ ├── cleaned_radiance.tif │ │ ├── quality_mask.tif │ │ └── stats_report.json │ └── 202312/ └── config.yaml # 存放路径、参数等配置 ``` 处理VIIRS夜间灯光数据，从令人望而生畏的压缩包到生成清晰的洞察图表，这个过程确实有诸多细节。但一旦用Python将流程固化下来，你会发现它变成了一个稳定、高效且可复用的数据管道。这套代码框架我已经在多个区域经济监测项目中实际应用过，它节省了大量的手动操作时间，让我能更专注于数据本身所揭示的故事——那些隐藏在光影变化背后的经济活力与城市脉搏。