专栏名称: GEE遥感训练营
专注GEE遥感算法，包括遥感影像下载、遥感影像制图、遥感GIS空间分析、遥感生态评价、遥感影像融合、遥感去干扰等多元遥感云计算
目录
相关文章推荐
51好读  ›  专栏  ›  GEE遥感训练营

Microsoft Planetary Computer环境配置和影像本地可视化

GEE遥感训练营  · 公众号  ·  · 2024-03-28 21:30

正文

请到「今天看啥」查看全文


Microsoft Planetary  Computer环境配置和影像本地可视化

引文

从2015年的第一款遥感云计算GEE面世以来，截止到2024年的今天，遥感云计算平台层出不穷（Microsoft Planetary  Computer、ai earth、星图地球智脑引擎、PIE-engine等）。

GEE是毫无疑问的行业老大，但是 哪一个平台对于中国用户 才会是最强大的平台了？我觉得是微软行星计算(Microsoft Planetary  Computer,以下简称MPC)。

Microsoft Planetary  Computer是微软开发的一个遥感云计算平台，于2021年底上线，旨在提供用户访问PB级地理空间数据和计算资源，用于地球观测和环境分析。相比其他平台，MPC有以下优点：

（1） 国内直连 ，无需借助梯子（优于GEE）；

（2）由Azure提供计算云服务， 算力更强 （优于PIE-engine）；

（3） 开放程度高 ，提供一系列的API接口以及开放的数据源，适合遥感工程类项目的开发（优于其他所有遥感云计算平台）

现在我想使用MPC去筛选哨兵影像，并在本地电脑进行显示。应该怎么做呢？注册账号、 配置环境 、 影像可视化 。

MPC平台注册

进入MPC官网https://planetarycomputer.microsoft.com/，点击 Request access

按照要求填写相关信息:

等待微软通过，我是在第二天就受到了申请通过通知。通过后，即可使用MPC平台。

MPC本地环境搭建

MPC的相关包需要使用python3.8以上环境，这里使用conda安装python3.10： 

conda create -n MPCevir python=3.10

然后再安装需要的包： 

conda install -c conda-forge pystac-client planetary-computer odc-stac matplotlib

最后安装编辑器环境jupyter: 

conda install -c anaconda jupyter

全部安装完毕后，新建一个python文件，填入下面的信息测试是否能运行成功。 

import pystac_client
import planetary_computer
import odc.stac
from pystac.extensions.eo import EOExtension as eo

上面代码能运行，环境安装成功

配置API密钥

首先登录MPC官网，在hub界面选中tkoen，网址为：https://pccompute.westeurope.cloudapp.azure.com/compute/hub/token，生成一个新的token。

在电脑终端设置环境，输入setx PC_SDK_SUBSCRIPTION_KEY "你的token"，待token保存成功，即本地可以调用MPC

影像可视化

从 pystac.extensions.eo 中导入 EOExtension 类并将其重命名为 eo 。导入 pystac_client 和 planetary_computer 。 

from pystac.extensions.eo import EOExtension as eo
import pystac_client
import planetary_computer

使用 pystac_client.Client.open() 方法打开STAC API，并指定修改器为 planetary_computer.sign_inplace 

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

使用字典定义一个名为 area_of_interest 的多边形区域，其中包含一系列坐标点: 

area_of_interest = {
    "type": "Polygon",
    "coordinates": [
        [
            [116.0737075805664, 39.63288230836462],
            [116.0737075805664, 40.19331309638433],
            [116.78659818115234, 40.19331309638433],
            [116.78659818115234, 39.63288230836462],
            [116.0737075805664, 39.63288230836462],
        ]
    ],
}

定义一个名为 time_of_interest 的时间范围，以过滤影像 

time_of_interest = "2024-01-01/2024-03-27"

使用数据集、区域和时间筛选影像: 

search = catalog.search(
    collections=["sentinel-2-l2a"],
    intersects=area_of_interest,
    datetime=time_of_interest,
    query={"eo:cloud_cover": {"lt": 10}},
)

# 获取搜索结果中云量最少的影像
least_cloudy_item = min(search.item_collection(), key=lambda item: eo.ext(item).cloud_cover)

数据集展示： 

# 获取影像的可视化URL
asset_href = least_cloudy_item.assets["visual"].href

# 使用rasterio读取影像数据
with rasterio.open(asset_href) as ds:
    # 获取影像数据的窗口
    aoi_bounds = features.bounds(area_of_interest)
    warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
    aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
    # 读取影像数据
    band_data = ds.read(window=aoi_window)

# 将数据从band-interleave格式转换为pixel-interleave格式，并下采样以便绘图
img = Image.fromarray(np.transpose(band_data, axes=[1, 2, 0]))
w = img.size[0]
h = img.size[1]
aspect = w / h
target_w = 800
target_h = int(target_w / aspect)
img = img.resize((target_w, target_h), Image.Resampling.BILINEAR)

# 使用matplotlib显示影像
plt.imshow(img)
plt.show()
效果展示

完整代码

sentinel加载影像的完整代码 

# 导入必要的库
from pystac.extensions.eo import EOExtension as eo
import pystac_client
import planetary_computer
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# 设置Planetary Computer STAC API的URL和修改器
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

# 定义感兴趣的时间范围
time_of_interest = "2024-01-01/2024-03-27"

# 定义感兴趣的地理区域（北京）
area_of_interest = {
    "type": "Polygon",
    "coordinates": [
        [
            [116.207275390625, 39.74009990142287],
            [116.598876953125, 39.74009990142287],
            [116.598876953125, 40.01969957176073],
            [116.207275390625, 40.01969957176073],
            [116.207275390625, 39.74009990142287],
        ]
    ],
}

# 搜索Sentinel-2卫星影像
search = catalog.search(
    collections=["sentinel-2-l2a"],
    intersects=area_of_interest,
    datetime=time_of_interest,
    query={"eo:cloud_cover": {"lt": 10}},
)

# 获取搜索结果中云量最少的影像
least_cloudy_item = min(search.item_collection(), key=lambda item: eo.ext(item).cloud_cover)

# 获取影像的可视化URL
asset_href = least_cloudy_item.assets["visual"].href

# 使用rasterio读取影像数据
with rasterio.open(asset_href) as ds:
    # 获取影像数据的窗口
    aoi_bounds = features.bounds(area_of_interest)
    warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
    aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
    # 读取影像数据
    band_data = ds.read(window=aoi_window)

# 将数据从band-interleave格式转换为pixel-interleave格式，并下采样以便绘图
img = Image.fromarray(np.transpose(band_data, axes=[1, 2, 0]))
w = img.size[0]
h = img.size[1]
aspect = w / h
target_w = 800
target_h = int(target_w / aspect)
img = img.resize((target_w, target_h), Image.Resampling.BILINEAR)

# 使用matplotlib显示影像
plt.imshow(img)
plt.show()

Landsat加载影像的完整代码 

import pystac_client
import planetary_computer
import odc.stac
import matplotlib.pyplot as plt
import matplotlib as mpl
# 设置字体为 SimHei
mpl.rcParams['font.sans-serif'] = ['SimHei']

# 从 pystac.extensions.eo 模块中导入 EOExtension 类，并将其重命名为 eo
from pystac.extensions.eo import EOExtension as eo

# 打开 Planetary Computer 的 STAC API 并创建一个客户端对象
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

# 定义一个感兴趣的地理边界框（bounding box），格式为 [西经，北纬，东经，]
bbox_of_interest = [113,39,114,49]

# 定义感兴趣的时间范围，格式为 "起始时间/结束时间"
time_of_interest = "2024-01-01/2024-03-30"

# 使用 catalog.search() 函数在 STAC API 中搜索符合条件的数据，包括感兴趣的地理边界框、时间范围和云覆盖率
search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox_of_interest,
    datetime=time_of_interest,
    query={"eo:cloud_cover": {"lt": 10}},
)

# 将搜索结果转换为一个项目（item）集合
items = search.item_collection()

# 打印搜索结果的项目数量
print(f"该地区指定时间范围内 共有{len(items)} 景影像")

# 选择云覆盖率最低的项目
selected_item = min(items, key=lambda item: eo.ext(item).cloud_cover)

# 打印所选项目的信息，包括项目 ID、日期和云覆盖率
print(
    f"从 {selected_item.id} 选中 {selected_item.datetime.date()}"
    + f" 云覆盖率为 {selected_item.properties['eo:cloud_cover']}% "
)

# 定义感兴趣的波段
bands_of_interest = ["nir08", "red", "green", "blue"]

# 使用 odc.stac.stac_load() 函数加载所选项目的数据，并指定感兴趣的波段和地理边界框
data = odc.stac.stac_load(
    [selected_item], bands=bands_of_interest, bbox=bbox_of_interest
).isel(time=0)

# 创建一个大小为 (10, 10) 的子图，并获取子图的轴对象
fig, ax = plt.subplots(figsize=(10, 10))

# 使用 data[["red", "green", "blue"]].to_array().plot.imshow() 函数绘制真彩色图像
data[["red", "green", "blue"]].to_array().plot.imshow(robust=True, ax=ax)

# 设置子图的标题
ax.set_title("2024年Landsat真彩色影像")

# 显示图像
plt.show()

# 使用 data[["red", "green", "blue"]].to_array().plot.imshow() 函数绘制假彩色图像
data[["nir08","red", "green"]].to_array().plot.imshow(robust=True, ax=ax)

# 设置子图的标题
ax.set_title("2024年Landsat假彩色影像")

# 显示图像
plt.show()

写在最后

MPC是优秀的、开放获取的遥感云计算平台，大家可以多试试。也可以把MPC可以视作添加了云计算的亚马逊开放平台（Open Data on AWS），稍微调整一下，可以和遥感的相关工程项目相结合，获取及时的影像和遥感产品 。

参考

MPC官网.https://planetarycomputer.microsoft.com/

官方代码示例.https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a#Example-Notebook








请到「今天看啥」查看全文


推荐文章
养花大全  ·  荐号 | 如何给花卉拍照？看看这个就知道了
8 年前
军武次位面  ·  “战利品”坦克主动防御系统，打RPG如拍苍蝇！
8 年前
史事挖掘机  ·  这个皇帝杀了15万贪官，陵墓600年无人敢盗
8 年前
教你看穿男人的心  ·  女人必懂的3个温柔交谈法，让老公从此乖乖听话 !
7 年前
PChouse家居APP  ·  1940年波尔图老公寓翻新，不变的是怀旧气质！
7 年前
Sov5搜索   ·   小百科   ·   今天看啥   ·   移动版
51好读 - 好文章就要读起来!