由于Sentinel影像的存储格式特殊,通常需要特殊的软件(如SNAP)来处理,有局限性。于是使用GDAL对其进行读取,并重新写为.tif格式,方便进一步处理。
影像数据用的是大气校正后的L2A级数据,当然L1C数据也可以。打开SAFE后缀的文件夹,可以看到里面的内容,jp2格式的图片数据是存储在GRANULE文件夹中,这里不需要知道它的存储规则,因为GDAL可以直接读取最下面的.xml文件(注意GDAL的版本,2.4以下的版本可能不支持)。
文件的打开方式与其他格式的遥感影像相同,用gdal.Open(filename)
即可读入。返回结果是一个list,list中的每个元素是一个tuple,每个tuple中包含了对数据集的路径,元数据等的描述信息。获取其中的子数据集,使用data.GetSubDatasets()
,这里面就包含了子数据集的路径,元数据等的描述信息。
ds_list = root_ds.GetSubDatasets() # 获取子数据集信息
然后用子数据集的路径打开,读取为ndarry格式。这里第1个子集存储的是10m分辨率的B2, B3, B4, B8这4个波段,如果需要其他波段,就改ds_list[0][0]换别的子集,如ds_list[1][0]是第2个子集的路径,ds_list[1][1]是第2个子集的描述信息。
visual_ds = gdal.Open(ds_list[0][0]) # 打开第1个数据子集的路径。
visual_arr = visual_ds.ReadAsArray() # 将数据集中的数据读取为ndarray
然后要写新的tif图像,包括图像数据和投影等,具体可以参考以下代码。
# -*- coding: utf-8 -*-
# @File : S2_read.py
# @Author: Freezinghot
# @Date : 2021/3/9
# @Desc : 读取Sentinel2并转为Geo_tiff
from osgeo import gdal
import os
import numpy as np
from osgeo import gdal, osr, ogr
import glob
# os.environ['CPL_ZIP_ENCODING'] = 'UTF-8'
def S2tif(filename):
# 打开栅格数据集
print(filename)
root_ds = gdal.Open(filename)
print(type(root_ds))
# 返回结果是一个list,list中的每个元素是一个tuple,每个tuple中包含了对数据集的路径,元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
ds_list = root_ds.GetSubDatasets() # 获取子数据集。该数据以数据集形式存储且以子数据集形式组织
visual_ds = gdal.Open(ds_list[0][0]) # 打开第1个数据子集的路径。ds_list有4个子集,内部前段是路径,后段是数据信息
visual_arr = visual_ds.ReadAsArray() # 将数据集中的数据读取为ndarray
# 创建.tif文件
band_count = visual_ds.RasterCount # 波段数
xsize = visual_ds.RasterXSize
ysize = visual_ds.RasterYSize
out_tif_name = filename.split(".SAFE")[0] + ".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])
for i in range(len(band[:])):
# 数据写出
out_tif.GetRasterBand(index + 1).WriteArray(band[i]) # 将每个波段的数据写入内存,此时没有写入硬盘
out_tif.FlushCache() # 最终将数据写入硬盘
out_tif = None # 注意必须关闭tif文件
if __name__ == "__main__":
from osgeo import gdal
SAFE_Path = (r'E:\RSDATA\Sentinel2\L2A')
data_list = glob.glob(SAFE_Path + "\\*.SAFE")
#filename = ('E:\\RSDATA\\Sentinel2\\L2A\\S2A_MSIL2A_20210220T024731_N9999_R132_T51STA_20210306T024402.SAFE\\MTD_MSIL2A.xml')
for i in range(len(data_list)):
data_path = data_list[i]
filename = data_path + "\\MTD_MSIL2A.xml"
S2tif(filename)
print(data_path + "-----转tif成功")
print("----转换结束----")
转换结束,在目标路径下生成同名的tif图像。如果有多幅图像,也可以实现批量处理。
可以直接拖到ENVI查看波段组合的效果,是保存了4个10m分辨率的波段。
10m分辨率真彩色:
—————————————————————————————————————————————
增加内容:
如果需要其他波段的数据又不知道数据结构的话,可以debug模式看一下文件的结构。下面是我找到的子数据集的结构,位于程序中ds_list中。如图所示,假如要获得20m分辨率的波段集合,只需要修改这一句:
visual_ds = gdal.Open(ds_list[1][0])
20m分辨率(B12,B11,B8A)组合:
欢迎在评论区交流,如果对您有帮助,请点个赞哦~