数据获取实战：用Python搞定Sentinel-2数据

说实话，做量化分析最头疼的事之一，就是数据获取。尤其是卫星数据，看着一堆API文档，头都大了。但别担心，今天我就带你手把手搞定Sentinel-2数据的下载和预处理。

我个人习惯用两个库：eodag 和 sentinelhub。前者适合快速搜索和下载，后者适合精细控制。你想想看，一个像淘宝搜商品，一个像直接去工厂定制，各有各的好。

为什么选Sentinel-2？

Sentinel-2是欧空局的明星卫星，10米空间分辨率，5天重访周期。说白了，就是既能看清田块细节，又能频繁更新数据。我在做大豆产量预测时，就靠它捕捉到了关键生长期的植被变化。

它的波段设计也很贴心：

波段 中心波长(nm) 分辨率(m) 主要用途
B2 (蓝) 490 10 水体识别
B3 (绿) 560 10 植被健康
B4 (红) 665 10 叶绿素吸收
B8 (近红外) 842 10 植被密度

环境准备：先把工具装好

嗯，这里要注意。别一上来就pip install，容易踩坑。我建议用conda创建独立环境：

conda create -n sentinel python=3.9
conda activate sentinel
pip install eodag sentinelhub rasterio numpy matplotlib

我曾经在Windows上装sentinelhub时遇到SSL证书问题，折腾了半天。后来发现加个--trusted-host参数就解决了。如果你也遇到类似问题，试试：

pip install sentinelhub --trusted-host pypi.org --trusted-host files.pythonhosted.org

方法一：用eodag快速上手

eodag的好处是配置简单。你只需要注册一个EarthData账号，拿到API Key就行。我个人觉得它特别适合快速验证想法。

先配置一下：

from eodag import EODataAccessGateway

dag = EODataAccessGateway()
dag.set_preferred_provider("cop_cloud")  # 选个靠谱的提供商

然后搜索数据：

# 搜索2023年7月，覆盖某个玉米产区的数据
products = dag.search(
    productType="S2_MSI_L2A",
    geom={"type": "Point", "coordinates": [116.4, 39.9]},
    start="2023-07-01",
    end="2023-07-31",
    cloudCover=20  # 云量小于20%
)

print(f"找到 {len(products)} 景影像")
for p in products[:3]:
    print(f"日期: {p.properties['startTimeFromAscendingNode']}, 云量: {p.properties['cloudCover']}%")

下载也很直接：

# 下载前3景
for product in products[:3]:
    product.download(output_dir="./sentinel_data")
小技巧： 搜索时加上cloudCover=20，能过滤掉大部分多云影像。我在做时间序列分析时，一般要求云量低于10%，否则插值误差太大。

方法二：用sentinelhub精细控制

如果你需要更精细的控制，比如只下载某个波段，或者指定投影方式，那sentinelhub是更好的选择。它需要你注册Sentinel Hub账号，拿到Instance ID。

先初始化：

from sentinelhub import (
    SHConfig, BBox, CRS, DataCollection,
    SentinelHubRequest, MimeType
)

config = SHConfig()
config.sh_client_id = "你的client_id"
config.sh_client_secret = "你的client_secret"

然后定义你要下载的区域和波段：

# 定义北京周边的一个矩形区域
bbox = BBox([116.3, 39.8, 116.5, 40.0], crs=CRS.WGS84)

# 定义要下载的波段：红、绿、蓝、近红外
evalscript = """
//VERSION=3
function setup() {
    return {
        input: ["B02", "B03", "B04", "B08"],
        output: { bands: 4 }
    };
}
function evaluatePixel(sample) {
    return [sample.B04, sample.B03, sample.B02, sample.B08];
}
"""

request = SentinelHubRequest(
    evalscript=evalscript,
    input_data=[
        SentinelHubRequest.input_data(
            data_collection=DataCollection.SENTINEL2_L2A,
            time_interval=("2023-07-15", "2023-07-16"),
            maxcc=20  # 最大云量20%
        )
    ],
    responses=[
        SentinelHubRequest.output_response("default", MimeType.TIFF)
    ],
    bbox=bbox,
    size=[512, 512],
    config=config
)

data = request.get_data()[0]
print(f"数据形状: {data.shape}")  # (512, 512, 4)
注意： sentinelhub的免费额度有限，每天大概能下载几百MB。如果你要批量下载，建议先算好预算，或者用eodag做初步筛选。

数据预处理：去云、裁剪、重采样

下载下来的数据，说白了就是一堆原始像素。直接拿去建模？那肯定不行。我踩过的坑太多了，这里给你总结三个关键步骤。

1. 去云：别让云朵骗了你的模型

云层会严重干扰植被指数计算。Sentinel-2自带一个云掩膜波段（QA60），我们可以用它来去云：

import rasterio
import numpy as np

def remove_clouds(tiff_path, output_path):
    with rasterio.open(tiff_path) as src:
        # 读取所有波段
        bands = src.read()
        # 读取云掩膜（通常在第11或12波段）
        cloud_mask = src.read(11)  # QA60波段
        
        # 云掩膜中，bit 10和bit 11分别代表云和卷云
        cloud_mask = (cloud_mask & (1 << 10)) | (cloud_mask & (1 << 11))
        cloud_mask = cloud_mask > 0
        
        # 将云区域设为NaN
        bands[:, cloud_mask] = np.nan
        
        # 写入新文件
        with rasterio.open(
            output_path, 'w',
            driver='GTiff',
            height=src.height,
            width=src.width,
            count=src.count,
            dtype=src.dtypes[0],
            crs=src.crs,
            transform=src.transform
        ) as dst:
            dst.write(bands)

# 使用
remove_clouds("S2_20230715.tiff", "S2_20230715_cloudfree.tiff")

我曾经犯过一个错误：直接用插值填充云区域。结果模型把云阴影误判为水体，预测偏差很大。后来我改用时间序列上的前后影像插值，效果才好起来。

2. 裁剪：只关注你的研究区

整景影像太大了，动辄几百MB。我们只需要研究区的那一小块：

import geopandas as gpd
from shapely.geometry import box

def clip_to_aoi(tiff_path, aoi_shapefile, output_path):
    # 读取研究区边界
    aoi = gpd.read_file(aoi_shapefile)
    
    with rasterio.open(tiff_path) as src:
        # 裁剪
        out_image, out_transform = rasterio.mask.mask(
            src, aoi.geometry, crop=True
        )
        
        # 保存
        with rasterio.open(
            output_path, 'w',
            driver='GTiff',
            height=out_image.shape[1],
            width=out_image.shape[2],
            count=src.count,
            dtype=out_image.dtype,
            crs=src.crs,
            transform=out_transform
        ) as dst:
            dst.write(out_image)

# 使用
clip_to_aoi("S2_20230715_cloudfree.tiff", "farm_boundary.shp", "S2_clipped.tiff")
避坑指南： 确保你的研究区shapefile和影像的坐标系一致。我曾经因为坐标系不匹配，裁剪出来一片空白，排查了半天才发现是CRS的问题。

3. 重采样：统一分辨率

Sentinel-2的不同波段分辨率不同（10m、20m、60m）。建模时，必须统一到同一个分辨率：

def resample_to_uniform(tiff_path, target_res=10, output_path="S2_resampled.tiff"):
    with rasterio.open(tiff_path) as src:
        # 计算目标分辨率下的尺寸
        scale = src.res[0] / target_res
        new_width = int(src.width * scale)
        new_height = int(src.height * scale)
        
        # 重采样
        data = src.read(
            out_shape=(src.count, new_height, new_width),
            resampling=rasterio.enums.Resampling.bilinear
        )
        
        # 更新变换矩阵
        transform = src.transform * src.transform.scale(
            (src.width / data.shape[-1]),
            (src.height / data.shape[-2])
        )
        
        # 保存
        with rasterio.open(
            output_path, 'w',
            driver='GTiff',
            height=new_height,
            width=new_width,
            count=src.count,
            dtype=data.dtype,
            crs=src.crs,
            transform=transform
        ) as dst:
            dst.write(data)

# 使用
resample_to_uniform("S2_clipped.tiff", target_res=10)

你想想看，如果不做重采样，10m和60m的波段混在一起建模，模型会优先学习分辨率高的波段，导致偏差。我一般统一到10m，因为植被指数（NDVI）主要依赖10m的红和近红外波段。

完整流程：从搜索到预处理

把上面所有步骤串起来，就是一个完整的pipeline：

def sentinel_pipeline(aoi_coords, start_date, end_date, output_dir):
    # 1. 搜索数据
    dag = EODataAccessGateway()
    products = dag.search(
        productType="S2_MSI_L2A",
        geom={"type": "Point", "coordinates": aoi_coords},
        start=start_date,
        end=end_date,
        cloudCover=20
    )
    
    for product in products[:5]:  # 处理前5景
        # 2. 下载
        product.download(output_dir=output_dir)
        
        # 3. 去云
        tiff_path = f"{output_dir}/{product.properties['title']}.tiff"
        cloudfree_path = tiff_path.replace(".tiff", "_cloudfree.tiff")
        remove_clouds(tiff_path, cloudfree_path)
        
        # 4. 裁剪
        clipped_path = cloudfree_path.replace("_cloudfree", "_clipped")
        clip_to_aoi(cloudfree_path, "aoi.shp", clipped_path)
        
        # 5. 重采样
        resampled_path = clipped_path.replace("_clipped", "_resampled")
        resample_to_uniform(clipped_path, target_res=10, output_path=resampled_path)
        
        print(f"处理完成: {resampled_path}")

# 运行
sentinel_pipeline(
    aoi_coords=[116.4, 39.9],
    start_date="2023-07-01",
    end_date="2023-07-31",
    output_dir="./sentinel_processed"
)

本章知识体系

下面这张图，是我自己总结的Sentinel-2数据处理流程，你保存下来，以后做项目时对照着看：

Sentinel-2数据获取与预处理流程 1. 数据搜索 eodag / sentinelhub 2. 数据下载 多线程 / 断点续传 3. 去云处理 QA60波段掩膜 4. 裁剪AOI Shapefile边界 5. 重采样 统一10m分辨率 6. 输出TIFF 可用于建模 关键参数说明 • 云量阈值：建议 ≤ 20%，时间序列分析建议 ≤ 10% • 重采样方法：双线性插值（连续数据）或最近邻（分类数据） • 输出格式：GeoTIFF，保留地理坐标信息 • 坐标系：统一为 WGS84 (EPSG:4326) 或 UTM 投影

嗯，到这里，数据获取和预处理的核心内容就讲完了。你可能会问：这些代码跑起来会不会很慢？说实话，如果数据量大，确实会慢。我一般用多线程下载，预处理时用rasterio的窗口读取，能快不少。

最后提醒一句：数据质量决定模型上限。预处理做得再精细都不为过。我见过太多人，数据没处理好就急着建模，结果花了几周调参，还不如别人花一天做数据清洗的效果好。

专注资料整理
← 上一章 📖 返回目录 下一章 →