使用 Amazon SageMaker 在 20 分钟内绘制地球植被地图 机器学习博客

使用 Amazon SageMaker 在 20 分钟内映射地球植被

重点摘要

在不断变化的世界中，监测地球植被健康变得尤为重要。本文介绍了如何利用 Amazon SageMaker 的地理空间能力，快速完成全球植被的映射，展示了其在环境监测和可持续发展方面的潜力。

当今世界，监测地球植被健康比以往任何时候都更加重要。植被在维持生态平衡、提供食物和作为碳汇方面发挥着至关重要的作用。传统的植被健康监测方法如实地调查和卫星数据人工分析，不仅耗时且资源需求大，同时也需要特定领域的专业知识。这些传统方法操作繁琐，常常导致数据收集和分析的延迟，使得对环境变化的快速跟踪和响应变得困难。此外，这些方法的高成本限制了它们的可及性和频率，妨碍了全球植被监测的全面性和持续性。鉴于这些挑战，我们开发了一种创新解决方案，以简化和提高全球植被监测的效率。

通过转变传统的劳动密集型植被健康监测方法，Amazon SageMaker 地理空间能力 提供了一种简便且经济有效的解决方案。Amazon SageMaker 支持地理空间机器学习ML能力，使数据科学家和 ML 工程师能够利用地理空间数据构建、训练和部署 ML 模型。这些地理空间能力为环境监测开辟了一种全新的可能性。使用 SageMaker，用户可以访问大量地理空间数据集，快速处理和丰富这些数据，并加快开发进度。过去需要数天或数周才能完成的任务，如今只需短短几分钟即可完成。

在本篇文章中，我们展示了 SageMaker 地理空间能力的强大之处，通过在 20 分钟内绘制全球植被。这一例子不仅强调了 SageMaker 的效率，也展示了地理空间 ML 如何被用于环境监测的可持续性和保护目的。

确定感兴趣区域

我们首先演示如何在全球范围内应用 SageMaker 分析地理空间数据。我们根据开始使用 Amazon SageMaker 地理空间能力 中的步骤进行操作。我们首先定义一个定义兴趣区域的地理坐标边界框。该边界框用作过滤器，仅选择覆盖地球陆地的相关卫星图像。

pythonimport osimport jsonimport timeimport boto3import geopandasfrom shapelygeometry import Polygonimport leafmapfoliumap as leafmapimport sagemakerimport sagemakergeospatialmap

session = boto3Session()executionrole = sagemakergetexecutionrole()sgclient = sessionclient(servicename=sagemakergeospatial)cooridinates =[ [179034845 55973798] [179371094 55973798] [179371094 83780085] [179034845 83780085] [179034845 55973798]] polygon = Polygon(cooridinates)worldgdf = geopandasGeoDataFrame(index=[0] crs=epsg4326 geometry=[polygon])m = leafmapMap(center=[37 119] zoom=4)maddbasemap(EsriWorldImagery)maddgdf(worldgdf layername=AOI style={color red})m

数据获取

SageMaker 的地理空间能力提供了对多种公共地理空间数据集的访问，包含Sentinel2、Landsat 8、Copernicus DEM 和 NAIP。在我们的植被映射项目中，我们选择了 Sentinel2，因为其覆盖范围广且更新频繁。Sentinel2卫星每5天获取一次10米分辨率的地球陆面图像。我们在本例中选择2023年12月的第一周。为了确保捕获更多的清晰可见地表，我们筛选出云覆盖度低于10的图像。这样，我们的分析基于清晰且可靠的图像。

pythonsearchrdcargs = { Arn arnawssagemakergeospatialuswest2378778860802rasterdatacollection/public/nmqj48dcu3g7ayw8 # sentinel2 L2A RasterDataCollectionQuery { AreaOfInterest { AreaOfInterestGeometry { PolygonGeometry { Coordinates [ [ [179034845 55973798] [179371094 55973798] [179371094 83780085] [179034845 83780085] [179034845 55973798] ] ] } } TimeRangeFilter { StartTime 20231201T000000Z EndTime 20231207T235959Z } PropertyFilters { Properties [{Property {EoCloudCover {LowerBound 0 UpperBound 10}}}] LogicalOperator AND } }}

s2items = []s2tileids = []s2geometries = { id [] geometry []}while searchrdcargsget(NextToken True) searchresult = sgclientsearchrasterdatacollection(searchrdcargs) for item in searchresult[Items] s2id = item[Id] s2tileid = s2idsplit()[1] # 筛选出覆盖相同区域的图像 if s2tileid not in s2tileids s2tileidsappend(s2tileid) s2geometries[id]append(s2id) s2geometries[geometry]append(Polygon(item[Geometry][Coordinates][0])) del item[DateTime] s2itemsappend(item)

searchrdcargs[NextToken] = searchresultget(NextToken)

print(f{len(s2items)} unique Sentinel2 images found)

通过利用 SageMaker 地理空间的 searchrasterdatacollection 函数，我们在 2023 年 12 月的第一周找到了 8581 个独特的 Sentinel2 图像。为了验证我们选择的准确性，我们在地图上绘制了这些图像的覆盖区域，确认我们获取的图像是正确的。

pythons2gdf = geopandasGeoDataFrame(s2geometries)m = leafmapMap(center=[37 119] zoom=4)maddbasemap(OpenStreetMap)maddgdf(s2gdf layername=Sentinel2 Tiles style={color blue})m

SageMaker 地理空间处理作业

在使用 SageMaker 地理空间能力查询数据时，我们获得了有关目标图像的全面信息，包括数据轮廓、光谱波段属性和直接访问的超链接。有了这些超链接，我们可以避免传统的先下载图像再本地处理的内存和存储密集型方法这项工作因数据集的大小和规模而变得尤为庞大，数据集总规模超过4 TB。每张图像的大小都很大，具有多个通道，单独大小约为 500 MB。在单机上处理数 TB 的数据将是耗时的。尽管设置处理集群是另一种选择，但这也带来了从数据分发到基础设施管理的一系列复杂性。SageMaker 地理空间通过 Amazon SageMaker 处理 简化了这一过程。我们使用专门构建的地理空间容器与 SageMaker 处理作业进行简化管理体验，创建并运行集群。只需几行代码，即可通过 SageMaker 处理作业扩展地理空间工作负载。您只需指定定义工作负载的脚本、地理空间数据的位置以及地理空间容器。SageMaker 处理根据您的城市、国家或大洲规模的地理空间 ML 工作负载为您提供集群资源。

在我们的项目中，我们使用 25 个集群，每个集群包含 20 个实例，以扩展我们的地理空间工作负载。接下来，我们将 8581 张图像拆分为 25 个批次，以便于高效处理。每个批次大约包含 340 张图像。这些批次随后平均分配到集群中的机器上。所有批次清单都上传到 Amazon S3，便于处理作业准备，这样每个部分就可以迅速高效地进行处理。

pythondef s2itemtorelativemetadataurl(item) parts = item[Assets][visual][Href]split(/) tileprefix = parts[41] return {}/{}jsonformat(/join(tileprefix) item[Id])

numjobs = 25numinstancesperjob = 20 # 最大20

manifestlist = {}for idx in range(numjobs) manifest = [{prefix s3//sentinelcogs/sentinels2l2acogs/}] manifestlist[idx] = manifest

为 N 个处理作业拆分清单

for idx item in enumerate(s2items) jobidx = idx numjobs manifestlist[jobidx]append(s2itemtorelativemetadataurl(item))

上传清单到 S3

sagemakersession = sagemakerSession()s3bucketname = sagemakersessiondefaultbucket()s3prefix = processingjobdemos3client = boto3client(s3)s3 = boto3resource(s3)

manifestdir = manifestsosmakedirs(manifestdir existok=True)

for jobidx manifest in manifestlistitems() manifestfile = f{manifestdir}/manifest{jobidx}json s3manifestkey = s3prefix / manifestfile with open(manifestfile w) as f jsondump(manifest f)

s3clientuploadfile(manifestfile s3bucketname s3manifestkey)print(Uploaded {} to {}format(manifestfile s3manifestkey))

准备好输入数据后，我们现在可以进行核心分析，这将揭示植被健康的见解，通过归一化差异植被指数NDVI。NDVI 通过近红外NIR与红色反射率之间的差异计算得出，并由其总和进行归一化，数值范围从 1 到 1。更高的 NDVI 值表示浓密而健康的植被，零值则表示没有植被，而负值通常指向水体。此指标是评估植被健康和分布的重要工具。以下是 NDVI 的示例。

pythonwritefile scripts/computevipy

import osimport rioxarrayimport jsonimport gcimport warnings

warningsfilterwarnings(ignore)

if name == main print(Starting processing)

inputpath = /opt/ml/processing/inputoutputpath = /opt/ml/processing/outputinputfiles = []items = []for currentpath subdirs files in oswalk(inputpath)    for file in files        if fileendswith(json)            fullfilepath = ospathjoin(inputpath currentpath file)            inputfilesappend(fullfilepath)            with open(fullfilepath r) as f                itemsappend(jsonload(f))print(Received {} input filesformat(len(inputfiles)))for item in items    print(Computing NDVI for {}format(item[id]))    redbandurl = item[assets][red][href]    nirbandurl = item[assets][nir][href]    sclmaskurl = item[assets][scl][href]    red = rioxarrayopenrasterio(redbandurl masked=True)    nir = rioxarrayopenrasterio(nirbandurl masked=True)    scl = rioxarrayopenrasterio(sclmaskurl masked=True)    sclinterp = sclinterp(        x=red[x] y=red[y]    )  # 将 SCL 插值为与红和 NIR 波段相同的分辨率    # 使用 SCL 遮蔽多云像素 (https//sentinelscopernicuseu/web/sentinel/technicalguides/sentinel2msi/level2a/algorithmoverview)    # 类别 8：中等概率云    # 类别 9：高概率云    # 类别 10：薄卷云    redcloudmasked = redwhere((sclinterp != 8) amp (sclinterp != 9) amp (sclinterp != 10))    nircloudmasked = nirwhere((sclinterp != 8) amp (sclinterp != 9) amp (sclinterp != 10))    ndvi = (nircloudmasked  redcloudmasked) / (nircloudmasked  redcloudmasked)    # 将 ndvi 保存为 geotiff    s2tileid = redbandurlsplit(/)[2]    filename = f{s2tileid}ndvitif    outputfilepath = f{outputpath}/{filename}    ndviriotoraster(outputfilepath)    print(Written output {}format(outputfilepath))    # 保持低内存使用    del red    del nir    del scl    del sclinterp    del redcloudmasked    del nircloudmasked    del ndvi    gccollect()

现在我们已经定义了计算逻辑，准备开始地理空间 SageMaker 处理作业。这包括简单的三步过程：设置计算集群、定义计算细节以及组织输入和输出信息。

首先，设置集群时，我们决定所需的实例数量和类型，确保它们适合地理空间数据处理。计算环境本身通过选择一个包含所有常用地理空间数据处理包的地理空间映像来准备。

接下来，就输入而言，我们使用之前创建的列出所有图像超链接的清单。我们还指定了一个 S3 位置来保存我们的结果。

配置好这些要素后，我们可以一次启动多个处理作业，使它们并行工作以提高效率。

pythonfrom multiprocessing import Processimport sagemakerimport boto3 from botocoreconfig import Configfrom sagemaker import getexecutionrolefrom sagemakersklearnprocessing import ScriptProcessorfrom sagemakerprocessing import ProcessingInput ProcessingOutput

role = getexecutionrole()geospatialimageuri = 081189585635dkrecruswest2amazonawscom/sagemakergeospatialv10latest

使用 boto3 的重试行为以避免限制问题

smboto = boto3client(sagemaker config=Config(connecttimeout=5 readtimeout=60 retries={maxattempts 20}))sagemakersession = sagemakerSession(sagemakerclient=smboto)

def runjob(jobidx) s3manifest = fs3//{s3bucketname}/{s3prefix}/{manifestdir}/manifest{jobidx}json s3output = fs3//{s3bucketname}/{s3prefix}/output scriptprocessor = ScriptProcessor( command=[python3] imageuri=geospatialimageuri role=role instancecount=numinstancesperjob instancetype=mlm5xlarge basejobname=fcas2nvdi{jobidx} sagemakersession=sagemakersession )

scriptprocessorrun(    code=scripts/computevipy    inputs=[        ProcessingInput(            source=s3manifest            destination=/opt/ml/processing/input/            s3datatype=ManifestFile            s3datadistributiontype=ShardedByS3Key        )    ]    outputs=[        ProcessingOutput(            source=/opt/ml/processing/output/            destination=s3output            s3uploadmode=Continuous        )    ])timesleep(2)

processes = []for idx in range(numjobs) p = Process(target=runjob args=(idx)) processesappend(p) pstart()

for p in processes pjoin()

启动