sagemaker-geospatial/mount-shasta-glacier-melting-monitoring/mount_shasta_glacier_melt_monitoring.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "09323c51-6cb3-4a7d-9c7a-ef3b41e63248",
   "metadata": {},
   "source": [
    "# Monitoring Mountain Glacier Melting with SageMaker Geospatial Capabilities\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",
    "![This us-west-2 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/us-west-2/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09323c51-6cb3-4a7d-9c7a-ef3b41e63248",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "This notebook runs with Kernel Geospatial 1.0. Note that the following policies need to be attached to the execution role that you used to run this notebook:\n",
    "\n",
    "- AmazonSageMakerFullAccess\n",
    "- AmazonSageMakerGeospatialFullAccess\n",
    "\n",
    "You can see the policies attached to the role in the IAM console under the permissions tab. If required, add the roles using the 'Add Permissions' button. \n",
    "\n",
    "In addition to these policies, ensure that the execution role's trust policy allows the SageMaker-GeoSpatial service to assume the role. This can be done by adding the following trust policy using the 'Trust relationships' tab:\n",
    "```\n",
    "{\n",
    "    \"Version\": \"2012-10-17\",\n",
    "    \"Statement\": [\n",
    "        {\n",
    "            \"Effect\": \"Allow\",\n",
    "            \"Principal\": {\n",
    "                \"Service\": [\n",
    "                    \"sagemaker.amazonaws.com\",\n",
    "                    \"sagemaker-geospatial.amazonaws.com\"\n",
    "                ]\n",
    "            },\n",
    "            \"Action\": \"sts:AssumeRole\"\n",
    "        }\n",
    "    ]\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04a08fa7-a784-40e8-a36d-cfa84830b66b",
   "metadata": {},
   "source": [
    "**Contents**\n",
    "* [Setup SageMaker geospatial capabilities](#1)\n",
    "* [Query Sentinel-2 Data](#2)\n",
    "* [Start an Earth Observation Job (EOJ) to identify the land cover types in the area of Mount Shasta](#3)\n",
    "* [Visualize EOJ inputs and outputs in FourSquare Studio](#4)\n",
    "* [Extract the land cover segmentation results](#5)\n",
    "* [Measure snow coverage of Mount Shasta](#6)\n",
    "* [Query Landsat 8 data](#7)\n",
    "* [Analyze the relationship between snow coverage and surface temperature](#8)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62e59d15-3a3e-4eae-8100-50b646425a12",
   "metadata": {},
   "source": [
    "<a id='1'></a>\n",
    "\n",
    "## Setup SageMaker geospatial capabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb593738-f41e-4bbe-bab3-12e4e5beab26",
   "metadata": {},
   "outputs": [],
   "source": [
    "import boto3\n",
    "import sagemaker\n",
    "import sagemaker_geospatial_map\n",
    "\n",
    "session = boto3.Session()\n",
    "execution_role = sagemaker.get_execution_role()\n",
    "sg_client = session.client(service_name=\"sagemaker-geospatial\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "767b519b-4b00-4702-9e99-34d91ba0c720",
   "metadata": {},
   "source": [
    "<a id='2'></a>\n",
    "\n",
    "## Query Sentinel-2 data\n",
    "\n",
    "Retrieve Sentinel-2 data over the Mount Shasta area by specifying the data location, time range, and property filters in the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a72b2ec9-f471-4f92-b442-78eff8ceba57",
   "metadata": {},
   "outputs": [],
   "source": [
    "search_rdc_args = {\n",
    "    \"Arn\": \"arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8\",  # sentinel-2 L2A COG\n",
    "    \"RasterDataCollectionQuery\": {\n",
    "        \"AreaOfInterest\": {\n",
    "            \"AreaOfInterestGeometry\": {\n",
    "                \"PolygonGeometry\": {\n",
    "                    \"Coordinates\": [\n",
    "                        [\n",
    "                            [-122.198, 41.407],\n",
    "                            [-122.191, 41.407],\n",
    "                            [-122.191, 41.411],\n",
    "                            [-122.198, 41.411],\n",
    "                            [-122.198, 41.407],\n",
    "                        ]\n",
    "                    ]\n",
    "                }\n",
    "            }\n",
    "        },\n",
    "        \"TimeRangeFilter\": {\n",
    "            \"StartTime\": \"2017-07-01T00:00:00Z\",\n",
    "            \"EndTime\": \"2022-07-10T23:59:59Z\",\n",
    "        },\n",
    "        \"PropertyFilters\": {\n",
    "            \"Properties\": [{\"Property\": {\"EoCloudCover\": {\"LowerBound\": 0, \"UpperBound\": 1}}}],\n",
    "            \"LogicalOperator\": \"AND\",\n",
    "        },\n",
    "        \"BandFilter\": [\"visual\"],\n",
    "    },\n",
    "}\n",
    "\n",
    "tci_urls = []\n",
    "data_manifests = []\n",
    "while search_rdc_args.get(\"NextToken\", True):\n",
    "    search_result = sg_client.search_raster_data_collection(**search_rdc_args)\n",
    "    if search_result.get(\"NextToken\"):\n",
    "        data_manifests.append(search_result)\n",
    "    for item in search_result[\"Items\"]:\n",
    "        tci_url = item[\"Assets\"][\"visual\"][\"Href\"]\n",
    "        # print(tci_url)\n",
    "        tci_urls.append(tci_url)\n",
    "\n",
    "    search_rdc_args[\"NextToken\"] = search_result.get(\"NextToken\")\n",
    "\n",
    "print(f\"{len(tci_urls)} images found.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23788ee9-bf5f-48b8-80f9-e0432048f679",
   "metadata": {},
   "source": [
    "<a id='3'></a>\n",
    "## Start an Earth Observation Job (EOJ) to identify the land cover types in the area of Mount Shasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a7f45ab-06fc-41cc-b7c2-9c169ea5b787",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform land cover segmentation on images returned from the sentinel dataset.\n",
    "eoj_input_config = {\n",
    "    \"RasterDataCollectionQuery\": {\n",
    "        \"RasterDataCollectionArn\": \"arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/nmqj48dcu3g7ayw8\",\n",
    "        \"AreaOfInterest\": {\n",
    "            \"AreaOfInterestGeometry\": {\n",
    "                \"PolygonGeometry\": {\n",
    "                    \"Coordinates\": [\n",
    "                        [\n",
    "                            [-122.198, 41.407],\n",
    "                            [-122.191, 41.407],\n",
    "                            [-122.191, 41.411],\n",
    "                            [-122.198, 41.411],\n",
    "                            [-122.198, 41.407],\n",
    "                        ]\n",
    "                    ]\n",
    "                }\n",
    "            }\n",
    "        },\n",
    "        \"TimeRangeFilter\": {\n",
    "            \"StartTime\": \"2017-07-01T00:00:00Z\",\n",
    "            \"EndTime\": \"2022-07-10T23:59:59Z\",\n",
    "        },\n",
    "        \"PropertyFilters\": {\n",
    "            \"Properties\": [{\"Property\": {\"EoCloudCover\": {\"LowerBound\": 0, \"UpperBound\": 1}}}],\n",
    "            \"LogicalOperator\": \"AND\",\n",
    "        },\n",
    "    }\n",
    "}\n",
    "eoj_config = {\"LandCoverSegmentationConfig\": {}}\n",
    "\n",
    "response = sg_client.start_earth_observation_job(\n",
    "    Name=\"mont-shasta-landcover\",\n",
    "    InputConfig=eoj_input_config,\n",
    "    JobConfig=eoj_config,\n",
    "    ExecutionRoleArn=execution_role,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47e00573-a20b-4264-8fb7-f80a1c92f4e3",
   "metadata": {},
   "source": [
    "Monitor the EOJ status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "309b74ca-883e-4b64-b597-0f0e154a23dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "eoj_arn = response[\"Arn\"]\n",
    "job_details = sg_client.get_earth_observation_job(Arn=eoj_arn)\n",
    "{k: v for k, v in job_details.items() if k in [\"Arn\", \"Status\", \"DurationInSeconds\"]}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec7662d-7528-4cb8-99e4-dedc4025edca",
   "metadata": {},
   "source": [
    "<a id='4'></a>\n",
    "\n",
    "## Visualize EOJ inputs and outputs in FourSquare Studio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2749fd0-a42b-43d6-9d32-39ed7d1dbc5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates an instance of the map to add EOJ input/ouput layer\n",
    "map = sagemaker_geospatial_map.create_map({\"is_raster\": True})\n",
    "map.set_sagemaker_geospatial_client(sg_client)\n",
    "# Render the map\n",
    "map.render()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eca43da-4bdc-4133-bea2-45bbf4fcc269",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize AOI\n",
    "config = {\"label\": \"Mount Shasta AOI\"}\n",
    "aoi_layer = map.visualize_eoj_aoi(Arn=eoj_arn, config=config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f82107b8-eefe-424e-8389-d11d31acc772",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize input after EOJ has completed.\n",
    "time_range_filter = {\"start_date\": \"2019-07-01T00:00:00Z\", \"end_date\": \"2019-07-10T23:59:59Z\"}\n",
    "config = {\"label\": \"Input\"}\n",
    "input_layer = map.visualize_eoj_input(\n",
    "    Arn=eoj_arn, config=config, time_range_filter=time_range_filter\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebf1b5d7-6e15-4c7f-9dcd-e98e02d42b34",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize output, EOJ needs to be in completed status.\n",
    "time_range_filter = {\"start_date\": \"2019-07-01T00:00:00Z\", \"end_date\": \"2019-07-10T23:59:59Z\"}\n",
    "config = {\"preset\": \"singleBand\", \"band_name\": \"mask\"}\n",
    "output_layer = map.visualize_eoj_output(\n",
    "    Arn=eoj_arn, config=config, time_range_filter=time_range_filter\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "321c260b-49fd-4973-9dd6-0f00f629fab6",
   "metadata": {},
   "source": [
    "<a id='5'></a>\n",
    "## Extract the land cover segmentation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d68d8237-c580-4164-b69c-6098eb76fd00",
   "metadata": {},
   "outputs": [],
   "source": [
    "sagemaker_session = sagemaker.Session()\n",
    "s3_bucket_name = sagemaker_session.default_bucket()  # Replace with your own bucket if needed\n",
    "s3_bucket = session.resource(\"s3\").Bucket(s3_bucket_name)\n",
    "prefix = \"eoj_montshasta\"  # Replace with the S3 prefix desired\n",
    "export_bucket_and_key = f\"s3://{s3_bucket_name}/{prefix}/\"\n",
    "\n",
    "eoj_output_config = {\"S3Data\": {\"S3Uri\": export_bucket_and_key}}\n",
    "export_response = sg_client.export_earth_observation_job(\n",
    "    Arn=eoj_arn,\n",
    "    ExecutionRoleArn=execution_role,\n",
    "    OutputConfig=eoj_output_config,\n",
    "    ExportSourceImages=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36eea163-528c-4471-9a4e-2cc71b55af70",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monitor the export job status\n",
    "export_job_details = sg_client.get_earth_observation_job(Arn=export_response[\"Arn\"])\n",
    "{k: v for k, v in export_job_details.items() if k in [\"Arn\", \"ExportStatus\"]}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a667db5f-ac1c-4218-9e28-48208d35d54a",
   "metadata": {},
   "source": [
    "<a id='6'></a>\n",
    "## Measure snow cover area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef9344c-9f0b-42a8-a1e9-7c77c57487c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from glob import glob\n",
    "import cv2\n",
    "import numpy as np\n",
    "import tifffile\n",
    "import matplotlib.pyplot as plt\n",
    "from urllib.parse import urlparse\n",
    "from botocore import UNSIGNED\n",
    "from botocore.config import Config\n",
    "\n",
    "# Download land cover masks\n",
    "mask_dir = \"./masks/mont_shasta\"\n",
    "os.makedirs(mask_dir, exist_ok=True)\n",
    "image_paths = []\n",
    "for s3_object in s3_bucket.objects.filter(Prefix=prefix).all():\n",
    "    path, filename = os.path.split(s3_object.key)\n",
    "    if (filename.endswith(\".tif\")) and (700 < int(filename.split(\"_\")[2][5:8]) < 708):\n",
    "        mask_name = mask_dir + \"/\" + filename\n",
    "        s3_bucket.download_file(s3_object.key, mask_name)\n",
    "        print(\"Downloaded mask: \" + mask_name)\n",
    "\n",
    "# Download source images for visualization\n",
    "image_dir = \"./images/mont_shasta\"\n",
    "os.makedirs(image_dir, exist_ok=True)\n",
    "for tci_url in tci_urls:\n",
    "    url_parts = urlparse(tci_url)\n",
    "    img_id = url_parts.path.split(\"/\")[-2]\n",
    "    # Only use data from the 1st week of July\n",
    "    if 700 < int(img_id.split(\"_\")[2][5:8]) < 708:\n",
    "        tci_download_path = image_dir + \"/\" + img_id + \"_TCI.tif\"\n",
    "        cogs_bucket = session.resource(\n",
    "            \"s3\", config=Config(signature_version=UNSIGNED, region_name=\"us-west-2\")\n",
    "        ).Bucket(url_parts.hostname.split(\".\")[0])\n",
    "        cogs_bucket.download_file(url_parts.path[1:], tci_download_path)\n",
    "        print(\"Downloaded image: \" + tci_download_path)\n",
    "\n",
    "print(\"Downloads complete.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1e76bc-6204-44b7-9ee4-4d618cfb2943",
   "metadata": {},
   "source": [
    "Extract the snow mask and measure the snow cover area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23111adf-e639-45eb-baa1-6f9c20ae47c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_files = glob(image_dir + \"/*.tif\")\n",
    "mask_files = glob(mask_dir + \"/*.tif\")\n",
    "image_files.sort(key=lambda x: x.split(\"_\")[3])\n",
    "mask_files.sort(key=lambda x: x.split(\"_\")[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bd68e7f-4dee-4f94-b699-f8ce91c1c76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "overlay_dir = \"./masks/mont_shasta_overlay\"\n",
    "os.makedirs(overlay_dir, exist_ok=True)\n",
    "snow_areas = []\n",
    "mask_dates = []\n",
    "\n",
    "for image_file, mask_file in zip(image_files, mask_files):\n",
    "    image_id = image_file.split(\"/\")[-1].split(\"_TCI\")[0]\n",
    "    mask_id = mask_file.split(\"/\")[-1].split(\".tif\")[0]\n",
    "    mask_date = mask_id.split(\"_\")[2]\n",
    "    mask_dates.append(mask_date)\n",
    "    assert image_id == mask_id\n",
    "    image = tifffile.imread(image_file)\n",
    "    image_ds = cv2.resize(image, (1830, 1830), interpolation=cv2.INTER_LINEAR)\n",
    "    snow_mount_img = image_ds[100:400, 1000:1250]  # a crop of the Mount Shasta area\n",
    "    mask = tifffile.imread(mask_file)\n",
    "    snow_mask = np.isin(mask, [11]).astype(np.uint8)  # extract snow mask\n",
    "    snow_mount_mask = snow_mask[100:400, 1000:1250]\n",
    "    snow_area = snow_mount_mask.sum() * 60 * 60 / (1000 * 1000)  # calculate the snow cover area\n",
    "    snow_areas.append(snow_area)\n",
    "    red_img = np.zeros(snow_mount_img.shape, snow_mount_img.dtype)\n",
    "    red_img[:, :] = (255, 0, 0)\n",
    "    red_mask = cv2.bitwise_and(red_img, red_img, mask=snow_mount_mask)\n",
    "    overlay_img = cv2.addWeighted(red_mask, 1, snow_mount_img, 1, 0, snow_mount_img)\n",
    "    cv2.putText(\n",
    "        overlay_img,\n",
    "        f\"{mask_date}\",\n",
    "        (0, 20),\n",
    "        cv2.FONT_HERSHEY_SIMPLEX,\n",
    "        0.5,\n",
    "        (255, 0, 0),\n",
    "        1,\n",
    "        cv2.LINE_AA,\n",
    "    )\n",
    "    cv2.putText(\n",
    "        overlay_img,\n",
    "        f\"{snow_area} [sq km]\",\n",
    "        (0, 40),\n",
    "        cv2.FONT_HERSHEY_SIMPLEX,\n",
    "        0.5,\n",
    "        (255, 0, 0),\n",
    "        1,\n",
    "        cv2.LINE_AA,\n",
    "    )\n",
    "    overlay_file = overlay_dir + \"/\" + mask_date + \".png\"\n",
    "    cv2.imwrite(overlay_file, cv2.cvtColor(overlay_img, cv2.COLOR_RGB2BGR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18fdc809-02bf-49fc-90b2-5d4bb1269983",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average snow area in July\n",
    "snow_area_date = {}\n",
    "for snow_area, mask_date in zip(snow_areas, mask_dates):\n",
    "    mask_date = mask_date[:6]  # year-month\n",
    "    if mask_date in snow_area_date:\n",
    "        snow_area_date[mask_date] += np.array([snow_area, 1])\n",
    "    else:\n",
    "        snow_area_date[mask_date] = np.array([snow_area, 1])\n",
    "dates = list(snow_area_date.keys())\n",
    "snow_area_year = np.zeros(len(dates))\n",
    "for i, snow_area in enumerate(snow_area_date.values()):\n",
    "    snow_area_year[i] = snow_area[0] / snow_area[1]\n",
    "\n",
    "# Plot snow cover area vs. time.\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.title(\"Mount Shasta snow area from 2017 to 2022.\", fontsize=20)\n",
    "plt.xticks(rotation=45)\n",
    "plt.ylabel(\"Snow cover area [sq km]\", fontsize=14)\n",
    "plt.plot(dates, snow_area_year, marker=\"o\")\n",
    "for i, v in enumerate(snow_area_year):\n",
    "    plt.text(i, v + 1.5, \"%d\" % v, ha=\"center\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320a3125-7622-4ef8-a674-ee220a32598c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio.v2 as imageio\n",
    "\n",
    "frames = []\n",
    "filenames = glob(overlay_dir + \"/*.png\")\n",
    "filenames.sort()\n",
    "\n",
    "for filename in filenames:\n",
    "    frames.append(imageio.imread(filename))\n",
    "imageio.mimsave(\"mont_shasta.gif\", frames, duration=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07542687-5288-4b6f-ae4e-483af2f130a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "\n",
    "HTML('<img src=\"./mont_shasta.gif\" width=\"400\">')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91a37c38-b257-4a52-a98f-312d6f21ec4f",
   "metadata": {},
   "source": [
    "<a id='7'></a>\n",
    "\n",
    "## Query Landsat 8 data\n",
    "Retrieve the Landsat surface temperature data to gain a better understanding of the changes that have taken place in the same area over the past 5 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b3a9a8f-1efb-4dd2-8d35-cd792ed83e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "search_rdc_args = {\n",
    "    \"Arn\": \"arn:aws:sagemaker-geospatial:us-west-2:378778860802:raster-data-collection/public/gmqa64dcu2g9ayx1\",  # Landsat 8 C2L2\n",
    "    \"RasterDataCollectionQuery\": {\n",
    "        \"AreaOfInterest\": {\n",
    "            \"AreaOfInterestGeometry\": {\n",
    "                \"PolygonGeometry\": {\n",
    "                    \"Coordinates\": [\n",
    "                        [\n",
    "                            [-122.198, 41.407],\n",
    "                            [-122.191, 41.407],\n",
    "                            [-122.191, 41.411],\n",
    "                            [-122.198, 41.411],\n",
    "                            [-122.198, 41.407],\n",
    "                        ]\n",
    "                    ]\n",
    "                }\n",
    "            }\n",
    "        },\n",
    "        \"TimeRangeFilter\": {\n",
    "            \"StartTime\": \"2017-07-01T00:00:00Z\",\n",
    "            \"EndTime\": \"2022-07-30T23:59:59Z\",\n",
    "        },\n",
    "        \"PropertyFilters\": {\n",
    "            \"Properties\": [{\"Property\": {\"EoCloudCover\": {\"LowerBound\": 0, \"UpperBound\": 1}}}],\n",
    "            \"LogicalOperator\": \"AND\",\n",
    "        },\n",
    "        \"BandFilter\": [\"lwir11\"],\n",
    "    },\n",
    "}\n",
    "\n",
    "lst_urls = []\n",
    "data_manifests = []\n",
    "while search_rdc_args.get(\"NextToken\", True):\n",
    "    search_result = sg_client.search_raster_data_collection(**search_rdc_args)\n",
    "    if search_result.get(\"NextToken\"):\n",
    "        data_manifests.append(search_result)\n",
    "    for item in search_result[\"Items\"]:\n",
    "        lst_url = item[\"Assets\"][\"lwir11\"][\"Href\"]\n",
    "        print(lst_url)\n",
    "        lst_urls.append(lst_url)\n",
    "\n",
    "    search_rdc_args[\"NextToken\"] = search_result.get(\"NextToken\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71bdffb7-87b6-4ed3-9987-3b93703f0d61",
   "metadata": {},
   "source": [
    "Download Landsat surface temperature data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0317bad-aeb0-4d61-8fd1-bd3ff8404b20",
   "metadata": {},
   "outputs": [],
   "source": [
    "s3 = boto3.client(\"s3\")\n",
    "lst_dir = \"./images/mont_shasta_lst\"\n",
    "os.makedirs(lst_dir, exist_ok=True)\n",
    "lst_paths = []\n",
    "\n",
    "for lst_url in lst_urls:\n",
    "    url_parts = urlparse(lst_url)\n",
    "    img_id = url_parts.path.split(\"/\")[-2]\n",
    "    if 700 < int(img_id[22:25]) < 731:\n",
    "        lst_download_path = lst_dir + \"/\" + img_id + \"_ST.tif\"\n",
    "        s3.download_file(\n",
    "            url_parts.netloc,\n",
    "            url_parts.path[1:],\n",
    "            lst_download_path,\n",
    "            ExtraArgs={\"RequestPayer\": \"requester\"},\n",
    "        )\n",
    "        print(\"Downloaded LST image: \" + img_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47b71c3a-5430-4fc7-a22d-1e100cbd5bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_files = glob(lst_dir + \"/*.tif\")\n",
    "lst_files.sort(key=lambda x: x.split(\"_\")[5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb160604-f2e4-4d82-bb9e-963593686215",
   "metadata": {},
   "source": [
    "Landsat surface temperature: https://www.usgs.gov/faqs/how-do-i-use-scale-factor-landsat-level-2-science-products\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8643f7f4-9cad-40ad-94c3-b4fef0e67e36",
   "metadata": {},
   "outputs": [],
   "source": [
    "mont_lsts = []\n",
    "lst_dates = []\n",
    "\n",
    "for lst_file in lst_files:\n",
    "    lst_date = lst_file.split(\"_\")[5]\n",
    "    lst_dates.append(lst_date)\n",
    "    lst = tifffile.imread(lst_file)\n",
    "    mont_lst = lst[5000:5500, 3900:4400]  # bounding box for Mount Shasta\n",
    "    mont_lst_c = (np.mean(mont_lst) * 0.00341802 + 149) - 273  # surface temperature in Celsius\n",
    "    mont_lsts.append(mont_lst_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf2900d2-44f5-4c7b-bb25-2a42f4800d58",
   "metadata": {},
   "source": [
    "<a id='8'></a>\n",
    "\n",
    "## Relationship between surface temperature and snow cover area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc07242a-9c49-4917-83b5-06686a16880e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average surface temperature in July\n",
    "mont_temp_date = {}\n",
    "for mont_lst, lst_date in zip(mont_lsts, lst_dates):\n",
    "    lst_date = lst_date[:6]  # year-month\n",
    "    if lst_date in mont_temp_date:\n",
    "        mont_temp_date[lst_date] += np.array([mont_lst, 1])\n",
    "    else:\n",
    "        mont_temp_date[lst_date] = np.array([mont_lst, 1])\n",
    "dates = list(mont_temp_date.keys())\n",
    "mont_temp_year = np.zeros(len(dates))\n",
    "for i, mont_temp in enumerate(mont_temp_date.values()):\n",
    "    mont_temp_year[i] = mont_temp[0] / mont_temp[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a26ab28d-16cb-4b39-8b1c-0ed8ab565054",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import pearsonr\n",
    "\n",
    "# Pearson's r to measure the strength and direction of the relationship between two variables\n",
    "res = pearsonr(snow_area_year, mont_temp_year)\n",
    "print(\"Pearson's correlction coefficient: %.2f\" % res.statistic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f8cdef3-28ad-4374-ba18-413f8393326a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(8, 4))\n",
    "ax2 = ax1.twinx()\n",
    "lns1 = ax1.plot(dates, snow_area_year, \"b-\", marker=\"o\", label=\"Snow cover area\")\n",
    "for i, v in enumerate(snow_area_year):\n",
    "    ax1.text(i, v - 5, \"%d\" % v, ha=\"center\")\n",
    "lns2 = ax2.plot(dates, mont_temp_year, \"r-\", marker=\"x\", label=\"Surface temperature\")\n",
    "for i, v in enumerate(mont_temp_year):\n",
    "    ax2.text(i, v + 0.2, \"%d\" % v, ha=\"center\")\n",
    "ax1.set_ylabel(\"Snow cover area [sq km]\", fontsize=12)\n",
    "ax2.set_ylabel(\"Surface temperature [\\N{DEGREE SIGN}C]\", fontsize=12)\n",
    "lns = lns1 + lns2\n",
    "labs = [l.get_label() for l in lns]\n",
    "ax1.legend(lns, labs, loc=7)\n",
    "plt.title(\"Pearsons correlation: %.2f\" % res.statistic)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook CI Test Results\n",
    "\n",
    "This notebook was tested in multiple regions. The test results are as follows, except for us-west-2 which is shown at the top of the notebook.\n",
    "\n",
    "![This us-east-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/us-east-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This us-east-2 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/us-east-2/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This us-west-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/us-west-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ca-central-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ca-central-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This sa-east-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/sa-east-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This eu-west-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/eu-west-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This eu-west-2 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/eu-west-2/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This eu-west-3 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/eu-west-3/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This eu-central-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/eu-central-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This eu-north-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/eu-north-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ap-southeast-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ap-southeast-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ap-southeast-2 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ap-southeast-2/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ap-northeast-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ap-northeast-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ap-northeast-2 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ap-northeast-2/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n",
    "\n",
    "![This ap-south-1 badge failed to load. Check your device's internet connectivity, otherwise the service is currently unavailable](https://prod.us-west-2.tcx-beacon.docs.aws.dev/sagemaker-nb/ap-south-1/sagemaker-geospatial|mount-shasta-glacier-melting-monitoring|mount_shasta_glacier_melt_monitoring.ipynb)\n"
   ]
  }
 ],
 "metadata": {
  "instance_type": "ml.m5.4xlarge",
  "kernelspec": {
   "display_name": "Python 3 (Geospatial 1.0)",
   "language": "python",
   "name": "python3__SAGEMAKER_INTERNAL__arn:aws:sagemaker:us-west-2:081189585635:image/sagemaker-geospatial-v1-0"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  },
  "toc-autonumbering": true,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}