利用Python自动化处理Sentinel2影像：从SAFE格式到GeoTIFF的高效转换

1. 为什么需要处理Sentinel2的SAFE格式数据

如果你接触过遥感影像处理，大概率听说过Sentinel2卫星数据。作为欧洲航天局哥白尼计划的重要成员，Sentinel2提供了全球覆盖的多光谱影像，广泛应用于农业监测、环境评估、灾害管理等领域。但当你兴冲冲下载了数据后，可能会遇到第一个难题：这些.SAFE格式的文件该怎么处理？

SAFE（Standard Archive Format for Europe）是欧空局专门为Sentinel系列卫星设计的数据格式。它本质上是一个文件夹，里面包含了XML元数据文件和多个JPEG2000格式的影像数据。这种设计虽然规范，但直接使用起来并不方便：

1. 专业软件限制：通常需要SNAP这类专业工具处理，但软件体积大、学习成本高
2. 格式兼容性问题：JPEG2000虽然压缩率高，但很多通用GIS软件支持有限
3. 批量处理困难：当需要处理数十甚至上百景影像时，手动操作效率极低

我在处理农业遥感项目时就深有体会：每次下载完数据都要花大量时间转换格式，直到开发了Python自动化脚本才彻底解决这个问题。下面我就分享如何用Python+GDAL实现高效转换。

2. 准备工作与环境配置

2.1 安装必要的Python库

这套方案的核心是GDAL（Geospatial Data Abstraction Library），它是一个开源的地理空间数据处理库。建议通过conda安装，能自动解决依赖问题：

conda create -n sentinel python=3.8 conda activate sentinel conda install -c conda-forge gdal numpy

如果不用conda，也可以用pip安装，但可能需要先安装GDAL的系统依赖：

# Ubuntu系统 sudo apt-get install libgdal-dev pip install GDAL==$(gdal-config --version) numpy # Windows系统 # 建议下载预编译的GDAL wheel文件

2.2 数据目录结构说明

下载的Sentinel2 L2A数据通常长这样：

S2A_MSIL2A_20230601T100031_N0509_R122_T33TUM_20230601T134559.SAFE/ ├── GRANULE/ │ └── L2A_T33TUM_A037018_20230601T100031/ │ ├── IMG_DATA/ │ │ ├── R10m/ │ │ │ ├── T33TUM_20230601T100031_B02_10m.jp2 │ │ │ ├── T33TUM_20230601T100031_B03_10m.jp2 │ │ │ └── ... │ │ └── R20m/ │ │ ├── T33TUM_20230601T100031_B05_20m.jp2 │ │ └── ... │ └── MTD_TL.xml └── MTD_MSIL2A.xml

关键文件是MTD_MSIL2A.xml，它相当于整个数据集的"说明书"，GDAL就是通过它来理解数据结构的。

3. 核心转换代码解析

3.1 读取SAFE格式数据

GDAL的神奇之处在于它能直接读取SAFE格式的XML文件，自动解析内部结构：

from osgeo import gdal def read_sentinel2(safe_path): # 构造XML文件路径 xml_path = f"{safe_path}/MTD_MSIL2A.xml" # 打开主元数据文件 root_ds = gdal.Open(xml_path) if not root_ds: raise ValueError(f"无法打开文件: {xml_path}") # 获取所有子数据集信息 sub_datasets = root_ds.GetSubDatasets() return sub_datasets

运行后会得到一个子数据集列表，每个元素是一个元组，格式如下：

[ ('SENTINEL2_L2A:/path/to/MTD_MSIL2A.xml:10m:EPSG_32632', 'Bands B2, B3, B4, B8 with 10m resolution'), ('SENTINEL2_L2A:/path/to/MTD_MSIL2A.xml:20m:EPSG_32632', 'Bands B5, B6, B7, B8A, B11, B12 with 20m resolution'), ... ]

3.2 选择需要的波段组合

Sentinel2数据按分辨率分组存储，常见的有：

• 10米分辨率：B2(蓝)、B3(绿)、B4(红)、B8(近红外)
• 20米分辨率：B5-B7、B8A、B11-B12
• 60米分辨率：B1、B9、B10

假设我们需要10米分辨率的自然色合成图像：

def get_visual_bands(sub_datasets): # 查找包含10m波段的子数据集 visual_dataset = None for desc in sub_datasets: if '10m' in desc[1] and 'B2' in desc[1]: visual_dataset = gdal.Open(desc[0]) break if not visual_dataset: raise ValueError("未找到10m分辨率波段数据") # 读取为numpy数组 band_data = visual_dataset.ReadAsArray() return band_data, visual_dataset

3.3 写入GeoTIFF文件

将读取的数据写入标准GeoTIFF格式：

def write_geotiff(output_path, data, template_ds): driver = gdal.GetDriverByName("GTiff") # 获取原始数据参数 bands, height, width = data.shape dtype = gdal.GDT_UInt16 # Sentinel2数据通常为16位无符号整型 # 创建新文件 out_ds = driver.Create( output_path, width, height, bands, dtype ) # 设置地理参考信息 out_ds.SetGeoTransform(template_ds.GetGeoTransform()) out_ds.SetProjection(template_ds.GetProjection()) # 写入每个波段 for i in range(bands): out_band = out_ds.GetRasterBand(i+1) out_band.WriteArray(data[i]) out_band.FlushCache() out_ds = None # 关闭文件

4. 批量处理与实用技巧

4.1 自动化批量转换

实际项目中我们往往需要处理大量数据，这个脚本可以实现全自动批量处理：

import glob import os def batch_convert(input_folder, output_folder): # 确保输出目录存在 os.makedirs(output_folder, exist_ok=True) # 查找所有SAFE文件夹 safe_files = glob.glob(os.path.join(input_folder, "*.SAFE")) for safe_path in safe_files: try: # 生成输出文件名 base_name = os.path.basename(safe_path).replace(".SAFE", ".tif") output_path = os.path.join(output_folder, base_name) # 执行转换 sub_datasets = read_sentinel2(safe_path) band_data, template_ds = get_visual_bands(sub_datasets) write_geotiff(output_path, band_data, template_ds) print(f"成功转换: {safe_path}") except Exception as e: print(f"处理{safe_path}时出错: {str(e)}")

4.2 分辨率处理技巧

如果需要不同分辨率的数据融合，比如将20m波段重采样到10m分辨率：

def resample_bands(sub_datasets, target_resolution="10m"): # 获取参考数据集（10m） ref_ds = None for desc in sub_datasets: if target_resolution in desc[1]: ref_ds = gdal.Open(desc[0]) break # 重采样其他波段 resampled_bands = [] for desc in sub_datasets: if target_resolution not in desc[1]: ds = gdal.Open(desc[0]) # 使用gdal.Warp进行重采样 resampled = gdal.Warp("", ds, format="MEM", xRes=10, yRes=10, resampleAlg="bilinear", dstSRS=ref_ds.GetProjection()) resampled_bands.append(resampled.ReadAsArray()) return np.concatenate(resampled_bands)

4.3 常见问题排查

1. GDAL版本问题：确保使用GDAL 2.4+版本，旧版本可能不支持SAFE格式
2. 内存不足：处理大区域数据时，可以使用分块读取：

   # 分块读取数据 chunk_size = 1024 data = visual_dataset.ReadAsArray( xoff=0, yoff=0, xsize=visual_dataset.RasterXSize, ysize=visual_dataset.RasterYSize, buf_xsize=chunk_size, buf_ysize=chunk_size )

3. 中文路径问题：如果路径包含中文，设置环境变量：

   os.environ['GDAL_FILENAME_IS_UTF8'] = 'YES'

5. 进阶应用示例

5.1 生成NDVI植被指数

转换为GeoTIFF后，可以方便地进行各种分析。比如计算NDVI（归一化植被指数）：

def calculate_ndvi(red_band_path, nir_band_path): # 读取红波段和近红外波段 red_ds = gdal.Open(red_band_path) nir_ds = gdal.Open(nir_band_path) red = red_ds.GetRasterBand(1).ReadAsArray().astype(float) nir = nir_ds.GetRasterBand(1).ReadAsArray().astype(float) # 计算NDVI ndvi = (nir - red) / (nir + red + 1e-10) # 避免除以零 # 保存结果 write_geotiff("ndvi.tif", ndvi[np.newaxis, ...], red_ds)

5.2 与Shapely结合进行区域统计

结合地理空间分析库，可以统计特定区域的数值：

from shapely.geometry import shape import rasterstats def zonal_stats(geojson_path, raster_path): # 读取GeoJSON中的几何图形 with open(geojson_path) as f: features = json.load(f)['features'] geometries = [shape(feat['geometry']) for feat in features] # 计算区域统计 stats = rasterstats.zonal_stats( geometries, raster_path, stats=['mean', 'min', 'max', 'median'] ) return stats

这套自动化处理方法已经在我参与的多个农业遥感项目中得到验证，平均处理效率比手动操作提升20倍以上。特别是在处理季节性连续监测数据时，优势更加明显——原本需要数天完成的工作，现在喝杯咖啡的时间就能搞定。