简介
哨兵卫星是欧洲空间局（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() # 获取子数据集列表
# ds_list[0]、ds_list[1]和ds_list[2]分别对应10m、20m和60m分辨率的波段
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)
代码分析
两个代码输出的文件大小差别过大，接近一倍：
数据类型
第一个代码片段使用了gdal.Create函数手动创建 GeoTIFF 文件，如果没有明确设置压缩选项，可能会生成较大的未压缩文件。第二个代码片段使用了gdal.Translate函数，该函数可能会在转换时默认应用某种压缩（我设置不压缩），并且它会尽量选择更高效的存储方式。 波段处理
第一个代码片段在迭代波段时，可能会有一些额外的数据处理或写入步骤，这可能导致生成的文件中包含一些额外的信息，影响文件大小。
我建议使用第二个代码，毕竟在没有影响数据质量的情况下，可以减少我们的硬盘开支。
简单的统计分析
转换完成我们可进行简单的统计分析：统计每个矢量面对应的影像像素均值和标准差。
# 像素统计
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