本文為英文版的機器翻譯版本，如內容有任何歧義或不一致之處，概以英文版為準。

存取 Sentinel-2 點陣資料集合並建立地球觀測作業，以執行土地分割

此 Python 型教學課程使用 SDK for Python （Boto3） 和 Amazon SageMaker Studio Classic 筆記本。若要成功完成此示範，請確定您擁有使用 SageMaker 地理空間所需的 AWS Identity and Access Management （IAM） 許可，且 Studio Classic. SageMaker geospatial 需要您擁有可存取 Studio Classic 的使用者、群組或角色。您還必須在其信任政策sagemaker-geospatial.amazonaws.com中具有指定 SageMaker 地理空間服務主體的 SageMaker 執行角色。

若要進一步了解這些需求，請參閱SageMaker 地理空間IAM角色 。

本教學課程示範如何使用 SageMaker 地理空間API來完成下列任務：

用 list_raster_data_collections 來尋找可用的資料集合

SageMaker 地理空間支援多個光柵資料收集。若要進一步了解可用的資料集合，請參閱資料集合。

此示範使用從 收集的衛星資料 Sentinel-2 雲端最佳化的地理TIFF衛星。這些衛星每五天提供覆蓋地球全體陸地表面的資料。除了收集地球的地面影像外，Sentinel-2 衛星還會收集各種不同光譜帶的資料。

若要搜尋關注區域 （AOI）ARN，您需要與 Sentinel-2 衛星資料相關聯的 。若要在您的 ARNs中找到可用的資料收集及其關聯 AWS 區域，請使用 list_raster_data_collectionsAPI操作。

由於回應可以分頁，因此必須使用 get_paginator 作業傳回所有相關資料：

import boto3
import sagemaker
import sagemaker_geospatial_map
import json 

## SageMaker Geospatial  is currently only avaialable in US-WEST-2  
session = boto3.Session(region_name='us-west-2')
execution_role = sagemaker.get_execution_role()

## Creates a SageMaker Geospatial client instance 
geospatial_client = session.client(service_name="sagemaker-geospatial")

# Creates a resusable Paginator for the list_raster_data_collections API operation 
paginator = geospatial_client.get_paginator("list_raster_data_collections")

# Create a PageIterator from the paginator class
page_iterator = paginator.paginate()

# Use the iterator to iterate throught the results of list_raster_data_collections
results = []
for page in page_iterator:
    results.append(page['RasterDataCollectionSummaries'])

print(results)

這是來自 list_raster_data_collectionsAPI操作的範例JSON回應。它被截斷，只包含資料收集 （Sentinel-2） 用於他的程式碼範例。如需有關特定點陣式資料集合的更多詳細資訊，請使用 get_raster_data_collection：

{
    "Arn": "arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8",
    "Description": "Sentinel-2a and Sentinel-2b imagery, processed to Level 2A (Surface Reflectance) and converted to Cloud-Optimized GeoTIFFs",
    "DescriptionPageUrl": "https://registry.opendata.aws/sentinel-2-l2a-cogs",
    "Name": "Sentinel 2 L2A COGs",
    "SupportedFilters": [
        {
            "Maximum": 100,
            "Minimum": 0,
            "Name": "EoCloudCover",
            "Type": "number"
        },
        {
            "Maximum": 90,
            "Minimum": 0,
            "Name": "ViewOffNadir",
            "Type": "number"
        },
        {
            "Name": "Platform",
            "Type": "string"
        }
    ],
    "Tags": {},
    "Type": "PUBLIC"
}

執行先前的程式碼範例後，您會取得 Sentinel-2 光柵資料收集ARN的 arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8。在下一節中，您可以使用 查詢 Sentinel-2 search_raster_data_collection 資料收集API。

搜尋 Sentinel-2 使用 的光柵資料收集 search_raster_data_collection

在上一節中，您曾list_raster_data_collections用來取得 ARN 的 Sentinel-2 資料收集。現在，您可以使用 ARN來搜尋指定關注區域的資料收集 （AOI）、時間範圍、屬性和可用的 UV 頻帶。

若要呼叫 search_raster_data_collectionAPI，您必須在 Python RasterDataCollectionQuery 參數的索引。此範例使用 AreaOfInterest、TimeRangeFilter、PropertyFilters 跟 BandFilter。為了方便起見，您可以使用變數 search_rdc_query 指定 Python 字典來保留搜尋查詢參數：

search_rdc_query = {
    "AreaOfInterest": {
        "AreaOfInterestGeometry": {
            "PolygonGeometry": {
                "Coordinates": [
                    [
                        # coordinates are input as longitute followed by latitude 
                        [-114.529, 36.142],
                        [-114.373, 36.142],
                        [-114.373, 36.411],
                        [-114.529, 36.411],
                        [-114.529, 36.142],
                    ]
                ]
            }
        }
    },
    "TimeRangeFilter": {
        "StartTime": "2022-01-01T00:00:00Z",
        "EndTime": "2022-07-10T23:59:59Z"
    },
    "PropertyFilters": {
        "Properties": [
            {
                "Property": {
                    "EoCloudCover": {
                        "LowerBound": 0,
                        "UpperBound": 1
                    }
                }
            }
        ],
        "LogicalOperator": "AND"
    },
    "BandFilter": [
        "visual"
    ]
}

在此範例中，您會查詢包含猶他州 Lake Mead 的一個 AreaOfInterest。此外 Sentinel-2 支援多個類型的影像帶。要測量水面的變化，你只需要 visual 帶。

建立查詢參數後，您可以使用 search_raster_data_collection API 提出請求。

下列程式碼範例會實作 search_raster_data_collectionAPI請求。這API不支援使用 get_paginator 分頁API。為了確保已收集完整的API回應，程式碼範例會使用while迴圈來檢查NextToken是否存在。然後，程式碼範例會使用 .extend()將衛星映像URLs和其他回應中繼資料附加至 items_list。

若要進一步了解 search_raster_data_collection，請參閱 Amazon 參考 SearchRasterDataCollection 中的 。 SageMaker API

search_rdc_response = sm_geo_client.search_raster_data_collection(
    Arn='arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8',
    RasterDataCollectionQuery=search_rdc_query
)


## items_list is the response from the API request. 
items_list = []

## Use the python .get() method to check that the 'NextToken' exists, if null returns None breaking the while loop 
while search_rdc_response.get('NextToken'):
    items_list.extend(search_rdc_response['Items'])
    search_rdc_response = sm_geo_client.search_raster_data_collection(
        Arn='arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8',
        RasterDataCollectionQuery=search_rdc_query, NextToken=search_rdc_response['NextToken']
    )

## Print the number of observation return based on the query
print (len(items_list))

以下是查詢的JSON回應。為保持清晰起見，已被截斷。在 Assets 鍵值對中僅會傳回請求中指定的 "BandFilter": ["visual"]：

{
    'Assets': {
        'visual': {
            'Href': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/UH/2022/6/S2A_15TUH_20220623_0_L2A/TCI.tif'
        }
    },
    'DateTime': datetime.datetime(2022, 6, 23, 17, 22, 5, 926000, tzinfo = tzlocal()),
    'Geometry': {
        'Coordinates': [
            [
                [-114.529, 36.142],
                [-114.373, 36.142],
                [-114.373, 36.411],
                [-114.529, 36.411],
                [-114.529, 36.142],
            ]
        ],
        'Type': 'Polygon'
    },
    'Id': 'S2A_15TUH_20220623_0_L2A',
    'Properties': {
        'EoCloudCover': 0.046519,
        'Platform': 'sentinel-2a'
    }
}

現在您有了查詢結果，在下一節中可以使用 matplotlib 來視覺化結果。這是為了驗證結果來自正確的地理區域。

使用 matplotlib 視覺化您的 search_raster_data_collection

在開始地球觀測任務 （EOJ） 之前，您可以使用 視覺化查詢的結果matplotlib。下列程式碼範例會從前一個程式碼範例中建立的 items_list 變數取得第一個項目 items_list[0]["Assets"]["visual"]["Href"]，並使用 matplotlib 列印影像。

# Visualize an example image.
import os
from urllib import request
import tifffile
import matplotlib.pyplot as plt

image_dir = "./images/lake_mead"
os.makedirs(image_dir, exist_ok=True)

image_dir = "./images/lake_mead"
os.makedirs(image_dir, exist_ok=True)

image_url = items_list[0]["Assets"]["visual"]["Href"]
img_id = image_url.split("/")[-2]
path_to_image = image_dir + "/" + img_id + "_TCI.tif"
response = request.urlretrieve(image_url, path_to_image)
print("Downloaded image: " + img_id)

tci = tifffile.imread(path_to_image)
plt.figure(figsize=(6, 6))
plt.imshow(tci)
plt.show()

檢查結果是否位於正確的地理區域之後，您可以在下一個步驟中啟動地球觀測任務 （EOJ）。您可以使用 EOJ來識別衛星影像中的水體，方法是使用稱為土地分割的程序。

啟動對一系列衛星影像執行土地分割的地球觀測任務 （EOJ）

SageMaker 地理空間提供多個預先訓練的模型，可用來處理來自光柵資料收集的地理空間資料。若要進一步了解可用的預先訓練模型和自訂操作，請參閱作業類型。

要計算水面積的變化，您需要識別影像中的哪些像素對應到水。地蓋分割是 支援的語意分割模型start_earth_observation_jobAPI。語意分割模型會將標籤與每個影像中的每個像素產生關聯。在結果中，每個像素都會根據模型的類別地圖指派一個標籤。以下是土地分割模型的類別地圖：

{
    0: "No_data",
    1: "Saturated_or_defective",
    2: "Dark_area_pixels",
    3: "Cloud_shadows",
    4: "Vegetation",
    5: "Not_vegetated",
    6: "Water",
    7: "Unclassified",
    8: "Cloud_medium_probability",
    9: "Cloud_high_probability",
    10: "Thin_cirrus",
    11: "Snow_ice"
}

若要啟動接地觀測任務，請使用 start_earth_observation_job API。當您提交請求時，您必須指定下列項目：

InputConfig 是 Python 字典。使用下列變數 eoj_input_config 來保留搜尋查詢參數。當您提出start_earth_observation_jobAPI請求時，請使用此變數。w.

# Perform land cover segmentation on images returned from the Sentinel-2 dataset.
eoj_input_config = {
    "RasterDataCollectionQuery": {
        "RasterDataCollectionArn": "arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8",
        "AreaOfInterest": {
            "AreaOfInterestGeometry": {
                "PolygonGeometry": {
                    "Coordinates":[
                        [
                            [-114.529, 36.142],
                            [-114.373, 36.142],
                            [-114.373, 36.411],
                            [-114.529, 36.411],
                            [-114.529, 36.142],
                        ]
                    ]
                }
            }
        },
        "TimeRangeFilter": {
            "StartTime": "2021-01-01T00:00:00Z",
            "EndTime": "2022-07-10T23:59:59Z",
        },
        "PropertyFilters": {
            "Properties": [{"Property": {"EoCloudCover": {"LowerBound": 0, "UpperBound": 1}}}],
            "LogicalOperator": "AND",
        },
    }
}

JobConfig 是 Python 用來指定要在資料上執行EOJ的操作的字典：

eoj_config = {"LandCoverSegmentationConfig": {}}

透過現在指定的字典元素，您可以使用下列程式碼範例提交start_earth_observation_jobAPI請求：

# Gets the execution role arn associated with current notebook instance 
execution_role_arn = sagemaker.get_execution_role()

# Starts an earth observation job
response = sm_geo_client.start_earth_observation_job(
    Name="lake-mead-landcover",
    InputConfig=eoj_input_config,
    JobConfig=eoj_config,
    ExecutionRoleArn=execution_role_arn,
)
            
print(response)

開始接地觀察任務會傳回 ARN和其他中繼資料。

若要取得所有進行中和目前地球觀測任務的清單，請使用 list_earth_observation_jobs API。若要監控單一地球觀測任務的狀態，請使用 get_earth_observation_job API。若要提出此請求，請使用提交EOJ請求後ARN建立的 。若要進一步了解，請參閱 Amazon 參考 GetEarthObservationJob 中的 。 SageMaker API

若要尋找與您的 ARNs相關聯的 ，EOJs請使用 list_earth_observation_jobsAPI操作。若要進一步了解，請參閱 Amazon 參考 ListEarthObservationJobs 中的 。 SageMaker API

# List all jobs in the account
sg_client.list_earth_observation_jobs()["EarthObservationJobSummaries"]

以下是JSON回應範例：

{
    'Arn': 'arn:aws:sagemaker-geospatial:us-west-2:111122223333:earth-observation-job/futg3vuq935t',
    'CreationTime': datetime.datetime(2023, 10, 19, 4, 33, 54, 21481, tzinfo = tzlocal()),
    'DurationInSeconds': 3493,
    'Name': 'lake-mead-landcover',
    'OperationType': 'LAND_COVER_SEGMENTATION',
    'Status': 'COMPLETED',
    'Tags': {}
}, {
    'Arn': 'arn:aws:sagemaker-geospatial:us-west-2:111122223333:earth-observation-job/wu8j9x42zw3d',
    'CreationTime': datetime.datetime(2023, 10, 20, 0, 3, 27, 270920, tzinfo = tzlocal()),
    'DurationInSeconds': 1,
    'Name': 'mt-shasta-landcover',
    'OperationType': 'LAND_COVER_SEGMENTATION',
    'Status': 'INITIALIZING',
    'Tags': {}
}

在您EOJ任務的狀態變更為 之後COMPLETED，請繼續下一個區段以計算 Lake 中的變更 Mead's 表面積。

在 Lake 中計算變更 Mead 表面積

若要計算 Lake Mead 表面積的變化，請先使用 將 的結果匯出EOJ至 Amazon S3export_earth_observation_job：

sagemaker_session = sagemaker.Session()
s3_bucket_name = sagemaker_session.default_bucket()  # Replace with your own bucket if needed
s3_bucket = session.resource("s3").Bucket(s3_bucket_name)
prefix = "export-lake-mead-eoj"  # Replace with the S3 prefix desired
export_bucket_and_key = f"s3://{s3_bucket_name}/{prefix}/"

eoj_output_config = {"S3Data": {"S3Uri": export_bucket_and_key}}
export_response = sm_geo_client.export_earth_observation_job(
    Arn="arn:aws:sagemaker-geospatial:us-west-2:111122223333:earth-observation-job/7xgwzijebynp",
    ExecutionRoleArn=execution_role_arn,
    OutputConfig=eoj_output_config,
    ExportSourceImages=False,
)

若要檢視匯出任務的狀態，請使用 get_earth_observation_job：

export_job_details = sm_geo_client.get_earth_observation_job(Arn=export_response["Arn"])

若要計算 Lake Mead 水位的變化，請將土地覆蓋遮罩下載至本機 SageMaker 筆記本執行個體，並從我們先前的查詢下載來源映像。在土地分割模型的分類地圖中，水的等級索引為 6。

從 中擷取水遮罩 Sentinel-2 映像，請遵循下列步驟。首先計算影像中標記為水 (類別索引 6) 的像素數目。其次，將計數乘以每個像素覆蓋的區域。頻帶在其空間解析度可能有所不同。針對土地覆蓋分割模型，所有頻帶都將降縮小抽樣取樣至等於 60 公尺的空間解析度。

import os
from glob import glob
import cv2
import numpy as np
import tifffile
import matplotlib.pyplot as plt
from urllib.parse import urlparse
from botocore import UNSIGNED
from botocore.config import Config

# Download land cover masks
mask_dir = "./masks/lake_mead"
os.makedirs(mask_dir, exist_ok=True)
image_paths = []
for s3_object in s3_bucket.objects.filter(Prefix=prefix).all():
    path, filename = os.path.split(s3_object.key)
    if "output" in path:
        mask_name = mask_dir + "/" + filename
        s3_bucket.download_file(s3_object.key, mask_name)
        print("Downloaded mask: " + mask_name)

# Download source images for visualization
for tci_url in tci_urls:
    url_parts = urlparse(tci_url)
    img_id = url_parts.path.split("/")[-2]
    tci_download_path = image_dir + "/" + img_id + "_TCI.tif"
    cogs_bucket = session.resource(
        "s3", config=Config(signature_version=UNSIGNED, region_name="us-west-2")
    ).Bucket(url_parts.hostname.split(".")[0])
    cogs_bucket.download_file(url_parts.path[1:], tci_download_path)
    print("Downloaded image: " + img_id)

print("Downloads complete.")

image_files = glob("images/lake_mead/*.tif")
mask_files = glob("masks/lake_mead/*.tif")
image_files.sort(key=lambda x: x.split("SQA_")[1])
mask_files.sort(key=lambda x: x.split("SQA_")[1])
overlay_dir = "./masks/lake_mead_overlay"
os.makedirs(overlay_dir, exist_ok=True)
lake_areas = []
mask_dates = []

for image_file, mask_file in zip(image_files, mask_files):
    image_id = image_file.split("/")[-1].split("_TCI")[0]
    mask_id = mask_file.split("/")[-1].split(".tif")[0]
    mask_date = mask_id.split("_")[2]
    mask_dates.append(mask_date)
    assert image_id == mask_id
    image = tifffile.imread(image_file)
    image_ds = cv2.resize(image, (1830, 1830), interpolation=cv2.INTER_LINEAR)
    mask = tifffile.imread(mask_file)
    water_mask = np.isin(mask, [6]).astype(np.uint8)  # water has a class index 6
    lake_mask = water_mask[1000:, :1100]
    lake_area = lake_mask.sum() * 60 * 60 / (1000 * 1000)  # calculate the surface area
    lake_areas.append(lake_area)
    contour, _ = cv2.findContours(water_mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    combined = cv2.drawContours(image_ds, contour, -1, (255, 0, 0), 4)
    lake_crop = combined[1000:, :1100]
    cv2.putText(lake_crop, f"{mask_date}", (10,50), cv2.FONT_HERSHEY_SIMPLEX, 1.5, (0, 0, 0), 3, cv2.LINE_AA)
    cv2.putText(lake_crop, f"{lake_area} [sq km]", (10,100), cv2.FONT_HERSHEY_SIMPLEX, 1.5, (0, 0, 0), 3, cv2.LINE_AA)
    overlay_file = overlay_dir + '/' + mask_date + '.png'
    cv2.imwrite(overlay_file, cv2.cvtColor(lake_crop, cv2.COLOR_RGB2BGR))

# Plot water surface area vs. time.
plt.figure(figsize=(20,10))
plt.title('Lake Mead surface area for the 2021.02 - 2022.07 period.', fontsize=20)
plt.xticks(rotation=45)
plt.ylabel('Water surface area [sq km]', fontsize=14)
plt.plot(mask_dates, lake_areas, marker='o')
plt.grid('on')
plt.ylim(240, 320)
for i, v in enumerate(lake_areas):
    plt.text(i, v+2, "%d" %v, ha='center')
plt.show()

您可以使用 matplotlib 以圖表視覺化結果。圖表顯示 Lake 的表面積 Mead 從 2021 年 1 月至 2022 年 7 月期間減少。