python landsatxplore下载Landsat8和Sentinel-2数据

说明：
1、通过landsatxplore进行query和download，需要在earthexplore网站注册
2、查询需要landsat8和sentinel-2的轨道文件，轨道文件与查询范围相交查询影像，landsat8的轨道文件网上较多，sentinel-2的轨道文件在另一篇博客中有提供下载地址（下载为kml，可转为shapefile文件）
3、参数：数据集名称（按规定填）、开始时间+结束时间（年-月-日格式）、云量、查询范围（shapefile文件，只查询第一个面要素）、存储文件夹（自定义，会在文件夹下按照年份再建一层文件夹，不想要可直接改掉）
4、代码逻辑：读取查询范围和轨道文件进行相交分析，得到轨道中心经纬度，加上数据集、时间、云量进行查询，查询结果的ID存储在列表中，根据ID列表一个个下载。
5、补充：虽然Sentinel-2查询输入的数据集是Sentinel2A，但是2A和2B都可查询得到。

import os,sys
from landsatxplore.api import API
from landsatxplore.earthexplorer import EarthExplorer
import geopandas as gpd
import warnings
warnings.filterwarnings(action="ignore")

username = 'username'  # 输入EE账号
password = 'password'  # 输入账号密码
# 初始化API接口获取key
api = API(username, password)

# 影像查询
def search_image(search_file,dataset,start_date,end_date,max_cloud_cover):
    # 输入要查询的边界
    data = gpd.read_file(search_file)
    # 输入轨道矢量
    if dataset.lower() == "sentinel_2a":
        grid_file = 'E:\*****\sentinel2_grid.shp'
    elif dataset.lower() == "landsat_8_c1":
        # 输入landsat轨道矢量
        grid_file = 'E:\*****\WRS2_descending.shp'
    wrs = gpd.GeoDataFrame.from_file(grid_file)
    # 查询边界覆盖的轨道中心坐标
    wrs_intersection = wrs[wrs.intersects(data.geometry[0])]
    longitude = (wrs_intersection.centroid).x.values
    latitude = (wrs_intersection.centroid).y.values
    # print(longitude, latitude)

    # 查询
    all_scene = []
    for i in range(len(latitude)):
        scenes = api.search(dataset=dataset,latitude=latitude[i],longitude=longitude[i],
        start_date=start_date, end_date=end_date, max_cloud_cover=max_cloud_cover)
        all_scene += scenes
    api.logout()
    return all_scene

# 下载影像数据
def Download_from_Landsatexplore(dataset,scene_list):
    if len(scene_list) > 0:
        # 根据ID下载影像
        for scene in scene_list:
            if dataset.lower() == "landsat_8_c1":
                output_dir = 'E:\***\Landsat8'# 输入下载路径
            elif dataset.lower() == "sentinel_2a":
                output_dir = 'E:\***\Sentinel2'  # 输入下载路径
            output_dir_demon = output_dir+'\\'+str(scene['acquisition_date'].year)
            if not os.path.isdir(output_dir_demon):
                os.mkdir(output_dir_demon)
            ee = EarthExplorer(username, password)
            print("Downloading: "+scene['display_id'])
            ee.download(identifier=scene['entity_id'], output_dir=output_dir_demon)
        ee.logout()

if __name__ == '__main__':

    """Landsat 5 TM Collection 1 Level 1-->landsat_tm_c1
        Landsat 5 TM Collection 2 Level 1-->landsat_tm_c2_l1
        Landsat 5 TM Collection 2 Level 2-->landsat_tm_c2_l2
        Landsat 7 ETM+ Collection 1 Level 1-->landsat_etm_c1
        Landsat 7 ETM+ Collection 2 Level 1-->landsat_etm_c2_l1
        Landsat 7 ETM+ Collection 2 Level 2-->landsat_etm_c2_l2
        Landsat 8 Collection 1 Level 1-->landsat_8_c1
        Landsat 8 Collection 2 Level 1-->landsat_ot_c2_l1
        Landsat 8 Collection 2 Level 2-->landsat_ot_c2_l2
        Sentinel 2A-->sentinel_2a"""
    # 输入查询条件
    dataset = sys.argv[1]#'landsat_8_c1' # 数据集
    start_date = sys.argv[2]#'2021-05-01' # 开始日期
    end_date = sys.argv[3]#'2021-05-20' # 结束日期
    cloud_cover = sys.argv[4]#15 # 云量（%）
    # 输入查询文件
    search_file = sys.argv[5]#r'E:\***\北京市.shp'  # 输入查询矢量

    # 查询数据
    search_list = search_image(search_file, dataset, start_date, end_date, cloud_cover)
    # print(search_list[0])
    # 下载数据
    Download_from_Landsatexplore(dataset, search_list)