前言

AI Earth提供了L2A级产品。此外，GEE PIE平台也可以进行检索下载。但是！我花了两天时间，这两个平台都不好用，很卡很卡，挂了梯子也卡，只能放弃。由于哨兵二号一些数据下架了，需要激活，而两次激活之间需要间隔一定时间，md，导致我批量下载的时候巨慢！

阿里云平台的AI Earth还是好用，下载速度快，下架产品激活快，缺点是不能批量勾选，很难想象这个平台没有批量勾选操作功能！算了算了，下载快就行了。

花了一天时间，写了Sentinel-2批处理的相关程序。哨兵二号下载的文件夹嵌套令人泪目啊。

1、批量解压，按日期存放

2、批量转换为tif格式，按日期、分辨率存放

3、批量融合波段，按日期、分辨率存放

4、批量拼接/裁剪，按日期、分辨率存放

写好这些代码后，后面有需求简单改改就行了。

批量解压，按日期存放

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#!/usr/bin/env python 
# -*- coding: utf-8 -*-
# @Time    : 2023/9/11 15:31
# @Author  : huangxiaohan
# @Site    : 
# @File    : batch_unzip_sentinel.py
# @Software: PyCharm

# https://zhuanlan.zhihu.com/p/531506915?utm_id=0


import os
import zipfile

"""
该模块功能：
只需要提供压缩文件（zip格式）所在的路径，代码会自动遍历该路径下
所有的zip文件，并将其解压缩，保存的指定的保存路径中
解压后按照日期存放（可能需要同一多幅影像进行拼接）
"""


def un_zip(filename, output_dir):
    """
    解压zip文件
    :param filename: 解压文件路径
    :param output_dir: 输出路径
    :return:
    """
    if not os.path.exists(filename):
        print("您输入的压缩包文件不存在，请检查路径！")
        return
    else:
        zip_file = zipfile.ZipFile(filename)
        members = zip_file.namelist()
        decompress_dir = output_dir
        for member in members:
            zip_file.extract(member, decompress_dir)
        zip_file.close()


if __name__ == '__main__':
    output_dir = 'F:/PoyangLake/poyanghu_data/test/unzip/'
    zipfile_path = 'F:/PoyangLake/poyanghu_data/test/origin/'

    # 建立所有日期的文件夹
    for file in os.listdir(zipfile_path):
        s2_date = file.split('_')[2][:8]
        out_dir = os.path.join(output_dir,s2_date)
        if not os.path.exists(out_dir):
            os.mkdir(out_dir)

    # 解压，按日期存放
    num = len(os.listdir(zipfile_path))
    n = 0
    for file in os.listdir(zipfile_path):
        n = n+1
        if os.path.splitext(os.path.basename(file))[1] == '.zip':
            print(file,'解压中...','进度:',n,'---',num)
            filepath = os.path.join(zipfile_path, file)
            # 截取日期  20220804
            s2_date = file.split('_')[2][:8]
            # 放到对应文件夹
            output_path = os.path.join(output_dir,s2_date)
            un_zip(filepath, output_path)

批量转换为tif格式，按日期、分辨率存放

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#!/usr/bin/env python 
# -*- coding: utf-8 -*-
# @Time    : 2023/9/11 15:46
# @Author  : huangxiaohan
# @Site    : 
# @File    : batch_jp2tif_sentinel.py
# @Software: PyCharm


'''
按日期文件夹进行读取
批量转换为tif格式后
在日期文件夹下建立/R10m /R20m /R60米子文件夹进行存放
'''
import os
from osgeo import gdal
import numpy as np

def jp2Totif(jp2_path,save_path,file_name):
    file_name = os.path.splitext(file_name)[0]
    save_file = os.path.join(save_path,file_name)
    save_file = save_file+'.tif'
    dataset = gdal.Open(jp2_path)
    rows = dataset.RasterYSize
    cols = dataset.RasterXSize
    projection = dataset.GetProjection()
    trans = dataset.GetGeoTransform()
    data = dataset.ReadAsArray()
    if data.dtype == 'uint16':
        driver = gdal.GetDriverByName('GTiff')
        out_dataset = driver.Create(save_file, cols, rows, 1, gdal.GDT_UInt16)
        out_dataset.SetProjection(projection)
        out_dataset.SetGeoTransform(trans)
        out_dataset.GetRasterBand(1).WriteArray(data)
        out_dataset.GetRasterBand(1).SetNoDataValue(0)
        out_dataset.FlushCache()
        del dataset, out_dataset
    elif data.dtype == 'uint8':
        driver = gdal.GetDriverByName('GTiff')
        out_dataset = driver.Create(save_file, cols, rows, 1, gdal.GDT_Byte)
        out_dataset.SetProjection(projection)
        out_dataset.SetGeoTransform(trans)
        out_dataset.GetRasterBand(1).WriteArray(data)
        out_dataset.GetRasterBand(1).SetNoDataValue(0)
        out_dataset.FlushCache()
        del dataset, out_dataset




if __name__ == "__main__":

    input_dir = 'F:/PoyangLake/poyanghu_data/test/unzip/'
    output_dir = 'F:/PoyangLake/poyanghu_data/test/tif/'
    # 建立日期文件夹
    for date in os.listdir(input_dir):
        date_dir = os.path.join(output_dir,date) # ././tif/20220804
        os.mkdir(date_dir)
        os.mkdir(date_dir + '/' + 'R10m')
        os.mkdir(date_dir + '/' + 'R20m')
        os.mkdir(date_dir + '/' + 'R60m')
    # 日期遍历
    for date in os.listdir(input_dir):
        date_dir = os.path.join(input_dir,date) # ./././20220804
        for subdir in os.listdir(date_dir):
            if os.path.isdir(os.path.join(date_dir, subdir)):
                root_subdir = os.path.join(date_dir, subdir)  # ./././20220804/S2文件夹
                subdir2 = root_subdir + '/'+'GRANULE' # ./././20220804/S2文件夹/GRANULE
                subdir3 = subdir2 + '/' + os.listdir(subdir2)[0]  # ./././20220804/S2文件夹/GRANULE/L2A_T50RLT_A037166_20220804T025911/
                subdir4 = subdir3 + '/' + 'IMG_DATA' # ./././20220804/S2文件夹/GRANULE/L2A_T50RLT_A037166_20220804T025911/IMG_DATA
                for k in os.listdir(subdir4):
                    # 此时的k为各个分辨率文件夹名R10m R20m R60m
                    subdir5 = subdir4 + '/' + k  # ./././R10m
                    for l in os.listdir(subdir5):
                        # l为各个jp2文件
                        band = l.split('_')[2]
                        # 需要转换的波段
                        if band == 'B04' or band == 'B03' or band == 'B02':
                            jp2_path = subdir5+'/'+ l
                            print('jp2_path:', jp2_path)
                            # 设置输出路径 tif/20220804/R10m
                            out_path = os.path.join(output_dir,date,k)
                            jp2Totif(jp2_path, out_path, l)

批量融合波段，按日期、分辨率存放

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#!/usr/bin/env python 
# -*- coding: utf-8 -*-
# @Time    : 2023/9/11 20:07
# @Author  : huangxiaohan
# @Site    : 
# @File    : batch_bandcom_sentinel.py
# @Software: PyCharm
'''
将相同分辨率的单波段影像进行融合

目录结构
/
---20220824/
        ---R10m/T50RLT_20220906T025541_B02_10m.tif、T50RMS_20220906T025541_B02_10m.tif
        ---R20m/
        ---R60m/

同一文件夹内，文件的前6个字符相同，进行波段组合
'''


import os
from osgeo import gdal
import numpy as np


def band_classify(band_names):
    '''
    按照区域进行分类
    :param band_name:
    :return:
    '''
    classified_strings = {}
    for string in band_names:
        key = string[:6]  # 假设按照前6个字符分类
        if key in classified_strings:
            classified_strings[key].append(string)
        else:
            classified_strings[key] = [string]
    return classified_strings

def bandcombination(input_path,out_path):
    for date in os.listdir(input_path):
        print('日期',date)
        r10m_dir = input_path + '/' + date + '/' + 'R10m'
        r20m_dir = input_path + '/' + date + '/' + 'R20m'
        r60m_dir = input_path + '/' + date + '/' + 'R60m'

        r10m_bands = os.listdir(r10m_dir)
        r20m_bands = os.listdir(r20m_dir)
        r60m_bands = os.listdir(r60m_dir)
        # 过滤下tif
        r10m_bands = [x for x in r10m_bands if os.path.splitext(x)[1] == '.tif']
        r20m_bands = [x for x in r20m_bands if os.path.splitext(x)[1] == '.tif']
        r60m_bands = [x for x in r60m_bands if os.path.splitext(x)[1] == '.tif']

        r10m_bands = band_classify(r10m_bands)
        r20m_bands = band_classify(r20m_bands)
        r60m_bands = band_classify(r60m_bands)
        # 过滤下tif文件


        os.makedirs(out_path + '/' + date + '/R10m')
        os.makedirs(out_path + '/' + date + '/R20m')
        os.makedirs(out_path + '/' + date + '/R60m')

        for key,bands in r10m_bands.items():
            print('10m---区域代码',key)
            print('10m---组合波段',bands)
            # 输出命名
            out_bandpath = out_path+'/'+date+'/R10m/'+key+'.tif'
            # 第一个band用来定义坐标系
            ds = gdal.Open(r10m_dir+'/'+bands[0])
            ds_width = ds.RasterXSize
            ds_height = ds.RasterYSize
            ds_geo = ds.GetGeoTransform()
            ds_prj = ds.GetProjection()
            driver = gdal.GetDriverByName('GTiff')
            ds_result = driver.Create(out_bandpath, ds_width, ds_height, bands=len(bands), eType=gdal.GDT_UInt16)
            ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
            ds_result.SetProjection(ds_prj)  # 导入投影信息
            del ds
            i = 1
            for band in bands:
                band_img = gdal.Open(r10m_dir+'/'+band)
                band_data = band_img.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height)
                ds_result.GetRasterBand(i).WriteArray(band_data)
                i = i+1
            del ds_result

        for key,bands in r20m_bands.items():
            print('20m---区域代码',key)
            print('20m---组合波段',bands)
            # 输出命名
            out_bandpath = out_path+'/'+date+'/R20m/'+key+'.tif'
            # 第一个band用来定义坐标系
            ds = gdal.Open(r20m_dir+'/'+bands[0])
            ds_width = ds.RasterXSize
            ds_height = ds.RasterYSize
            ds_geo = ds.GetGeoTransform()
            ds_prj = ds.GetProjection()
            driver = gdal.GetDriverByName('GTiff')
            ds_result = driver.Create(out_bandpath, ds_width, ds_height, bands=len(bands), eType=gdal.GDT_UInt16)
            ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
            ds_result.SetProjection(ds_prj)  # 导入投影信息
            del ds
            i = 1
            for band in bands:
                band_img = gdal.Open(r20m_dir+'/'+band)
                band_data = band_img.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height)
                ds_result.GetRasterBand(i).WriteArray(band_data)
                i = i+1
            del ds_result

        for key,bands in r60m_bands.items():
            print('60m---区域代码',key)
            print('60m---组合波段',bands)
            # 输出命名
            out_bandpath = out_path+'/'+date+'/R60m/'+key+'.tif'
            # 第一个band用来定义坐标系
            ds = gdal.Open(r60m_dir+'/'+bands[0])
            ds_width = ds.RasterXSize
            ds_height = ds.RasterYSize
            ds_geo = ds.GetGeoTransform()
            ds_prj = ds.GetProjection()
            driver = gdal.GetDriverByName('GTiff')
            ds_result = driver.Create(out_bandpath, ds_width, ds_height, bands=len(bands), eType=gdal.GDT_UInt16)
            ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
            ds_result.SetProjection(ds_prj)  # 导入投影信息
            del ds
            i = 1
            for band in bands:
                band_img = gdal.Open(r60m_dir+'/'+band)
                band_data = band_img.GetRasterBand(1).ReadAsArray(0, 0, ds_width, ds_height)
                ds_result.GetRasterBand(i).WriteArray(band_data)
                i = i+1
            del ds_result

if __name__ == '__main__':
    input_path = 'F:/PoyangLake/poyanghu_data/test/tif'
    out_path = 'F:/PoyangLake/poyanghu_data/test/bandcombination'
    bandcombination(input_path, out_path)

批量拼接/裁剪，按日期、分辨率存放