本教程根据卫星图像的时间序列构建无云影像图。我们将看到以下内容：
-
查找地球上特定点的时间序列图像
-
将这些图像堆叠成一个数组
-
通过取中值计算无云马赛克
-
可视化结果
此示例使用Sentinel-2 Level-2A数据。这里使用的技术同样适用于其他遥感数据集。
先看以下数据介绍：
我们可以通过在线的影像查看sentinel数据的云量状态：
Planetary Computer
在前哨2方案提供以10m-60M分辨率和大约5天，重访时间十三谱带全球图像。该数据集代表了从 2016 年到现在的全球 Sentinel-2 档案，使用Sen2Cor处理为 L2A（大气层底部）并转换为云优化的 GeoTIFF格式。
STAC系列：https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-2-l2a
供应商：ESA （生产商、许可方）| Esri （处理器）| 微软 （主机、处理器）
执照:哥白尼哨兵数据条款
波段信息
数据属性：
导入所需要的库：
import numpy as np
import xarray as xr
import rasterio.features
import stackstac
import pystac_client
import planetary_computer
import xrspatial.multispectral as ms
from dask_gateway import GatewayCluster创建 Dask 集群
我们将处理大量数据。为了减少执行时间，我们将使用 Dask 集群并行执行计算，自适应扩展以根据需要添加和删除工作人员。有关使用 Dask的更多信息，请参阅使用 Dask 扩展。
cluster = GatewayCluster() # Creates the Dask Scheduler. Might take a minute.
client = cluster.get_client()
cluster.adapt(minimum=4, maximum=24)
print(cluster.dashboard_link)https://pccompute.westeurope.cloudapp.azure.com/compute/services/dask-gateway/clusters/prod.bf38bb60fa564246b994e51950921e92/status
选择研究区
在此示例中，我们将感兴趣的区域定义为 GeoJSON 对象。它在华盛顿州雷德蒙德附近。
area_of_interest = {
"type": "Polygon",
"coordinates": [
[
[-122.27508544921875, 47.54687159892238],
[-121.96128845214844, 47.54687159892238],
[-121.96128845214844, 47.745787772920934],
[-122.27508544921875, 47.745787772920934],
[-122.27508544921875, 47.54687159892238],
]
],
}
bbox = rasterio.features.bounds(area_of_interest)用pystac_client我们可以在行星计算机的 STAC 端点中搜索与我们的查询参数匹配的项目。
#通过代理来打开，以后这个代码经常要用，因为要获取此平台的数据
stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
#设置筛选条件
search = stac.search(
bbox=bbox,
datetime="2016-01-01/2020-12-31",
collections=["sentinel-2-l2a"],
#这里需要我们进行数据集自己输入，因为官方暂时没有给出明确的数据导入名称的规范格式，相对于GEE有一点差
limit=500, # 以500个为一个批次获取限制
query={"eo:cloud_cover": {"lt": 25}},#这里设置一个云量筛选，在官方给出的文档中没有明确数据属性中可以选择的地方，所以这里eo:cloud_cover还没有明确的给出，我们只能参考，期待官方给出相应的详细的数据介绍
)
#获取筛选后的影像打印
items = list(search.get_items())
print(len(items))这里结果输出显示为138，即符合条目的影像总量，为了创建我们的无云影像，我们将使用stackstac将数据加载到xarray DataArray 中，然后将图像的时间序列减少到单个图像。
signed_items = [planetary_computer.sign(item).to_dict() for item in items]#这里是获取
data = (
stackstac.stack(
signed_items,
assets=["B04", "B03", "B02"], # 获取三个波段red, green, blue
chunksize=4096,
resolution=100,#分辨率100米
)
.where(lambda x: x > 0, other=np.nan) # sentinel-2 uses 0 as nodata
.assign_coords(band=lambda x: x.common_name.rename("band")) # use common names
)
data
由于与我们的查询匹配的数据不会太大，我们可以将其保存在分布式内存中。一旦进入内存，后续操作就会快得多。
data = data.persist()中位数影像合成
使用普通的 xarray 操作，我们可以计算时间维度上的中值。在云是瞬态的假设下，合成不应包含（许多）云，因为它们不应该是许多图像上那个点的中值像素值。
median = data.median(dim="time").compute()为了可视化数据，我们将使用 xarray-spatial 的true_color方法转换为红/绿/蓝值。
image = ms.true_color(*median) # 通过真彩色影像来展示#导入我们的画图软件
import matplotlib.pyplot as plt
#设置图像大小
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_axis_off()
#展示图像
image.plot.imshow(ax=ax);
每月综合图的展示
现在假设我们不想组合一年中不同时间段的图像（例如，我们可能不想组合通常包含雪的 1 月图像和 7 月图像）。再次使用标准的 xarray 语法，我们可以通过按月分组然后取中位数来创建一组每月组合。
monthly = data.groupby("time.month").median().compute()让我们将这些数组中的每一个影像转换为真彩色图像并将结果绘制为有序的图像。
images = [ms.true_color(*x) for x in monthly]
images = xr.concat(images, dim="time")
g = images.plot.imshow(x="x", y="y", rgb="band", col="time", col_wrap=3, figsize=(6, 8))
for ax in g.axes.flat:
ax.set_axis_off()
plt.tight_layout()最后的结果：
用python的好处就是可以直接利用matplotlib.pyplot完成图的创作。
