使用 Python 查询 Sentinel-1 数据，绘制查询结果范围并下载数据

from gma import io
import requests

xian = io.ReadVector('T2024年初县级.shp')

# 从县级里筛选出 宜阳县，也可以根据需求自行筛选或打开其他矢量文件
ly = xian.Select('"县级" = "宜阳县"')

## 获取筛选图层的范围，生成筛选范围 wkt 字符串
area_ft = ly.Extent().GetFeature()
ploygon = area_ft.ExportGeom()

## 配置查询参数
filters = ["Collection/Name eq 'SENTINEL-1'", # 数据集名称
           f"OData.CSC.Intersects(area=geography'SRID=4326;{ploygon}')",  # 查询范围，我们要找的宜阳县边界，必须为 WGS 84（EPSG:4326）坐标系
           "Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' and att/OData.CSC.StringAttribute/Value eq 'IW_SLC__1S')", # 查询产品 IW_SLC__1S 
           "ContentDate/Start gt 2025-09-01T00:00:00.000Z", # 数据开始时间
           "ContentDate/End lt 2026-01-01T00:00:00.000Z", # 数据结束时间
          ]
query_res = requests.get("https://catalogue.dataspace.copernicus.eu/odata/v1/Products?"
                         f"$filter={' and '.join(filters)}&"
                         "$orderby=ContentDate/Start&"# 按开始时间排序
                         "$top=1000").json()['value']

from gma.map import inres, plot
import numpy as np

## 读取省、市图层
shi = io.ReadVector('T2024年初地级.shp')
sheng = io.ReadVector('T2024年初省级.shp')

## 将影像范围 wkt 字符串转换为 gma 要素(Feature)
fts = [io.CreateFeature(d['Footprint'].split(';')[-1]) for d in query_res]

## 根据查询的得到的数据范围，设置一个绘图范围
bous = np.array([ft.Boundary for ft in fts])
extent = bous[:, 0].min() - 2, bous[:, 1].min() - 1, bous[:, 2].max() + 2, bous[:, 3].max() + 1

## 绘制查询结果范围
MapF = plot.MapFrame(BaseMapProj = 3857, Extent = extent)

### 添加所有查询结果
for ft in fts:
    mft = MapF.AddFeature(ft, FaceColor = 'lightgreen', alpha = 0.1)

### 添加省市县图层   
mly1 = MapF.AddLayer(xian, FaceColor = 'none', LineColor = 'lightgray', LineWidth = 0.1)
mly2 = MapF.AddLayer(shi, FaceColor = 'none', LineColor = 'darkgray', LineWidth = 0.2)
mly3 = MapF.AddLayer(sheng, FaceColor = 'none', LineColor = 'black', LineWidth = 0.3)

mly4 = MapF.AddLayer(ly, FaceColor = 'none', LineColor = 'blue', LineWidth = 0.3)

### 添加查询输入图层
mly4.AddLabel("宜阳县", Font = 'SimSun', FontSize = 5, FontColor = 'blue')

### 添加经纬网和地图框
lines = MapF.AddGridLines(LONRange = (100, 130, 2), LATRange = (26, 46, 2))
frame = MapF.SetFrame()

## 根据用户账户获取 token
def get_access_token(username, password):
    token_https = 'https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token'
    data = {"username":username,
            "password":password, 
            "grant_type":"password", 
            "client_id":"cdse-public"}
    
    r = requests.post(token_https, data = data)
    r.raise_for_status()
    return r.json()['access_token']

import tqdm

## 下载数据
### 下载地址的基础 url 模版
url = "https://download.dataspace.copernicus.eu/odata/v1/Products(%s)/$value"
out_path = r"D:\下载"# 文件保存路径
chunk_size = 8192# 下载时分块大小

for i, d in enumerate(query_res):

    Id = d["Id"]
    name = d["Name"]

    out_file = f"{out_path}\\{name[:-5]}.zip"
    if os.path.exists(out_file):
        print(f'跳过已存在的下载（{i + 1}/{len(query_res)}）：{out_file}')
        continue

    ## token会失效，每次获取一个避免失效（反正数据很大每一个数据也要下载很长时间）
    access_token = get_access_token("你的用户名", "你的密码")
    
    headers = {"Authorization": f"Bearer {access_token}"}
    session = requests.Session()
    session.headers.update(headers)
    
    data_url = url % Id
    response = session.get(data_url, stream=True)
    ##

    if response.status_code == 200:
        ## 创建进度条
        print(f'当前下载（{i + 1}/{len(query_res)}）：{out_file}')
        total_size = int(response.headers.get('content-length', 0)) # 数据总大小
        total = int(np.ceil(total_size / chunk_size)) # 分块数量
    
        with open(out_file, "wb") as file: # 分块下载
             for chunk in tqdm.tqdm(response.iter_content(chunk_size = chunk_size), total = total):
                 if chunk:
                     file.write(chunk)
    else:
        print(f"Error: {response.status_code}. {response.text}")