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

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

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

これらの要件の詳細については、「SageMaker geospatial IAM roles」を参照してください。

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

list_raster_data_collections を使用して利用可能なデータコレクションを検索する

SageMaker 地理空間は、複数のラスターデータコレクションをサポートしています。利用可能なデータコレクションの詳細については、「データコレクション」を参照してください。

このデモでは、Sentinel-2Cloud-Optimized GeoTIFF 衛星から収集された衛星データを使用します。これらの衛星は、5 日ごとに世界中をカバーする地球の地表を提供します。地球の表面画像を収集することに加えて、Sentinel-2 衛星はさまざまなスペクトルバンドのデータも収集します。

対象エリア (AOI) を検索するには、Sentinel-2 衛星データに関連付けられている ARN が必要です。で使用可能なデータ収集とそれに関連する ARNs を検索するには AWS リージョン、 list_raster_data_collections API オペレーションを使用します。

レスポンスはページ分割される場合があるため、 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)

これは API list_raster_data_collections オペレーションからの 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 が得られます。次のセクション では、search_raster_data_collection API を使用して Sentinel-2 データコレクションのクエリを実行できます。

search_raster_data_collection を使用して Sentinel-2 ラスターデータコレクションを検索する

前のセクションでは、 list_raster_data_collections を使用して Sentinel-2 データコレクションの ARN を取得しました。そのためその ARN を使用して、特定の対象エリア (AOI)、時間範囲、プロパティ、使用可能な UV バンドのデータコレクションを検索できるようになりました。

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

の詳細についてはsearch_raster_data_collection、Amazon SageMaker AI 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 を使用して結果を視覚化できます。これは、結果が正しい地域からもたらされたものであることを検証するためです。

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_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 リクエストを行うときは、この変数を使用します。

# 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 は、データに対して実行する EOJ オペレーションを指定するために使用する Python ディクショナリです。

eoj_config = {"LandCoverSegmentationConfig": {}}

ディクショナリ要素を指定したら、次のコードサンプルを使用して start_earth_observation_job API リクエストを送信できます。

# 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 GetEarthObservationJob」を参照してください。

EOJ に関連付けられている ARN を検索するには、list_earth_observation_jobs API オペレーションを使用します。詳細については、Amazon SageMaker AI 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 の表面積の変化を計算するには、まず export_earth_observation_job を使用して EOJ の結果を Amazon S3 にエクスポートします。

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 月にかけて減少したことを示しています。