SageMaker 地理空間機能　植生画像を出力してみた

## 目次

1. [はじめに](#はじめに)
1. [NDVIとは](#ndviとは)
1. [NDVI画像の取得](#ndvi画像の取得)
1. [まとめ](#まとめ)

## はじめに

以前の記事で、SageMakerの地理空間機能を使用して衛星データを取得しました。今回は植生指数を出力しようと思います。植生指数はNDVIと呼ばれ(以降NDVIと記述)、植生の有無や活性度を把握するために使用されます。

SageMakergeospatialcapabilities の start_earth_observation_job を使用することで、NDVIを自動で算出し画像データとしてS3に保存することができます。

今回は、start_earth_observation_job を実行してNDVIの画像データを取得しました。

https://qiita.com/tosimitsu104/private/ad3754869788271e65fb

https://docs.aws.amazon.com/ja_jp/sagemaker/latest/dg/geospatial-eoj-ntb.html

https://boto3.amazonaws.com/v1/documentation/api/latest/reference/services/sagemaker-geospatial.html

## NDVIとは

NDVIとはNormalized Difference Vegetation Indexの略で、植物の活性度を表す指数です。光学衛星の可視域赤バンド(RED)と近赤外域バンド(NIR)を使用した以下の式から算出されます。

```math
NDVI = \frac{(NIR-RED)}{(NIR+RED)}
```

数値は-1から1の範囲で表され、1に近いほど植生が良い状態であることを示します。

https://sorabatake.jp/22019/

## NDVI画像の取得

まずはSageMaker geospatialのクライアントなど必要な諸々を準備します。

```python
import boto3
import sagemaker
import sagemaker_geospatial_map
import json

# SageMaker Geospatialのクライアントを作成
session = boto3.Session(region_name='us-west-2')
execution_role = sagemaker.get_execution_role()
geospatial_client = session.client(service_name="sagemaker-geospatial")
paginator = geospatial_client.get_paginator("list_raster_data_collections")
```

続いて、データセット一覧を出力、選択します。

```python
# 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'])
```
```python
# 利用可能なデータセットを出力
for collection in results[0]:
    print(f"Name: {collection['Name']}, ARN: {collection['Arn']}")
```

出力結果（Arnは一部伏せてます）

```bash
Name: Copernicus DEM GLO-30, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
Name: Copernicus DEM GLO-90, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
Name: Landsat Collection 2 Level-2 Science Products, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
Name: National Agriculture Imagery Program, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
Name: Sentinel 1 GRD, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
Name: Sentinel 2 L2A COGs, ARN: arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・
```
各種条件を指定します。
```python
# 衛星データのArn
SATELLITE_ARN = "arn:aws:sagemaker-geospatial:us-west-2:012345678910:raster-data-collection/public/・・・"
# 出力先S3のURI
OUTPUT_URI = "s3://sagemaker-studio-012345678910-・・・"
# 雲カバー率の条件
cloud_cover_property = {
    "EoCloudCover": {
        "LowerBound": 0,
        "UpperBound": 20
    },
}
# オレゴン州の領域
OREGON = [
     [-123.9150942493107, 46.183629329503916],
     [-124.1897524462052, 42.085988872609015],
     [-117.14751627783, 42.17561276955665],
     [-117.0376529990722, 45.801994745667486],
     [-123.9150942493107, 46.183629329503916]
]
# 観測期間
START_TIME = "2024-06-01T00:00:00Z"
END_TIME = "2024-08-31T23:59:59Z"
```

設定した条件から start_earth_observation_job で使用する設定項目を作成します。

job_config の PredefinedIndices には「NDVI」を指定します。

```python
input_config={
    'RasterDataCollectionQuery': {
        "AreaOfInterest": {
            "AreaOfInterestGeometry": {
                "PolygonGeometry": {
                    "Coordinates": [OREGON]
                }
            }
        },
        'RasterDataCollectionArn': SATELLITE_ARN,
        "TimeRangeFilter": {
            "StartTime": START_TIME,
            "EndTime": END_TIME
        },
        "PropertyFilters": {
            "Properties": [
                {
                    "Property": cloud_cover_property
                }
            ],
            "LogicalOperator": "AND"
        }
    }
}

job_config={
    'BandMathConfig': {
        'PredefinedIndices': ['NDVI']  # NDVIを指定
    }
}

name='ndvi-job-with-bandmath'
```

start_earth_observation_job を開始します。

```python
response = geospatial_client.start_earth_observation_job(
    ExecutionRoleArn=execution_role,
    InputConfig=input_config,
    JobConfig=job_config,
    Name=name
)

print("ジョブが開始されました:", response['Arn'])
```

```bash
# 出力結果
ジョブが開始されました: arn:aws:sagemaker-geospatial:us-west-2:012345678910:earth-observation-job/・・・
```

また、get_earth_observation_job を実行することでステータスを確認することができます。

```python
# ジョブステータスを確認する
import time

job_arn = response['Arn']

# ジョブのステータスを確認する関数
def check_job_status(client, job_arn):
    response = client.get_earth_observation_job(Arn=job_arn)
    return response['Status']

# ステータスの監視
status = check_job_status(geospatial_client, job_arn)
while status in ['IN_PROGRESS', 'STARTING']:
    print(f"ジョブの進行中: {status}")
    time.sleep(30)  # 30秒ごとにステータスを確認
    status = check_job_status(geospatial_client, job_arn)

print(f"ジョブが完了しました: {status}")
```

```bash
# 出力結果
ジョブの進行中: IN_PROGRESS
ジョブの進行中: IN_PROGRESS
ジョブの進行中: IN_PROGRESS
ジョブの進行中: IN_PROGRESS
...
ジョブの進行中: IN_PROGRESS
ジョブの進行中: IN_PROGRESS
ジョブの進行中: IN_PROGRESS
ジョブが完了しました: COMPLETED
```

start_earth_observation_job の実行結果をS3へ出力するため export_earth_observation_job を実行します。

```python
export_response = geospatial_client.export_earth_observation_job(
    Arn=job_arn,
    ExecutionRoleArn=execution_role_arn,
    OutputConfig={
        'S3Data': {
            'S3Uri': OUTPUT_URI
        }
    }
)

print("出力ジョブが開始されました")
```

また、export_earth_observation_job のステータスも確認できます。

```python
# export_earth_observation_jobのArn
export_job_arn = export_response['Arn']
# ステータスの監視
status = check_job_status(geospatial_client, export_job_arn)
while status in ['IN_PROGRESS', 'STARTING']:
    print(f"ジョブの進行中: {status}")
    time.sleep(30)  # 30秒ごとにステータスを確認
    status = check_job_status(geospatial_client, export_job_arn)

print(f"ジョブが完了しました: {status}")

# 出力結果のS3の場所を確認
output_uri = export_response['OutputConfig']['S3Data']['S3Uri']
print(f"出力データは次の場所に保存されました: {output_uri}")
print(f"作成時刻：{export_response['CreationTime']}")
```
```bash
# 出力結果
ジョブが完了しました: COMPLETED
出力データは次の場所に保存されました: s3://sagemaker-studio-012345678910-・・・
作成時刻：2024-10-22 03:18:42.167552+00:00
```

start_earth_observation_job により作成されたNDVI画像がS3に保存されていることを確認しましょう。

## まとめ

今回は、SageMakergeospatialcapabilities を使用してNDVI画像を作成、出力しました。可視域赤バンドと近赤外域バンドを取得して算出するコードを実装する手間が省け、楽にNDVIを取得できます。

SageMakerの自動化実行を行えば農業プロダクトに利用することもできそうです。

衛星データの活用の際は、SageMakerの使用も選択肢の一つとして誰かの参考になれば幸いです。