python-data/exercises/ex09b_satellite_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "73b81a5a-4fc6-4c33-849b-3b717a43b1c8",
   "metadata": {},
   "source": [
    "# Exercise 9b: Working with Satellite Data\n",
    "\n",
    "## Aim: Use python tools to search for, download, and manipulate satellite data\n",
    "\n",
    "### Issues covered:\n",
    "- Search for and request data from a public STAC catalogue of satellite imagery\n",
    "- Download satellite imagery as raster data \n",
    "- Read rasters into python using the rioxarray package\n",
    "- Visualise single/multi-band raster data\n",
    "\n",
    "### Introduction\n",
    "\n",
    "A number of satellites take snapshots of the Earth’s surface from space. The images recorded by these remote sensors represent a very precious data source for any activity that involves monitoring changes on Earth. Satellite imagery is typically provided in the form of geospatial raster data, with the measurements in each grid cell (“pixel”) being associated to accurate geographic coordinate information.\n",
    "\n",
    "In this notebook exercise we will explore how to access open satellite data using Python. In particular, we will consider [the Sentinel-2 data collection that is hosted on AWS](https://registry.opendata.aws/sentinel-2-l2a-cogs). This dataset consists of multi-band optical images acquired by the two satellites of [the Sentinel-2 mission](https://sentinel.esa.int/web/sentinel/missions/sentinel-2) and it is continuously updated with new images.\n",
    "\n",
    "\n",
    "# 1. Search for satellite imagery\n",
    "\n",
    "**The SpatioTemporal Asset Catalog (STAC) specification**\n",
    "\n",
    "Current sensor resolutions and satellite revisit periods are such that terabytes of data products are added daily to the corresponding collections. Such datasets cannot be made accessible to users via full-catalog download. Space agencies and other data providers often offer access to their data catalogs through interactive Graphical User Interfaces (GUIs), see for instance the [Copernicus Open Access Hub portal](https://scihub.copernicus.eu/dhus/#/home) for the Sentinel missions. Accessing data via a GUI is a nice way to explore a catalog and get familiar with its content, but it represents a heavy and error-prone task that should be avoided if carried out systematically to retrieve data.\n",
    "\n",
    "A service that offers programmatic access to the data enables users to reach the desired data in a more reliable, scalable and reproducible manner. An important element in the software interface exposed to the users, which is generally called the Application Programming Interface (API), is the use of standards. Standards, in fact, can significantly facilitate the reusability of tools and scripts across datasets and applications.\n",
    "\n",
    "The SpatioTemporal Asset Catalog (STAC) specification is an emerging standard for describing geospatial data. By organizing metadata in a form that adheres to the STAC specifications, data providers make it possible for users to access data from different missions, instruments and collections using the same set of tools.\n",
    "\n",
    "\n",
    "![Views of the STAC browser](https://carpentries-incubator.github.io/geospatial-python/fig/E05/STAC-browser.jpg)\n",
    "Views of the radiant earth STAC browser\n",
    "\n",
    "## More Resources on STAC\n",
    "- [STAC specification](https://github.com/radiantearth/stac-spec#readme)\n",
    "- [Tools based on STAC](https://stacindex.org/ecosystem)\n",
    "- [STAC catalogs](https://stacindex.org/catalogs)\n",
    "\n",
    "## Search a STAC catalog\n",
    "\n",
    "The [STAC browser](https://radiantearth.github.io/stac-browser/#/) is a good starting point to discover available datasets, as it provides an up-to-date list of existing STAC catalogs. From the list, let's click on the \"Earth Search\" catalog, i.e. the access point to search the archive of Sentinel-2 images hosted on AWS.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb95dfdf-8721-4d47-af20-211e2f7bd491",
   "metadata": {},
   "source": [
    "## Install some packages we will need\n",
    "\n",
    "We need to install some additional python packages which unfortunately aren't (yet) on Jaspy. To do this, we run:\n",
    "\n",
    "`!pip install --user pystac_client`\n",
    "\n",
    "Which will install these python packages into your local python path, so we can use them with your account alongside all the packages in Jaspy. __NOTE: this command may take some time__ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4ab589c7-f462-4c6d-aaa6-b1dc630e5cf6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Type the pip command here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "517be10e-1c03-433c-b6b9-3722cc0d15b9",
   "metadata": {},
   "source": [
    "## **Exercise:** Discover a STAC catalog\n",
    "Let's take a moment to explore the Earth Search STAC catalog, which is the catalog indexing the Sentinel-2 collection\n",
    "that is hosted on AWS. We can interactively browse this catalog using the STAC browser at [this link](https://radiantearth.github.io/stac-browser/#/external/earth-search.aws.element84.com/v1).\n",
    "\n",
    "1. Open the link in your web browser. Which (sub-)catalogs are available?\n",
    "2. Open the Sentinel-2 Level 2A collection, and select one item from the list. Each item corresponds to a satellite\n",
    "\"scene\", i.e. a portion of the footage recorded by the satellite at a given time. Have a look at the metadata fields\n",
    "and the list of assets. What kind of data do the assets represent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c869fc6-9581-4ad2-9a37-72d1d61b84d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try something in here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "462ea61f-3cc0-4793-be7a-b7f9886e1484",
   "metadata": {},
   "source": [
    "When opening a catalog with the STAC browser, you can access the API URL by clicking on the \"Source\" button on the top\n",
    "right of the page. By using this URL, we have access to the catalog content and, if supported by the catalog, to the\n",
    "functionality of searching its items. For the Earth Search STAC catalog the API URL is:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca27bf8b-05c4-4d6b-b8ee-d152487e6f06",
   "metadata": {},
   "outputs": [],
   "source": [
    "api_url = \"https://earth-search.aws.element84.com/v1\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d298328b-2f37-442e-996f-ea61c68eb039",
   "metadata": {},
   "source": [
    "You can query a STAC API endpoint from Python using the `pystac_client` library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cac729fb-8cda-484b-9fd4-da68a2c8267c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystac_client import Client\n",
    "\n",
    "client = Client.open(api_url)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81f1cc52-814e-4af0-908f-6d4aa7cc10fe",
   "metadata": {},
   "source": [
    "In the following, we ask for scenes belonging to the `sentinel-2-l2a` collection. This dataset includes Sentinel-2 data products pre-processed at level 2A (bottom-of-atmosphere reflectance) and saved in Cloud Optimized GeoTIFF (COG) format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e7e407d-721e-4187-8a47-e8e9634262c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "collection = \"sentinel-2-l2a\"  # Sentinel-2, Level 2A, Cloud Optimized GeoTiffs (COGs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f5d59b-f60b-45d0-b191-f5e4d5e8939e",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## A note on cloud-optimized GeoTIFFs\n",
    "\n",
    "Cloud Optimized GeoTIFFs (COGs) are regular GeoTIFF files with some additional features that make them ideal to be employed in the context of cloud computing and other web-based services. This format builds on the widely-employed GeoTIFF format. In short, a GeoTIFF is a standard .tif image format with additional spatial (georeferencing) information embedded in the file as tags. These tags should include the following raster metadata:\n",
    "- Extent\n",
    "- Resolution\n",
    "- Coordinate Reference System (CRS)\n",
    "- Values that represent missing data (NoDataValue)\n",
    "\n",
    "COGs, by extension, are regular GeoTIFF files with a special internal structure. One of the features of COGs is that data is organized in \"blocks\" that can be accessed remotely via independent HTTP requests. Data users can thus access the only blocks of a GeoTIFF that are relevant for their analysis, without having to download the full file. In addition, COGs typically include multiple lower-resolution versions of the original image, called \"overviews\", which can also be accessed independently. By providing this \"pyramidal\" structure, users that are not interested in the details provided by a high-resolution raster can directly access the lower-resolution versions of the same image, significantly saving on the downloading time. More information on the COG format can be found [here](https://www.cogeo.org).\n",
    "\n",
    "---\n",
    "\n",
    "We also ask for scenes intersecting a geometry defined using the `shapely` library (in this case, a point):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b58e48e4-2609-4e86-b434-ed664dafa6f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geometry import Point\n",
    "point = Point(4.89, 52.37)  # AMS coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b597f5-db31-4bad-b858-59e9f2961d92",
   "metadata": {},
   "source": [
    "Note: at this stage, we are only dealing with metadata, so no image is going to be downloaded yet. But even metadata can be quite bulky if a large number of scenes match our search! For this reason, we limit the search result to 10 items:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08d64392-75f9-442c-bcc5-e173e1448683",
   "metadata": {},
   "outputs": [],
   "source": [
    "search = client.search(\n",
    "    collections=[collection],\n",
    "    intersects=point,\n",
    "    max_items=10,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71b0aa02-8e8c-44c2-974e-b49ef84c7163",
   "metadata": {},
   "source": [
    "We submit the query and find out how many scenes match our search criteria (please note that this output can be different as more data is added to the catalog):\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd8d411-a898-452a-be69-f19a8cdb920c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(search.matched())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa2aeb2-0f3a-4a5c-bd2a-c8e1e5b704e3",
   "metadata": {},
   "source": [
    "Finally, we retrieve the metadata of the search results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2790d193-02a1-42e2-a51c-42adfc17041c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "items = search.item_collection()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "500a6faf-0d21-421e-b8b3-60b0d0df6e64",
   "metadata": {},
   "source": [
    "The variable `items` is an `ItemCollection` object. We can check its size by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb1f4c0-104a-4c38-953d-f5b1f2117643",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(items))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2ccac65-d1aa-4e66-8243-b4d27968d638",
   "metadata": {},
   "source": [
    "which is consistent with the maximum number of items that we have set in the search criteria. We can iterate over the returned items and print these to show their IDs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c94ff82-849e-4164-89b2-8fb6bec6d622",
   "metadata": {},
   "outputs": [],
   "source": [
    "for item in items:\n",
    "    print(item)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbddf824-83d3-4789-a354-f1989936c438",
   "metadata": {},
   "source": [
    "Each of the items contains information about the scene geometry, its acquisition time, and other metadata that can be accessed as a dictionary from the `properties` attribute.\n",
    "\n",
    "Let's inspect the metadata associated with the first item of the search results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f32b169d-31df-4081-ba4d-9219f4efeb15",
   "metadata": {},
   "outputs": [],
   "source": [
    "item = items[0]\n",
    "print(item.datetime)\n",
    "print(item.geometry)\n",
    "print(item.properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad300420-d4d9-4765-8d2c-1aa961ae77d3",
   "metadata": {},
   "source": [
    "## **Exercise**: Search satellite scenes using metadata filters\n",
    "Search for all the available Sentinel-2 scenes in the `Sentinel-2-Level-2A` collection that satisfy the following criteria:\n",
    "- Intersect a provided bounding box, use ±0.01 deg in lat/lon from the previously defined point (hint: use `.buffer().bounds` methods on the shapely `point` object we created above). Give this as the `bbox` argument of `client.search()`.\n",
    "- Have been recorded between 20 March 2020 and 30 March 2020 (hint: use `client.search()` with an argument `datetime=\"date1/date2\"`);\n",
    "- Have a cloud coverage smaller than 15% (hint: use the `query` argument of `client.search` - info can be found [here](https://pystac-client.readthedocs.io/en/latest/usage.html#query-extension)).\n",
    "- Once you have added these as arguments to `client.search()` print out `.matched()` as we did early to see our number of matches.\n",
    "- Retrieve the metadata of the item as we did earlier.\n",
    "\n",
    "How many scenes are available? Save the search results in GeoJSON format as `search.json`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7b4d17e-2746-4546-9c83-a7c43342ce9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try something in here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b0566d9-fc33-48b9-9ca0-afa8a889add0",
   "metadata": {},
   "source": [
    "## Access the assets\n",
    "\n",
    "So far we have only discussed metadata - but how can one get to the actual images of a satellite scene (the \"assets\" in the STAC nomenclature)? These can be reached via links that are made available through the item's attribute `assets`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18778b46-71fd-4f21-872f-4112a5656439",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "assets = items[0].assets  # first item's asset dictionary\n",
    "\n",
    "# Have a look at the keys\n",
    "print(assets.keys())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "664aa928-bf54-4003-833e-a7e93bab27a7",
   "metadata": {
    "tags": []
   },
   "source": [
    "We can print a minimal description of the available assets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4f2cb8-c392-4d43-b298-824529df3131",
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, asset in assets.items():\n",
    "    print(f\"{key}: {asset.title}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9012da82-cb57-455e-96ac-5c6df9eeb3ae",
   "metadata": {},
   "source": [
    "Among the others, assets include multiple raster data files (one per optical band, as acquired by the multi-spectral instrument), a thumbnail, a true-color image (\"visual\"), instrument metadata and scene-classification information (\"SCL\"). Let's get the URL links to the actual asset:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20a37bef-2d6e-41d8-acbf-3eb25fa88581",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(assets[\"thumbnail\"].href)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae6b336c-18b6-48cb-8ded-309b0287bdf3",
   "metadata": {
    "tags": []
   },
   "source": [
    "This can be used to download the corresponding file:\n",
    "\n",
    "![Overview of the true-colour image](../images/sentinel-image.png)\n",
    "\n",
    "###### Overview of the true-colour image (\"thumbnail\")\n",
    "\n",
    "\n",
    "Remote raster data can be directly opened via the `rioxarray` library. Raster data is a type of digital data representation used primary in GIS and remote sensing. It represents spatial information as a grid of cells/pixels.\n",
    "\n",
    "We will\n",
    "learn more about this library in the next part of the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4b25f88-f1a8-4e53-8c1e-823e83cb4549",
   "metadata": {},
   "source": [
    "`rioxarray` is an extension of the `xarray` library that provides geospatial raster data capabilities. It bridges the gap between `xarray` and `rasterio` allowing users to work with multi-dimensional arrays that include geospatial metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a78e28a3-2c80-41b8-be30-47941de2186b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rioxarray\n",
    "# Get the URL for the Near Infrared (NIR) band.\n",
    "nir_href = assets[\"nir\"].href\n",
    "# Open the raster file and load it into an xarray data array\n",
    "nir = rioxarray.open_rasterio(nir_href)\n",
    "print(nir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf3b34f-be9b-48ed-8647-0a592763311b",
   "metadata": {},
   "source": [
    "We can then save the data to disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45965a50-8df3-40f1-b3d4-b23906c84fd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save whole image to disk as a geotiff file\n",
    "# NOTE: This might take a while\n",
    "nir.rio.to_raster(\"data/nir.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52993398-3216-4d12-808e-1b357b12f684",
   "metadata": {},
   "source": [
    "Since that might take a while, given there are over 10000 x 10000 = a hundred million pixels in the 10 meter NIR band, you can take a smaller subset before downloading it. Because the raster is a COG, we can download just what we need!\n",
    "\n",
    "Here, we specify that we want to download the first (and only) band in the tif file, and a slice of the width and height dimensions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa06ed25-cdcc-4efe-bee0-cbe6f37af9b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save portion of an image to disk\n",
    "nir[0,1500:2200,1500:2200].rio.to_raster(\"data/nir_subset.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5b47a49-3ade-40fe-899e-3b328cf78e0f",
   "metadata": {},
   "source": [
    "The difference is 155 Megabytes for the large image vs about 1 Megabyte for the subset.\n",
    "\n",
    "\n",
    "## **Exercise:** Downloading Landsat 8 Assets\n",
    "In this exercise we put in practice all the skills we have learned thusfar to retrieve images from a different mission: [Landsat 8](https://www.usgs.gov/landsat-missions/landsat-8). In particular, we browse images from the [Harmonized Landsat Sentinel-2 (HLS) project](https://lpdaac.usgs.gov/products/hlsl30v002/), which provides images from NASA's Landsat 8 and ESA's Sentinel-2 that have been made consistent with each other. The HLS catalog is indexed in the NASA Common Metadata Repository (CMR) and it can be accessed from the STAC API endpoint at the following URL:\n",
    "`https://cmr.earthdata.nasa.gov/stac/LPCLOUD`.\n",
    "\n",
    "1. Using `pystac_client`, search for all assets of the Landsat 8 collection (`HLSL30.v2.0`) from February to March\n",
    "  2021, intersecting the point with longitude/latitute coordinates (-73.97, 40.78) deg. (Hint: start by connecting to the STAC API endpoint, create a `Client` object, do a `client.search` with the `collections`, `intersects` and `datetime` arguments, retrieve the search results with `.item_collection()` and then print out the length of the resulting data.\n",
    "2. You can sort and select by cloud cover using the following, where `items` is your item collection:\n",
    "```\n",
    "items_sorted = sorted(items, key=lambda x: x.properties[\"eo:cloud_cover\"]) # sorting and then selecting by cloud cover\n",
    "item = items_sorted[0]\n",
    "print(item)\n",
    "```\n",
    "This line is using a lamdba function in the `key`. This is a shorthand way of defining a function without giving it a name, useful for a simple function that we only want to use once.\n",
    "3. Visualize an item's thumbnail using `item.assets[\"browse\"].href`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87816ef7-b896-47ea-a254-9b066434cf4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try something in here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ad7784e-5320-4bc1-a61f-dc714eb6e6d6",
   "metadata": {},
   "source": [
    "## Public catalogs, protected data\n",
    "\n",
    "Publicly accessible catalogs and STAC endpoints do not necessarily imply publicly accessible data. Data providers, in\n",
    "fact, may limit data access to specific infrastructures and/or require authentication. For instance, the NASA CMR STAC\n",
    "endpoint considered in the last exercise offers publicly accessible metadata for the HLS collection, but most of the\n",
    "linked assets are available only for registered users (the thumbnail is publicly accessible).\n",
    "\n",
    "The authentication procedure for datasets with restricted access might differ depending on the data provider. For the\n",
    "NASA CMR, follow these steps in order to access data using Python:\n",
    "\n",
    "* Create a NASA Earthdata login account [here](https://urs.earthdata.nasa.gov);\n",
    "* Set up a netrc file with your credentials, e.g. by using [this script](https://git.earthdata.nasa.gov/projects/LPDUR/repos/daac_data_download_python/browse/EarthdataLoginSetup.py);\n",
    "* Define the following environment variables:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20eb8446-450a-4769-9e94-df0ff0f20373",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"GDAL_HTTP_COOKIEFILE\"] = \"./cookies.txt\"\n",
    "os.environ[\"GDAL_HTTP_COOKIEJAR\"] = \"./cookies.txt\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77abe4c0-fe83-4161-9269-13c9a6591ff4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Key takeaways:\n",
    "\n",
    "Accessing satellite images via the providers' API enables a more reliable and scalable data retrieval.\n",
    "\n",
    " - STAC catalogs can be browsed and searched using the same tools and scripts.\n",
    " - `rioxarray` allows you to open and download remote raster files.\n",
    " \n",
    "---\n",
    "\n",
    "# 2. Read and visualise raster data\n",
    "\n",
    "Next, we introduce the fundamental principles, packages and metadata/raster attributes for working with raster data (pixelated/gridded data) in Python. We will also explore how Python handles missing and bad data values.\n",
    "\n",
    "[`rioxarray`](https://corteva.github.io/rioxarray/stable/) is the Python package we will use throughout the rest of this notebook to work with raster data. It is based on the popular [`rasterio`](https://rasterio.readthedocs.io/en/latest/) package for working with rasters and [`xarray`](https://xarray.pydata.org/en/stable/) for working with multi-dimensional arrays.\n",
    "`rioxarray` extends `xarray` by providing top-level functions (e.g. the `open_rasterio` function to open raster datasets) and by adding a set of methods to the main objects of the `xarray` package (the `Dataset` and the `DataArray`). These additional methods are made available via the `rio` accessor and become available from `xarray` objects after importing `rioxarray`.\n",
    "\n",
    "We will also use the [`pystac`](https://github.com/stac-utils/pystac) package (library for working with STAC in Python) to load rasters from the search results we created in the previous section.\n",
    "\n",
    "### About Raster Data\n",
    "\n",
    "Raster data is any pixelated (or gridded) data where each pixel is associated\n",
    "with a specific geographic location. The value of a pixel can be\n",
    "continuous (e.g. elevation) or categorical (e.g. land use). If this sounds\n",
    "familiar, it is because this data structure is very common: it's how\n",
    "we represent any digital image. A geospatial raster is only different\n",
    "from a digital photo in that it is accompanied by spatial information\n",
    "that connects the data to a particular location. This includes the\n",
    "raster's extent and cell size, the number of rows and columns, and\n",
    "its coordinate reference system (or CRS).\n",
    "\n",
    "![raster-concept](https://carpentries-incubator.github.io/geospatial-python/fig/E01/raster_concept.png)\n",
    "###### Raster Concept (Source: National Ecological Observatory Network (NEON))\n",
    "\n",
    "Some examples of continuous rasters include:\n",
    "\n",
    "1. Precipitation maps.\n",
    "2. Maps of tree height derived from LiDAR data.\n",
    "3. Elevation values for a region.\n",
    "\n",
    "A map of elevation for Harvard Forest derived from the [NEON AOP LiDAR sensor](https://www.neonscience.org/data-collection/airborne-remote-sensing)\n",
    "is below. Elevation is represented as a continuous numeric variable in this map. The legend\n",
    "shows the continuous range of values in the data from around 300 to 420 meters.\n",
    "\n",
    "![elevation plot](https://carpentries-incubator.github.io/geospatial-python/fig/E01/continuous-elevation-HARV-plot-01.png)\n",
    "###### Continuous Elevation Map: HARV Field Site\n",
    "\n",
    "\n",
    "## Load a Raster and View Attributes\n",
    "In the previous episode, we searched for Sentinel-2 images, and then saved the search results to a file: `data/search.json`. This contains the information on where and how to access the target images from a remote repository. We can use the function `pystac.ItemCollection.from_file()` to load the search results as an `Item` list.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "727637ca-374f-47cf-884d-d2b2766742c1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pystac\n",
    "items = pystac.ItemCollection.from_file(\"data/search.json\")\n",
    "items"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94c70f55-7825-4083-af2d-c943b76f9cc2",
   "metadata": {},
   "source": [
    "In the search results, we have 7 `Item` type objects, corresponding to several Sentinel-2 scenes from March 21th and 28th in 2020. We will focus on the scene `S2A_31UFU_20200328_0_L2A`, and load band `nir09` (central wavelength 945 nm). We can load this band using the function `rioxarray.open_rasterio()`, via the Hypertext Reference `href` (commonly referred to as a URL):\n",
    "\n",
    "## **Exercise:** finding the right item and asset\n",
    "How do we go about selecting the correct item and asset from our ItemCollection we just loaded?\n",
    "1. Find the item corresponding to scene S2A_31UFU_20200328_0_L2A. Use the following loop:\n",
    "```\n",
    "for item in items:\n",
    "    if item.id == \"S2A_31UFU_20200328_0_L2A\":\n",
    "        break\n",
    "item\n",
    "```\n",
    "This loops through each item in the `items` collection we just retrieved. It stops (because of the `break` when the item.id matches the one we want which means `item` is then equal to the corresponding item. \n",
    "\n",
    "2. Find the asset `href` for the `nir09` band in the item's asset dictionary. (Hint: use the `item.assets[].href` syntax).\n",
    "3. Load it using rioxarray's `open_rasterio` method into a variable called `raster_ams_b9`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "720cdf4b-6472-483c-9e48-2e7e6817392f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Try something in here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91e2911c-fa5d-473c-8d72-3430607da2c0",
   "metadata": {},
   "source": [
    "By calling the variable name in the jupyter notebook we can get a quick look at the shape and attributes of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eedd0705-06ee-4490-9aa0-9677e495a390",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_ams_b9"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ba8bbc2-b368-4429-8bd9-32088e3b7e9a",
   "metadata": {},
   "source": [
    "The first call to `rioxarray.open_rasterio()` opens the file from remote or local storage, and then returns a `xarray.DataArray` object. The object is stored in a variable, i.e. `raster_ams_b9`. Reading in the data with `xarray` instead of `rioxarray` also returns a `xarray.DataArray`, but the output will not contain the geospatial metadata (such as projection information). You can use numpy functions or built-in Python math operators on a `xarray.DataArray` just like a numpy array. Calling the variable name of the `DataArray` also prints out all of its metadata information.\n",
    "\n",
    "The output tells us that we are looking at an `xarray.DataArray`, with `1` band, `1830` rows, and `1830` columns. We can also see the number of pixel values in the `DataArray`, and the type of those pixel values, which is unsigned integer (or `uint16`). The `DataArray` also stores different values for the coordinates of the `DataArray`. When using `rioxarray`, the term coordinates refers to spatial coordinates like `x` and `y` but also the `band` coordinate. Each of these sequences of values has its own data type, like `float64` for the spatial coordinates and `int64` for the `band` coordinate.\n",
    "\n",
    "This `DataArray` object also has a couple of attributes that are accessed like `.rio.crs`, `.rio.nodata`, and `.rio.bounds()`, which contain the metadata for the file we opened. Note that many of the metadata are accessed as attributes without `()`, but `bounds()` is a method (i.e. a function in an object) and needs parentheses.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28143c46-f7dd-464b-94b5-431cf7cc81d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(raster_ams_b9.rio.crs)\n",
    "print(raster_ams_b9.rio.nodata)\n",
    "print(raster_ams_b9.rio.bounds())\n",
    "print(raster_ams_b9.rio.width)\n",
    "print(raster_ams_b9.rio.height)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93f6f2eb-ccfd-4512-a2a1-4ca420946b75",
   "metadata": {},
   "source": [
    "The Coordinate Reference System, or `raster_ams_b9.rio.crs`, is reported as the string `EPSG:32631`. The `nodata` value is encoded as 0 and the bounding box corners of our raster are represented by the output of `.bounds()` as a `tuple` (like a list but you can't edit it). The height and width match what we saw when we printed the `DataArray`, but by using `.rio.width` and `.rio.height` we can access these values if we need them in calculations.\n",
    "\n",
    "We will be exploring this data throughout this episode. By the end of this episode, you will be able to understand and explain the metadata output.\n",
    "\n",
    "\n",
    "<span style=\"color: blue;\">*TIP: To improve code readability, file and object names should be used that make it clear what is in the file. The data for this episode covers Amsterdam, and is from Band 9, so we'll use a naming convention of `raster_ams_b9` for the variable name.*</span>\n",
    "\n",
    "\n",
    "## Visualize a Raster\n",
    "\n",
    "After viewing the attributes of our raster, we can examine the raw values of the array with `.values`:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a8ac587-b7ab-4e59-9262-702790888e43",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_b9.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bb2fe6f-ea90-43d8-ba98-709860e2caed",
   "metadata": {},
   "source": [
    "This can give us a quick view of the values of our array, but only at the corners. Since our raster is loaded in python as a `DataArray` type, we can plot this in one line similar to a pandas `DataFrame` with `DataArray.plot()`.\n",
    "\n",
    "__Exercise: plot our raster file using the plot() method__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "959a6f86-4db2-4d14-bd10-f870d9796160",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_ams_b9..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6fae531-5581-49be-b40f-7264ce5f1844",
   "metadata": {},
   "source": [
    "Notice that `rioxarray` helpfully allows us to plot this raster with spatial coordinates on the x and y axis (this is not the default in many cases with other functions or libraries).\n",
    "\n",
    "This plot shows the satellite measurement of the spectral band `nir09` for an area that covers part of the Netherlands. According to the [Sentinel-2 documentaion](https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/msi-instrument), this is a band with the central wavelength of 945nm, which is sensitive to water vapor. It has a spatial resolution of 60m. Note that the `band=1` in the image title refers to the ordering of all the bands in the  `DataArray`, not the Sentinel-2 band number `09` that we saw in the pystac search results.\n",
    "\n",
    "With a quick view of the image, we notice that half of the image is blank, no data is captured. We also see that the cloudy pixels at the top have high reflectance values, while the contrast of everything else is quite low. This is expected because this band is sensitive to the water vapor. However if one would like to have a better color contrast, one can add the option `robust=True`, which displays values between the 2nd and 98th percentile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2f2b9d-ed77-4325-8771-90cdcae636d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_b9.plot(robust=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad3e3a1b-651d-4400-9be3-f246178bbf8e",
   "metadata": {},
   "source": [
    "Now the color limit is set in a way fitting most of the values in the image. We have a better view of the ground pixels.\n",
    "\n",
    "---\n",
    "\n",
    "*NOTE: The option `robust=True` always forces displaying values between the 2nd and 98th percentile. Of course, this will not work for every case. For a customized displaying range, you can also manually specifying the keywords `vmin` and `vmax`. For example ploting between `100` and `7000`:*\n",
    "\n",
    "__Exercise: plot the raster with vmin and vmax arguments__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a93139e6-8380-4b36-a45d-e6d8d4da5531",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_b9 ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "211a5bc3-749a-467e-a787-79dac2562851",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## View Raster Coordinate Reference System (CRS) in Python\n",
    "Another information that we're interested in is the CRS, and it can be accessed with `.rio.crs`. To find out more about CRS look at [the earlier\n",
    "episode](https://carpentries-incubator.github.io/geospatial-python/instructor/03-crs.html) in the software carpentry course.\n",
    "Now we will see how features of the CRS appear in our data file and what\n",
    "meanings they have. We can view the CRS string associated with our DataArray's `rio` object using the `crs`\n",
    "attribute.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2024c3e-ceb2-4b7b-9420-25e727f2047e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(raster_ams_b9.rio.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a06dc157-4d79-4e58-bdea-27c63c3c0ee8",
   "metadata": {},
   "source": [
    "To print the EPSG code number (unique identifier used to represent coordinate systems) as an `int`, we use the `.to_epsg()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86c0e9eb-5587-4e00-8ec2-09eb09d3248d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_ams_b9.rio.crs.to_epsg()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d65ca6a-2d92-4aee-95ab-2e2db3f4762a",
   "metadata": {},
   "source": [
    "EPSG codes are great for succinctly representing a particular coordinate reference system. But what if we want to see more details about the CRS, like the units? For that, we can use `pyproj`, a library for representing and working with coordinate reference systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5803b0ea-3ae6-41dc-95a4-da0d3231d604",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyproj import CRS\n",
    "epsg = raster_ams_b9.rio.crs.to_epsg()\n",
    "crs = CRS(epsg)\n",
    "crs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30c00eb9-a206-43ea-a467-59b9ef161b2d",
   "metadata": {},
   "source": [
    "The `CRS` class from the `pyproj` library allows us to create a `CRS` object with methods and attributes for accessing specific information about a CRS, or the detailed summary shown above.\n",
    "\n",
    "A particularly useful attribute is `area_of_use`, which shows the geographic bounds that the CRS is intended to be used.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5caeb4d-cfde-4726-b265-f316cef2896a",
   "metadata": {},
   "outputs": [],
   "source": [
    "crs.area_of_use"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "396b9eac-2847-4073-9233-0746d578eb37",
   "metadata": {},
   "source": [
    "## **Exercise**: find the axes units of the CRS\n",
    "What units are our data in? See if you can find a method to examine this information using `help(crs)` or `dir(crs)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d1f760a-6fd9-496e-bfdc-531e7f6a1f95",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try something in here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f573adc5-b5de-4f8a-9441-d123722446d0",
   "metadata": {},
   "source": [
    "Let's break down the pieces of the `pyproj` CRS summary. The string contains all of the individual CRS elements that Python or another GIS might need, separated into distinct sections, and datum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "108a05cf-5faa-41fd-8451-03922edfe641",
   "metadata": {},
   "outputs": [],
   "source": [
    "crs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1adeed3-3871-4240-8e1f-468dc952ffce",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "* **Name** of the projection is UTM zone 31N (UTM has 60 zones, each 6-degrees of longitude in width). The underlying datum is WGS84. Universal Transverse Mercator (UTM) is a map projection system for assigning coordinates to locations.\n",
    "* **Axis Info**: the CRS shows a Cartesian system with two axes, easting and northing, in meter units.\n",
    "* **Area of Use**: the projection is used for a particular range of longitudes `0°E to 6°E` in the northern hemisphere (`0.0°N to 84.0°N`)\n",
    "* **Coordinate Operation**: the operation to project the coordinates (if it is projected) onto a cartesian (x, y) plane. Transverse Mercator is accurate for areas with longitudinal widths of a few degrees, hence the distinct UTM zones.\n",
    "* **Datum**: Details about the datum, or the reference point for coordinates. `WGS 84` and `NAD 1983` are common datums. `NAD 1983` is [set to be replaced in 2022](https://en.wikipedia.org/wiki/Datum_of_2022).\n",
    "\n",
    "Note that the zone is unique to the UTM projection. Not all CRSs will have a\n",
    "zone. Below is a simplified view of US UTM zones.\n",
    "\n",
    "![UTMZones](https://upload.wikimedia.org/wikipedia/commons/thumb/8/8d/Utm-zones-USA.svg/1920px-Utm-zones-USA.svg.png)\n",
    "###### The UTM zones across the continental United States (Chrismurf at English Wikipedia, via [Wikimedia Commons](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#/media/File:Utm-zones-USA.svg) (CC-BY))\n",
    "\n",
    "## Calculate Raster Statistics\n",
    "\n",
    "It is useful to know the minimum or maximum values of a raster dataset. __Exercise: compute these and other descriptive statistics with `min`, `max`, `mean`, and `std`.__\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c46bd232-2850-4cca-963d-5cf571d109d1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13dd8ed9-18a3-4466-b519-64546c66b7c9",
   "metadata": {},
   "source": [
    "The information above includes a report of the min, max, mean, and standard deviation values, along with the data type. If we want to see specific quantiles, we can use xarray's `.quantile()` method. For example for the 25% and 75% quantiles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bd4db6c-374b-48d1-a1c8-db4f7472fb33",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(raster_ams_b9.quantile([0.25, 0.75]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05138d80-9b26-4b92-8749-df2d69b1473b",
   "metadata": {},
   "source": [
    "---\n",
    "*NOTE: You could also get each of these values one by one using `numpy`.*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe927ff4-6467-4488-a671-05e59a884d12",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "print(numpy.percentile(raster_ams_b9, 25))\n",
    "print(numpy.percentile(raster_ams_b9, 75))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d64fa87d-3e80-4fb7-b973-a6595d37b8f9",
   "metadata": {},
   "source": [
    "You may notice that `raster_ams_b9.quantile` and `numpy.percentile` didn't require an argument specifying the axis or dimension along which to compute the quantile. This is because `axis=None` is the default for most numpy functions, and therefore `dim=None` is the default for most xarray methods. It's always good to check out the docs on a function to see what the default arguments are, particularly when working with multi-dimensional image data. To do so, we can use`help(raster_ams_b9.quantile)` (or `?raster_ams_b9.percentile` in jupyter notebook), e.g.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25334bd8-d927-4df3-987a-b2f0fc2fd443",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "?raster_ams_b9.quantile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f8a771b-b8c4-4096-b887-f23bf5fb7b15",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Dealing with Missing Data\n",
    "So far, we have visualized a band of a Sentinel-2 scene and calculated its statistics. However, we need to take missing data into account. Raster data often has a \"no data value\" associated with it and for raster datasets read in by `rioxarray`. This value is referred to as `nodata`. This is a value assigned to pixels where data is missing or no data were collected. There can be different cases that cause missing data, and it's common for other values in a raster to represent different cases. The most common example is missing data at the edges of rasters.\n",
    "\n",
    "By default the shape of a raster is always rectangular. So if we have a dataset that has a shape that isn't rectangular, some pixels at the edge of the raster will have no data values. This often happens when the data were collected by a sensor which only flew over some part of a defined region.\n",
    "\n",
    "As we have seen above, the `nodata` value of this dataset (`raster_ams_b9.rio.nodata`) is 0. When we have plotted the band data, or calculated statistics, the missing value was not distinguished from other values. Missing data may cause some unexpected results. For example, the 25th percentile we just calculated was 0, probably reflecting the presence of a lot of missing data in the raster.\n",
    "\n",
    "To distinguish missing data from real data, one possible way is to use `NaN` to represent them. This can be done by specifying `masked=True` when loading the raster:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aefc091-a57f-48c3-aba0-aa415a516495",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "raster_ams_b9 = rioxarray.open_rasterio(items[0].assets[\"nir09\"].href, masked=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1e9160d-db4c-4c49-89c8-1d5181bfe723",
   "metadata": {},
   "source": [
    "One can also use the `where` function to select all the pixels which are different from the `nodata` value of the raster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9704ddf0-8692-496d-82f5-d0dd53f0b069",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_b9.where(raster_ams_b9!=raster_ams_b9.rio.nodata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caf54238-6a6c-49f4-9813-918c33a13f93",
   "metadata": {},
   "source": [
    "Either way will change the `nodata` value from 0 to `nan`. Now if we compute the statistics again the missing data will not be considered.\n",
    "\n",
    "__Exercise: Compute the statistics (`min`, `max`, `mean`, `std`) again:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08403853-b38c-49e8-b5f1-06ff6d885cdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)\n",
    "print(raster_ams_b9...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97bd2c01-819b-4e52-8e70-f0ec67a1f8a8",
   "metadata": {},
   "source": [
    "And if we plot the image, the `nodata` pixels are not shown because they are not 0 anymore. \n",
    "\n",
    "__Exercise: plot the masked image with `robust` set to true__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a827b027-9b97-4c4b-808d-0e87a9fa6b18",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_b9..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3fd6acc-b9cf-4a87-abcb-86d0e1af6d39",
   "metadata": {},
   "source": [
    "One should notice that there is a side effect of using `nan` instead of `0` to represent the missing data: the data type of the `DataArray` was changed from integers to float. This need to be taken into consideration when the data type matters in your application.\n",
    "\n",
    "## Raster Bands\n",
    "So far we looked into a single band raster, i.e. the `nir09` band of a Sentinel-2 scene. However, to get a smaller, non georeferenced version of the scene, one may also want to visualize the true-color overview of the region. This is provided as a multi-band raster -- a raster dataset that contains more than one band.\n",
    "\n",
    "![Sketch of a multi-band raster image](https://carpentries-incubator.github.io/geospatial-python/fig/E06/single_multi_raster.png)\n",
    "###### Sketch of a multi-band raster image\n",
    "\n",
    "The `overview` asset in the Sentinel-2 scene is a multiband asset. Similar to `nir09`, we can load it by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a63e401b-e438-4291-94ab-c45aa378732a",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_overview = rioxarray.open_rasterio(items[0].assets['visual'].href, overview_level=3)\n",
    "raster_ams_overview\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ac1d52-a434-46a2-bbbb-5869e554a851",
   "metadata": {},
   "source": [
    "The band number comes first when GeoTiffs are read with the `.open_rasterio()` function. As we can see in the `xarray.DataArray` object, the shape is now `(band: 3, y: 687, x: 687)`, with three bands in the `band` dimension. It's always a good idea to examine the shape of the raster array you are working with and make sure it's what you expect. Many functions, especially the ones that plot images, expect a raster array to have a particular shape. One can also check the shape using the `.shape` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13291ec-3701-4bbd-96cf-b33c8b76f165",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_overview.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5097c1d-f1ef-45da-9ac6-c2b84e408c99",
   "metadata": {},
   "source": [
    "One can visualize the multi-band data with the `DataArray.plot.imshow()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91c685c2-4375-4822-9829-bc6ae6a01cf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_ams_overview.plot.imshow()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37ea12ba-fd14-4f5d-bf65-915334df90a9",
   "metadata": {},
   "source": [
    "Note that the `DataArray.plot.imshow()` function makes assumptions about the shape of the input DataArray, that since it has three channels, the correct colormap for these channels is RGB. It does not work directly on image arrays with more than 3 channels. One can replace one of the RGB channels with another band, to make a false-color image.\n",
    "\n",
    "## **Exercise**: set the plotting aspect ratio\n",
    "As seen in the figure above, the true-color image is stretched. Visualize it with the right aspect ratio. You can use the [documentation](https://xarray.pydata.org/en/stable/generated/xarray.DataArray.plot.imshow.html) of `DataArray.plot.imshow()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26d2ac75-1bd5-4594-8b97-538742660295",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try something in here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b15a23d3-8fd3-40e5-9846-a5ba9f4c7175",
   "metadata": {},
   "source": [
    "## Key takeaways:\n",
    "- `rioxarray` and `xarray` are for working with multidimensional arrays like pandas is for working with tabular data.\n",
    "- `rioxarray` stores CRS information as a CRS object that can be converted to an EPSG code or PROJ4 string.\n",
    "- Missing raster data are filled with nodata values, which should be handled with care for statistics and visualization."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 + Jaspy",
   "language": "python",
   "name": "jaspy"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}