哨兵（Sentinel-2）数据转.tif文件

简介

哨兵卫星是欧洲空间局（ESA）推出的一系列地球观测卫星，为全球提供了高分辨率、多光谱的遥感数据。但是原始数据格式不方便我进行后续分析，所以想要获取.tif格式的数据。

并且.tif格式的数据有以下几个优势：

• 广泛支持： .tif（Tagged Image File Format）是一种通用的图像文件格式，得到广泛支持。几乎所有的地理信息系统（GIS）软件和遥感分析工具都能够处理.tif文件。这使得.tif成为共享、存储和传输遥感数据的常用格式。这是重点，因为ArcMap和ENVI都可打开这个格式的数据。
• 保留地理信息： .tif文件支持嵌入地理信息，如投影坐标、仿射变换参数等。这样的信息对于空间分析和地图制图非常重要。在您的代码中，通过设置输出.tif文件的投影和仿射变换参数，确保了输出文件保留了原始数据的地理参考。
• 多波段数据合并： 遥感数据通常包含多个波段，每个波段对应于不同的光谱范围。将多个波段合并到一个.tif文件中，方便了数据的管理和分析。这在地学、生态学和农业等领域的研究中尤为重要。
• 易于处理： .tif格式的图像文件易于处理，可以使用各种图像处理工具进行操作。它支持压缩和无损压缩，使得文件大小合理，同时保持图像质量。
• 与其他数据集的兼容性： .tif格式在许多科学和工程应用中都是一种标准格式。使用.tif文件格式使得遥感数据更容易与其他地理、气象、和环境数据集集成，进而进行跨学科的分析。

在下载好的Sentinel-2数据中有个名为MTD_MSIL1C.xml的文件，这是元数据，包含了有关卫星数据集的详细信息。这个XML文件通常包括以下内容：

• 产品元数据： 包括卫星数据产品的标识符、生成日期、卫星名称等。
• 影像元数据： 包括波段信息、分辨率、影像格式等。
• 几何元数据： 包括地理参考信息、投影信息、仿射变换参数等。
• 辐射度校准元数据： 包括辐射度校准参数等。
• 质量控制信息： 包括云覆盖率、数据质量标志等。

转换代码也是基于这个元数据进行。

代码

代码融合了分辨率为10m的波段，可以根据需要修改。

代码1

参考：https://blog.csdn.net/KilllerQueen/article/details/114637970

from osgeo import gdal
import numpy as np

def convert_to_tiff(filename, output_path):
    # 打开栅格数据集
    root_ds = gdal.Open(filename)
    ds_list = root_ds.GetSubDatasets()  # 获取子数据集
    # print(df_list)  # 查看各个子集对应的波段详解
    visual_ds = gdal.Open(ds_list[0][0])  # 打开第一个数据子集的路径，ds_list[i][0]表示不同的子集，也就对应不同分辨率的波段。
    visual_arr = visual_ds.ReadAsArray()  # 将数据集中的数据读取为ndarray

    # 创建.tif文件
    band_count = visual_ds.RasterCount  # 波段数
    xsize = visual_ds.RasterXSize  # 列数
    ysize = visual_ds.RasterYSize  # 行数
    out_tif_name = output_path + "/" + filename.split(".SAFE")[0].split("/")[-2] + ".tif"  # 输出文件名
    driver = gdal.GetDriverByName("GTiff")  # 创建输出文件
    out_tif = driver.Create(out_tif_name, xsize, ysize, band_count, gdal.GDT_Float32)  # 创建输出文件
    out_tif.SetProjection(visual_ds.GetProjection())  # 设置投影坐标
    out_tif.SetGeoTransform(visual_ds.GetGeoTransform())  # 设置仿射变换参数

    for index, band in enumerate(visual_arr):  # 遍历每个波段
        band = np.array([band])  # 将每个波段的数据转换为ndarray
        for i in range(len(band[:])):
            # 数据写出
            out_tif.GetRasterBand(index + 1).WriteArray(band[i])

    out_tif.FlushCache()  # 最终将数据写入硬盘
    out_tif = None  # 注意必须关闭tif文件

# 调用示例
path = '遥感影像/S2B_MSIL1C_20191230T024119_N0208_R089_T51RUQ_20191230T042547.SAFE/'
output_path = '遥感影像/data processing'
filename = path + 'MTD_MSIL1C.xml'
convert_to_tiff(filename, output_path)  # 调用转换函数
print("转换完成")

代码2

from osgeo import gdal

def convert_to_tiff(input_file, output_file):
    ifile_name = 'MTD_MSIL1C.xml'  # 构造 MTD 文件名
    #  ofile_name = input_file.split(".SAFE")[0].split("/")[-1] + "_1.tif"  # 输出文件名
    ofile_name = input_file.split(".SAFE")[0].split("/")[-1] + ".tif"  # 输出文件名
    
    root_ds = gdal.Open(input_file+ifile_name)  # 打开 Sentinel-2 数据集
    ds_list = root_ds.GetSubDatasets()  # 获取子数据集列表
    vrt_ds = gdal.BuildVRT('merged.vrt', ds_list[0][0])  # 创建VRT文件
    gdal.Translate(output_file+ofile_name, vrt_ds, options=['-co', 'COMPRESS=NONE'])  # 将VRT文件转换为GeoTIFF文件，不压缩数据
    # 清理 - 关闭数据集
    vrt_ds = None
    root_ds = None

# 示例调用
input_file_path = '../Z_GIS_Data/遥感影像/Sentinel-2/S2B_MSIL1C_20191230T024119_N0208_R089_T51RUQ_20191230T042547.SAFE/'
output_file_path = '../Z_GIS_Data/遥感影像/Sentinel-2/'
convert_to_tiff(input_file_path, output_file_path)

代码分析

两个代码输出的文件大小差别过大，接近一倍：

1. 数据类型
   • 第一个代码片段使用了gdal.Create函数手动创建 GeoTIFF 文件，如果没有明确设置压缩选项，可能会生成较大的未压缩文件。
   • 第二个代码片段使用了gdal.Translate函数，该函数可能会在转换时默认应用某种压缩（我设置不压缩），并且它会尽量选择更高效的存储方式。
2. 波段处理
   • 第一个代码片段在迭代波段时，可能会有一些额外的数据处理或写入步骤，这可能导致生成的文件中包含一些额外的信息，影响文件大小。

我建议使用第二个代码，毕竟在没有影响数据质量的情况下，可以减少我们的硬盘开支。

简单的统计分析

转换完成我们可进行简单的统计分析：统计每个矢量面对应的影像像素均值和标准差。

# 像素统计
import geopandas as gpd
import rasterstats
import pandas as pd

# 设置文件路径
tif_file = '../Z_GIS_Data/遥感影像/Sentinel-2/S2B_MSIL1C_20191230T024119_N0208_R089_T51RUQ_20191230T042547.tif'
tif_file1 = '../Z_GIS_Data/遥感影像/Sentinel-2/S2B_MSIL1C_20191230T024119_N0208_R089_T51RUQ_20191230T042547_1.tif'
shp_file = '../Z_GIS_Data/遥感影像/Sentinel-2/徐汇AOI_fc.shp'

# 读取 Shapefile
shapefile = gpd.read_file(shp_file)

# 定义统计量函数
def calculate_stats(image, geometry):
    stats = rasterstats.zonal_stats(
        geometry,
        image,
        stats=["mean", "std"],
        nodata=0,  # 指定 NoData 值
        #affine=image.transform,  # 使用图像的仿射变换信息
        all_touched=True
    )
    return stats

# 循环处理每个要素
for idx, feature in shapefile.iterrows():
    geometry = feature.geometry
    stats = calculate_stats(tif_file, geometry)
    for band_stats in stats:
        mean = band_stats['mean']
        std = band_stats['std']
        shapefile.loc[idx, 'Mean'] = mean
        shapefile.loc[idx, 'Std'] = std

# 循环处理每个要素
for idx, feature in shapefile.iterrows():
    geometry = feature.geometry
    stats = calculate_stats(tif_file1, geometry)
    for band_stats in stats:
        mean = band_stats['mean']
        std = band_stats['std']
        shapefile.loc[idx, 'Mean1'] = mean
        shapefile.loc[idx, 'Std1'] = std

shapefile.to_file('../Z_GIS_Data/遥感影像/Sentinel-2/徐汇AOI_fc1.shp')
gdf = gpd.read_file('../Z_GIS_Data/遥感影像/Sentinel-2/徐汇AOI_fc1.shp')
gdf