{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "# Digital Farming with Amazon SageMaker Geospatial Capabilities - Part I\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "This notebook's CI test result for us-west-2 is as follows. CI test results in other regions can be found at the end of the notebook. \n", "\n", "\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "\n", "In this notebook, we will explore some of most common tasks for processing geospatial data in the Digital Farming domain, by working with Amazon SageMaker geospatial capabilities.\n", "\n", "----" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Environment Set-Up\n", "\n", "We will start by importing a few libraries required, the necessary dependencies are already pre-installed in the Geospatial kernel." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "import boto3\n", "import sagemaker\n", "import sagemaker_geospatial_map\n", "\n", "import json\n", "from datetime import datetime\n", "import rasterio\n", "from rasterio.plot import show\n", "from matplotlib import pyplot as plt\n", "import pandas as pd\n", "import time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now define a few variables for which we need to create sessions in the SageMaker and Boto3 SDKs.\n", "\n", "We will also create the client for SageMaker geospatial capabilities by creating a Botocore session..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "sagemaker_session = sagemaker.Session()\n", "bucket = sagemaker_session.default_bucket() ### Replace with your own bucket if needed\n", "role = sagemaker.get_execution_role(sagemaker_session)\n", "prefix = \"sm-geospatial-e2e\" ### Replace with the S3 prefix desired\n", "region = boto3.Session().region_name\n", "print(f\"S3 bucket: {bucket}\")\n", "print(f\"Role: {role}\")\n", "print(f\"Region: {region}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note, at the time of writting this notebook sagemaker-geospatial is only available in the region \"us-west-2\".**\n", "\n", "Make sure you have the proper policy allowed and a trust relationship added to your role for \"sagemaker-geospatial\", as specified in the [Get Started with Amazon SageMaker Geospatial Capabilties](https://docs.aws.amazon.com/sagemaker/latest/dg/geospatial-getting-started.html) documentation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "gsClient = boto3.client(\"sagemaker-geospatial\", region_name=region)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "\n", "## Common geospatial processing tasks for Digital Farming\n", "\n", "Now that we have the client setup, we are ready to work with SageMaker Geospatial.\n", "\n", "We can in example check the raster collections available..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "list_raster_data_collections_resp = gsClient.list_raster_data_collections()\n", "list_raster_data_collections_resp[\"RasterDataCollectionSummaries\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will save the ARN of our collection of interest. In our example we will work with satellite imagery data from the [Sentinel-2-L2A](https://registry.opendata.aws/sentinel-2-l2a-cogs/) collection..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "data_collection_arn = \"arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8\"\n", "### Replace with the ARN of the collection of your choice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we will define the input configuration with the polygon of coordinates for our area of interest and the time range we are interested on." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "### Replace with the coordinates for the polygon of your area of interest...\n", "coordinates = [\n", " [9.742977, 53.615875],\n", " [9.742977, 53.597119],\n", " [9.773620, 53.597119],\n", " [9.773620, 53.615875],\n", " [9.742977, 53.615875],\n", "]\n", "### Replace with the time-range of interest...\n", "time_start = \"2022-03-01T12:00:00Z\"\n", "time_end = \"2022-03-31T12:00:00Z\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could graphically visualize this area by using the SageMaker Studio Geospatial extension, but in our case we will use an open-source library for programatically embedding a map in the notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "def get_lat_from_coordinates(coord):\n", " lat = []\n", " for i, j in enumerate(coord):\n", " lat.append(coord[i][1])\n", " return lat\n", "\n", "\n", "def get_lon_from_coordinates(coord):\n", " lon = []\n", " for i, j in enumerate(coord):\n", " lon.append(coord[i][0])\n", " return lon" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Instantiate a new map and display it\n", "Map = sagemaker_geospatial_map.create_map()\n", "Map.set_sagemaker_geospatial_client(gsClient)\n", "\n", "Map.render()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "dataset = Map.add_dataset(\n", " {\n", " \"data\": pd.DataFrame(\n", " {\n", " \"Latitude\": get_lat_from_coordinates(coordinates),\n", " \"Longitude\": get_lon_from_coordinates(coordinates),\n", " }\n", " )\n", " }, \n", " auto_create_layers=True\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Typically, we are interested on working with images that are not covered by much clouds over our area of interest. For exploring this in our notebook, we will define some additional parameters like e.g. the ranges for cloud cover we want to consider (less than 2% in our example)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "eoj_input_config = {\n", " \"RasterDataCollectionQuery\": {\n", " \"AreaOfInterest\": {\n", " \"AreaOfInterestGeometry\": {\"PolygonGeometry\": {\"Coordinates\": [coordinates]}}\n", " },\n", " \"TimeRangeFilter\": {\"StartTime\": time_start, \"EndTime\": time_end},\n", " \"PropertyFilters\": {\n", " \"Properties\": [{\"Property\": {\"EoCloudCover\": {\"LowerBound\": 0, \"UpperBound\": 2}}}],\n", " },\n", " }\n", "}\n", "\n", "data_collection_arn = \"arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now search for satellite images of the area of interest and visualize, as an example, the True Color Image (TCI) tiles provided by SageMaker, or any other band of interest.\n", "\n", "Note these images are directly provided by the built-in data sources in SageMaker geospatial capabilities, like e.g. AWS Data Exchange in our case, and the query will return the locations of these image tiles in Amazon S3." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "response = gsClient.search_raster_data_collection(**eoj_input_config, Arn=data_collection_arn)\n", "response" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "for i, j in enumerate(response[\"Items\"]):\n", " print(response[\"Items\"][i][\"Assets\"][\"visual\"][\"Href\"])\n", " with rasterio.open(response[\"Items\"][i][\"Assets\"][\"visual\"][\"Href\"]) as src:\n", " arr = src.read(out_shape=(src.height // 40, src.width // 40))\n", " show(arr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cloud Gap Filling - Earth Observation Job\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Currently, SageMaker Geospatial supports different types of built-in processing Earth Observation Jobs (EOJs), such as:\n", "\n", "* Cloud removal\n", "* Temporal or zonal statistics\n", "* Stacking of bands\n", "* Geo-Mosaic combinations\n", "* Re-sampling\n", "* Segmentation with built-in models\n", "\n", "etc.\n", "\n", "Let us explore some of these techniques with an example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some cases we might be interested on removing the cloudy pixels of the images, this could be e.g. if we don't have any clear observation available in the data for our area and time-range of interest. We can do this through a **Cloud Removal** EOJ, which would replace or fill these cloudy pixels using any of the techniques available.\n", "\n", "For this we will setup the cloud removal parameters. This can e.g. be set to a fixed value to run your own gap filling algorithm/logic, or alternatively use a built-in interpolation type supported in the service." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "eoj_config = {\n", " \"CloudRemovalConfig\": {\"AlgorithmName\": \"INTERPOLATION\", \"InterpolationValue\": \"-9999\"}\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we are ready to start the EOJ with the given config. We will define a function for this purpose that will be handy for the other tasks in this notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "eoj_input_config[\"RasterDataCollectionQuery\"][\"RasterDataCollectionArn\"] = data_collection_arn\n", "\n", "\n", "def start_earth_observation_job(eoj_name, role, eoj_input_config, eoj_config):\n", " # Start EOJ...\n", " response = gsClient.start_earth_observation_job(\n", " Name=eoj_name,\n", " ExecutionRoleArn=role,\n", " InputConfig=eoj_input_config,\n", " JobConfig=eoj_config,\n", " )\n", " eoj_arn = response[\"Arn\"]\n", " print(f\"{datetime.now()} - Started EOJ: {eoj_arn}\")\n", "\n", " # Wait for EOJ to complete... check status every minute\n", " gs_get_eoj_resp = {\"Status\": \"IN_PROGRESS\"}\n", " while gs_get_eoj_resp[\"Status\"] == \"IN_PROGRESS\":\n", " time.sleep(60)\n", " gs_get_eoj_resp = gsClient.get_earth_observation_job(Arn=eoj_arn)\n", " print(f'{datetime.now()} - Current EOJ status: {gs_get_eoj_resp[\"Status\"]}')\n", " return eoj_arn" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "cr_eoj_arn = start_earth_observation_job(\n", " f'cloudremovaljob-{datetime.now().strftime(\"%Y-%m-%d-%H-%M\")}',\n", " role,\n", " eoj_input_config,\n", " eoj_config,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note the EOJ processing takes some minutes.** We can check the status programatically by getting the EOJ with the SageMaker Geospatial client, or graphically by using the Geospatial extension for SageMaker Studio." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once completed, note the resulting imagery is stored in a service bucket and it can be used for chaining this to another EOJ. Also, for visualizing these images we can, if desired, copy the results of the EOJ to our own S3 location in our account by exporting it.\n", "\n", "Check the \"Exporting the Results\" section of this notebook for details." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Geo Mosaic - Earth Observation Job\n", "\n", "A common process when working with tiles of Satellite images is combining these geographically, for having a consolidated view of the whole area of interest. We can perform this through the use of a Geo Mosaic job, applicable to the output of most of the examples provided with other processing EOJs in this notebook.\n", "\n", "For this purpose we can set the configuration as follows." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "eoj_config = {\"GeoMosaicConfig\": {\"AlgorithmName\": \"NEAR\"}}\n", "\n", "gm_eoj_arn = start_earth_observation_job(\n", " f'geomosaicjob-{datetime.now().strftime(\"%Y-%m-%d-%H-%M\")}', role, eoj_input_config, eoj_config\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note the EOJ processing takes some minutes.** We can check the status programatically by getting the EOJ with the SageMaker Geospatial client, or graphically by using the Geospatial extension for SageMaker Studio." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Band-Math - Earth Observation Job\n", "\n", "Follwing our example, we will now perform a band-math over predefined vegetation and soil indices supported by SageMaker geospatial capabilities. \n", "Note SageMaker geospatial capabilities can perform the band-math over one or many of these indices, and supports the use of equations.\n", "\n", "For our case let's obtain the **moisture index** by using the equation below." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "eoj_config = {\n", " \"BandMathConfig\": {\n", " \"CustomIndices\": {\n", " \"Operations\": [{\"Name\": \"moisture\", \"Equation\": \"(nir08 - swir16) / (nir08 + swir16)\"}]\n", " }\n", " }\n", "}\n", "\n", "bm_eoj_arn = start_earth_observation_job(\n", " f'bandmathjob-{datetime.now().strftime(\"%Y-%m-%d-%H-%M\")}', role, eoj_input_config, eoj_config\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note the EOJ processing takes some minutes.** We can check the status programatically by getting the EOJ with the SageMaker Geospatial client, or graphically by using the Geospatial extension for SageMaker Studio." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----\n", "\n", "## Exporting the Results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As mentioned before, the results of our EOJs are stored in the service and are available for chaining as input for another EOJ, but can also export these to Amazon S3 for visualizing the imagery directly.\n", "\n", "We will define a function for exporting the results of our EOJs." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def export_earth_observation_job(eoj_arn, role, bucket, prefix, task_suffix):\n", " # Export EOJ results to S3...\n", " response = gsClient.export_earth_observation_job(\n", " Arn=eoj_arn,\n", " ExecutionRoleArn=role,\n", " OutputConfig={\n", " \"S3Data\": {\"S3Uri\": f\"s3://{bucket}/{prefix}/{task_suffix}/\", \"KmsKeyId\": \"\"}\n", " },\n", " )\n", " export_arn = response[\"Arn\"]\n", " print(f\"{datetime.now()} - Exporting with ARN: {export_arn}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's go through the EOJs created before for checking it's status and exporting accordingly. Keep in mind each EOJ takes some minutes to complete, so we will add a check on the status every 30 seconds..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Check status of EOJs...\n", "EOJs = [cr_eoj_arn, gm_eoj_arn, bm_eoj_arn]\n", "eoj_suffixes = [\"cloud_removal\", \"geo_mosaic\", \"band_math\"]\n", "\n", "eoj_status = [\"\"] * len(EOJs)\n", "while not all(i == \"Exported\" for i in eoj_status):\n", " # Wait for EOJs to complete and export... check status every 30 seconds\n", " for j, eoj in enumerate(EOJs):\n", " gs_get_eoj_resp = gsClient.get_earth_observation_job(Arn=eoj)\n", " if gs_get_eoj_resp[\"Status\"] == \"COMPLETED\":\n", " # EOJ completed, exporting...\n", " if not \"ExportStatus\" in gs_get_eoj_resp:\n", " export_earth_observation_job(eoj, role, bucket, prefix, eoj_suffixes[j])\n", " elif gs_get_eoj_resp[\"ExportStatus\"] == \"IN_PROGRESS\":\n", " eoj_status[j] = \"Exporting\"\n", " elif gs_get_eoj_resp[\"ExportStatus\"] == \"SUCCEEDED\":\n", " eoj_status[j] = \"Exported\"\n", " else:\n", " raise Exception(\"Error exporting\")\n", " elif gs_get_eoj_resp[\"Status\"] == \"IN_PROGRESS\":\n", " # EOJ still in progress, keep waiting...\n", " eoj_status[j] = \"In progress\"\n", " else:\n", " raise Exception(\"Error with the EOJ\")\n", " print(f\"{datetime.now()} - EOJ: {eoj} Status: {eoj_status[j]}\")\n", " if all(i == \"Exported\" for i in eoj_status):\n", " break\n", " time.sleep(30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have all our EOJs exported, let's visualize a few of the images obtained in S3.\n", "\n", "For this we will use the open library \"rasterio\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "s3 = boto3.resource(\"s3\")\n", "my_bucket = s3.Bucket(bucket)\n", "\n", "\n", "def visualize_cogt(task, eoj_arn, band, number):\n", " gs_get_eoj_resp = gsClient.get_earth_observation_job(Arn=eoj_arn)\n", " if gs_get_eoj_resp[\"ExportStatus\"] == \"SUCCEEDED\":\n", " i = 0\n", " for index, image in enumerate(\n", " my_bucket.objects.filter(\n", " Prefix=f'{prefix}/{task}/{eoj_arn.split(\"/\",1)[1]}/output/consolidated/'\n", " )\n", " ):\n", " if f\"{band}.tif\" in image.key:\n", " i = i + 1\n", " tif = f\"s3://{bucket}/{image.key}\"\n", " with rasterio.open(tif) as src:\n", " arr = src.read(out_shape=(src.height // 20, src.width // 20))\n", " if band != \"visual\":\n", " # Sentinel-2 images are stored as uint16 for optimizing storage\n", " # but these need to be reslaced (by dividing each pixel value by 10000)\n", " # to get the true reflectance values. This is a common \u201ccompression\u201d\n", " # technique when storing satellite images...\n", " arr = arr / 10000\n", " # As a result of the transformation, there might be some pixel values\n", " # over 1 in the RGB, so we need to replace those by 1...\n", " arr[arr > 1] = 1\n", " show(arr)\n", " print(tif)\n", " if i == number:\n", " break\n", " else:\n", " print(\n", " f'Export of job with ARN:\\n{eoj_arn}\\nis in ExportStatus: {gs_get_eoj_resp[\"ExportStatus\"]}'\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For our Cloud Removal, we can check in example some of the images for the \"blue\" band" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "visualize_cogt(\"cloud_removal\", cr_eoj_arn, \"blue\", 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the Geo Mosaic, let's visualize the consolidated image shown as \"visual\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "visualize_cogt(\"geo_mosaic\", gm_eoj_arn, \"visual\", 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the Band Math, let's visualize a few of the output images for the Moisture Index." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "visualize_cogt(\"band_math\", bm_eoj_arn, \"moisture\", 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, take into account the legend for the resulting moisture index would be similar to the below.\n", "\n", "