Accédez à la collection de données matricielles Sentinel-2 et créez une tâche d'observation de la Terre pour effectuer la segmentation des terres

Ce didacticiel basé sur Python utilise le SDK pour Python (Boto3) et un bloc-notes Amazon Studio Classic. SageMaker Pour mener à bien cette démonstration, assurez-vous que vous disposez des autorisations AWS Identity and Access Management (IAM) requises pour utiliser SageMaker Geospatial et Studio Classic. SageMaker geospatial nécessite que vous disposiez d'un utilisateur, d'un groupe ou d'un rôle pouvant accéder à Studio Classic. Vous devez également avoir un rôle SageMaker d'exécution qui spécifie le principal du service SageMaker géospatial, sagemaker-geospatial.amazonaws.com dans sa politique de confiance.

Pour en savoir plus sur ces exigences, consultez la section Rôles IAM SageMaker géospatiaux.

Ce didacticiel explique comment utiliser l'API SageMaker géospatiale pour effectuer les tâches suivantes :

• Trouvez les collections de données raster disponibles aveclist_raster_data_collections.

• Recherchez une collection de données raster spécifiée à l'aide desearch_raster_data_collection.

• Créez une tâche d'observation de la Terre (EOJ) en utilisantstart_earth_observation_job.

Utilisation list_raster_data_collections pour trouver les collections de données disponibles

SageMaker geospatial prend en charge plusieurs collections de données matricielles. Pour en savoir plus sur les collectes de données disponibles, voirCollections de données.

Cette démo utilise des données satellites collectées à partir de satellites GeoTIFF Sentinel-2 optimisés pour le cloud. Ces satellites fournissent une couverture mondiale de la surface terrestre de la Terre tous les cinq jours. En plus de collecter des images de surface de la Terre, les satellites Sentinel-2 collectent également des données sur diverses bandes spectrales.

Pour rechercher une zone d'intérêt (AOI), vous avez besoin de l'ARN associé aux données du satellite Sentinel-2. Pour trouver les collections de données disponibles et leurs ARN associés dans votre répertoire Région AWS, utilisez l'opération list_raster_data_collections API.

Comme la réponse peut être paginée, vous devez utiliser l'get_paginatoropération pour renvoyer toutes les données pertinentes :

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)

Il s'agit d'un exemple de réponse JSON provenant de l'opération d'list_raster_data_collectionsAPI. Il est tronqué pour inclure uniquement la collection de données (Sentinel-2) utilisée dans cet exemple de code. Pour plus de détails sur une collecte de données raster spécifique, utilisez 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"
}

Après avoir exécuté l'exemple de code précédent, vous obtenez l'ARN de la collection de données raster Sentinel-2,. arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8 Dans la section suivante, vous pouvez interroger la collecte de données Sentinel-2 à l'aide de l'search_raster_data_collectionAPI.

Recherche dans la collection de données Sentinel-2 matricielles à l'aide de search_raster_data_collection

Dans la section précédente, vous avez utilisé list_raster_data_collections pour obtenir l'ARN pour la collecte de Sentinel-2 données. Vous pouvez désormais utiliser cet ARN pour rechercher la collecte de données sur une zone d'intérêt (AOI) donnée, une plage de temps, des propriétés et les bandes UV disponibles.

Pour appeler l'search_raster_data_collectionAPI, vous devez transmettre Python un dictionnaire au RasterDataCollectionQuery paramètre. Cet exemple utilise AreaOfInterestTimeRangeFilter,PropertyFilters, etBandFilter. Pour plus de facilité, vous pouvez spécifier le dictionnaire Python à l'aide de la variable search_rdc_query pour contenir les paramètres de la requête de recherche :

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"
    ]
}

Dans cet exemple, vous recherchez une annonce AreaOfInterest qui inclut Lake Mead dans l'Utah. En outre, Sentinel-2 prend en charge plusieurs types de bandes d'images. Pour mesurer l'évolution de la surface de l'eau, vous n'avez besoin que du visual bracelet.

Après avoir créé les paramètres de requête, vous pouvez utiliser l'search_raster_data_collectionAPI pour effectuer la demande.

L'exemple de code suivant implémente une demande d'search_raster_data_collectionAPI. Cette API ne prend pas en charge la pagination à l'aide de l'get_paginatorAPI. Pour s'assurer que la réponse complète de l'API a été collectée, l'exemple de code utilise une while boucle pour vérifier qu'NextTokenelle existe. L'exemple de code est ensuite utilisé .extend() pour ajouter les URL des images satellites et d'autres métadonnées de réponse au. items_list

Pour en savoir plus à ce sujetsearch_raster_data_collection, consultez SearchRasterDataCollectionle Amazon SageMaker API Reference.

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))

Ce qui suit est une réponse JSON à votre requête. Il a été tronqué pour des raisons de clarté. Seul le paramètre "BandFilter": ["visual"] spécifié dans la demande est renvoyé dans la paire Assets clé-valeur :

{
    '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'
    }
}

Maintenant que vous avez obtenu les résultats de votre requête, dans la section suivante, vous pouvez visualiser les résultats en utilisantmatplotlib. Cela permet de vérifier que les résultats proviennent de la bonne région géographique.

Visualisation de votre utilisation search_raster_data_collectionmatplotlib

Avant de commencer le travail d'observation de la Terre (EOJ), vous pouvez visualiser le résultat de notre requête avec matplotlib. L'exemple de code suivant prend le premier élément de la items_list variable créée dans l'exemple de code précédent et imprime une image en utilisantmatplotlib. items_list[0]["Assets"]["visual"]["Href"]

# 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()

Après avoir vérifié que les résultats se situent dans la bonne région géographique, vous pouvez démarrer le Earth Observation Job (EOJ) à l'étape suivante. Vous utilisez l'EOJ pour identifier les plans d'eau à partir des images satellites à l'aide d'un processus appelé segmentation des terres.

Démarrage d'une tâche d'observation de la Terre (EOJ) qui effectue une segmentation du sol sur une série d'images satellites

SageMaker geospatial fournit plusieurs modèles préentraînés que vous pouvez utiliser pour traiter les données géospatiales issues de collections de données matricielles. Pour en savoir plus sur les modèles préentraînés et les opérations personnalisées disponibles, consultezTypes d'opérations.

Pour calculer la variation de la surface de l'eau, vous devez identifier les pixels des images qui correspondent à l'eau. La segmentation de la couverture terrestre est un modèle de segmentation sémantique soutenu par l'start_earth_observation_jobAPI. Les modèles de segmentation sémantique associent une étiquette à chaque pixel de chaque image. Dans les résultats, chaque pixel se voit attribuer une étiquette basée sur la carte de classes du modèle. Voici la carte des classes pour le modèle de segmentation des terres :

{
    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"
}

Pour démarrer une tâche d'observation de la Terre, utilisez l'start_earth_observation_jobAPI. Lorsque vous soumettez votre demande, vous devez spécifier les éléments suivants :

• InputConfig(dict) — Utilisé pour spécifier les coordonnées de la zone que vous souhaitez rechercher, ainsi que les autres métadonnées associées à votre recherche.

• JobConfig(dict) — Utilisé pour spécifier le type d'opération EOJ que vous avez effectuée sur les données. Cet exemple utilise LandCoverSegmentationConfig.

• ExecutionRoleArn(string) — L'ARN du rôle SageMaker d'exécution disposant des autorisations nécessaires pour exécuter la tâche.

• Name(chaîne) : nom de la tâche d'observation de la Terre.

InputConfigC'est un Python dictionnaire. Utilisez la variable suivante eoj_input_config pour contenir les paramètres de la requête de recherche. Utilisez cette variable lorsque vous faites la demande d'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",
        },
    }
}

JobConfigIl s'agit d'un Python dictionnaire utilisé pour spécifier l'opération EOJ que vous souhaitez effectuer sur vos données :

eoj_config = {"LandCoverSegmentationConfig": {}}

Les éléments du dictionnaire étant désormais spécifiés, vous pouvez envoyer votre demande d'start_earth_observation_jobAPI à l'aide de l'exemple de code suivant :

# 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)

Le démarrage d'une tâche d'observation de la Terre renvoie un ARN ainsi que d'autres métadonnées.

Pour obtenir une liste de toutes les tâches d'observation de la Terre en cours et en cours, utilisez l'list_earth_observation_jobsAPI. Pour surveiller l'état d'une seule tâche d'observation de la Terre, utilisez l'get_earth_observation_jobAPI. Pour effectuer cette demande, utilisez l'ARN créé après avoir soumis votre demande EOJ. Pour en savoir plus, consultez GetEarthObservationJoble Amazon SageMaker API Reference.

Pour trouver les ARN associés à vos EOJ, utilisez l'opération list_earth_observation_jobs API. Pour en savoir plus, consultez ListEarthObservationJobsle Amazon SageMaker API Reference.

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

Voici un exemple de réponse 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': {}
}

Une fois que le statut de votre tâche EOJ est passé àCOMPLETED, passez à la section suivante pour calculer la variation de la Mead's surface du lac.

Calcul de l'évolution de la Mead surface du lac

Pour calculer l'évolution de la superficie du lac Mead, exportez d'abord les résultats de l'EOJ vers Amazon S3 en utilisant : 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,
)

Pour connaître le statut de votre tâche d'exportation, utilisez get_earth_observation_job :

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

Pour calculer les variations du niveau d'eau du lac Mead, téléchargez les masques de couverture terrestre sur l'instance de SageMaker bloc-notes locale et téléchargez les images sources de notre requête précédente. Dans la carte des classes du modèle de segmentation des terres, l'indice de classe de l'eau est 6.

Pour extraire le masque à eau d'une Sentinel-2 image, procédez comme suit. Commencez par compter le nombre de pixels marqués comme de l'eau (indice de classe 6) dans l'image. Ensuite, multipliez le nombre par la zone couverte par chaque pixel. Les bandes peuvent avoir une résolution spatiale différente. Pour le modèle de segmentation de la couverture terrestre, toutes les bandes sont sous-échantillonnées à une résolution spatiale égale à 60 mètres.

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()

À l'aide dematplotlib, vous pouvez visualiser les résultats sous forme de graphique. Le graphique montre que la superficie du lac Mead a diminué entre janvier 2021 et juillet 2022.