翻訳は機械翻訳により提供されています。提供された翻訳内容と英語版の間で齟齬、不一致または矛盾がある場合、英語版が優先します。

Sentinel-2 ラスターのデータ収集にアクセスし、土地セグメンテーションを実行する地球観測ジョブを作成する

この Python ベースのチュートリアルでは、 SDK for Python (Boto3) と Amazon SageMaker Studio Classic ノートブックを使用します。このデモを正常に完了するには、 SageMaker 地理空間を使用するために必要な AWS Identity and Access Management （IAM) アクセス許可があることを確認してください。Studio Classic. SageMaker 地理空間には、Studio Classic にアクセスできるユーザー、グループ、またはロールが必要です。また、信頼ポリシーsagemaker-geospatial.amazonaws.comで SageMaker 地理空間サービスプリンシパルを指定する SageMaker 実行ロールも必要です。

これらの要件の詳細については、SageMaker 「地理空間IAMロール」を参照してください。

このチュートリアルでは、 SageMaker 地理空間を使用して次のタスクAPIを完了する方法を示します。

list_raster_data_collections を使用して利用可能なデータ収集を検索する

SageMaker geospatial は、複数のラスターデータ収集をサポートしています。使用可能なデータ収集の詳細については、「」を参照してくださいデータ収集。

このデモでは、 から収集された衛星データを使用します。 Sentinel-2 クラウド最適化 GeoTIFF 衛星。これらの衛星は、5 日ごとに地球の地表をグローバルにカバーします。Earth の表面画像を収集することに加えて、Sentinel-2 衛星はさまざまなスペクトルバンドにわたってデータを収集します。

関心領域 (AOI) を検索するには、Sentinel-2 衛星データARNに関連付けられている が必要です。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

前のセクションでは、 ARNの を取得list_raster_data_collectionsするために を使用していました。Sentinel-2 データ収集。これを使用してARN、特定の関心領域 (AOI）、時間範囲、プロパティ、使用可能な UV バンドでデータ収集を検索できるようになりました。

を呼び出すにはsearch_raster_data_collectionAPI、Python dic は RasterDataCollectionQueryパラメータに対してオプションです。この例ではAreaOfInterest、、TimeRangeFilter、PropertyFilters、および を使用しますBandFilter。簡単にするために、 変数を使用して Python ディクショナリを指定search_rdc_queryして、検索クエリパラメータを保持できます。

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リクエストを実装します。これは、 get_paginator を使用したページ分割をサポートAPIしていませんAPI。完全なAPIレスポンスが収集されていることを確認するために、コードサンプルはwhileループを使用してNextToken存在することを確認します。次に、コードサンプルは .extend() を使用して、衛星画像URLsやその他のレスポンスメタデータを に追加しますitems_list。

の詳細についてはsearch_raster_data_collection、「Amazon SageMaker API リファレンスSearchRasterDataCollection」の「」を参照してください。

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レスポンスです。わかりやすくするために切り捨てられました。リクエストで"BandFilter": ["visual"]指定された のみがAssetsキーと値のペアで返されます。

{
    '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。これは、結果が正しい地理的リージョンからのものであることを検証するためです。

search_raster_data_collection を使用した の視覚化 matplotlib

地球観測ジョブ (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 geospatial には、ラスターデータ収集からの地理空間データを処理するために使用できる事前トレーニング済みのモデルが複数用意されています。利用可能な事前トレーニング済みモデルとカスタムオペレーションの詳細については、「」を参照してください操作タイプ。

水面面積の変化を計算するには、画像内のどのピクセルが水に対応するかを特定する必要があります。ランドカバーセグメンテーションは、 でサポートされているセマンティックセグメンテーションモデルstart_earth_observation_jobですAPI。セマンティックセグメンテーションモデルは、各画像内のすべてのピクセルにラベルを関連付けます。結果では、各ピクセルに、モデルのクラスマップに基づくラベルが割り当てられます。土地セグメンテーションモデルのクラスマップを次に示します。

{
    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_job API リクエストを行うときは、この変数を使用します。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 SageMaker API リファレンスGetEarthObservationJob」の「」を参照してください。

ARNs に関連付けられている を検索するには、 list_earth_observation_jobsAPIオペレーションEOJsを使用します。詳細については、「Amazon SageMaker API リファレンスListEarthObservationJobs」の「」を参照してください。

# 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 の表面積の変化を計算するには、まず を使用して の結果を Amazon S3 EOJにエクスポートしますexport_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) としてマークされたピクセル数をカウントします。2 つ目は、カウントに各ピクセルがカバーする領域を掛けます。バンドの空間解像度は異なる場合があります。ランドカバーセグメンテーションモデルでは、すべてのバンドが 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 月まで減少しました。