{ "cells": [ { "cell_type": "markdown", "id": "7d794375-5d42-42b3-8780-f4259ed87cf3", "metadata": {}, "source": [ "# Exploring PACE OCI L2 Surface Reflectance " ] }, { "cell_type": "markdown", "id": "41157cb0-b9ce-4123-b1ec-4c9ffb3fa912", "metadata": {}, "source": [ "In this tutorial we will work with Level 2 data from the PACE's Ocean Color Instrument (OCI), specificaclly the Surface Reflectance data product. For information on how this data product has been created please visit: [PACE OCI Surface Reflectance (SFREFL) Algorithm Theortical Basis Document (ATBD)]\n", "\n", "The OCI instrument is a hyperspectral imaging radiometer, collecting hyperspectral measurements from 340 nm - 895 nm (Ultra Violet (UV) to near-infrared-NIR), and multi-spectral measurements from 940-2260 nm (NIR to shortwave-infrared). \n", "\n", "This tutorial leverages code from [Ocean Color OB DAAC Training Tutorials] and [NASA LP DAAC VITALS tutorials] \n", "\n", "[Ocean Color OB DAAC Training Tutorials]: https://oceancolor.gsfc.nasa.gov/resources/docs/tutorials/ \n", "[NASA LP DAAC VITALS tutorials]: nasa.github.io/VITALS/\n", "[PACE OCI Surface Reflectance (SFREFL) Algorithm Theortical Basis Document (ATBD)]: https://www.earthdata.nasa.gov/apt/documents/sfrefl/v1.0#doc-header\n" ] }, { "cell_type": "markdown", "id": "24e5d0e8-3e08-4fec-b008-b769325ebfeb", "metadata": {}, "source": [ "### Learning Objectives\n", "\n", "This tutorial will guide you through working with PACE OCI L2 SFREL data on a cloud platform. Upon completion of this tutorial, you will be able to:\n", "\n", "- Access and search for available granules by time and location.\n", "- Explore data structure and format.\n", "- Visualize data by plotting individual spectral bands.\n", "- Apply cloud and water flags to effectively mask out unwanted pixels.\n", "- Create interactive RGB maps.\n", "- Display, and export spectral signatures for selected locations.\n", "- Calculate advanced hyperspectral indices using these L2 datasets.\n", " " ] }, { "cell_type": "markdown", "id": "feed5245-6397-46f1-8753-acd710280c39", "metadata": {}, "source": [ "### Contents\n", "* [Setup](#setup)\n", "* [Searching for collections and granules](#searching-for-collections)\n", "* [Read PACE L2 Surface Reflectance data](#reading-dataset) \n", "* [Mask clouds and water from images](#masking-clouds)\n", "* [Display spectral signatures for selected locations](#interactive-spatial-and-spectral-plots)\n", "* [Calculate unique spectral indices](#calculate-unique-indices)" ] }, { "cell_type": "markdown", "id": "990e9b50-5862-4a54-a1be-62d04b10451a", "metadata": {}, "source": [ "#### Setup\n" ] }, { "cell_type": "code", "execution_count": null, "id": "bb5daa5b-3525-4fb0-8b95-b9174057b4b6", "metadata": {}, "outputs": [], "source": [ "!conda install cf_xarray -y" ] }, { "cell_type": "code", "execution_count": null, "id": "dcf11f93-1bbf-4be3-ae1c-af2fdf885ae6", "metadata": {}, "outputs": [], "source": [ "# Import required libraries\n", "import csv\n", "import math\n", "from io import BytesIO\n", "\n", "import cf_xarray\n", "import earthaccess\n", "import rasterio as rio\n", "import xarray as xr\n", "import cartopy\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import matplotlib.pyplot as plt\n", "import holoviews as hv\n", "import geoviews as gv\n", "import hvplot.xarray\n" ] }, { "cell_type": "markdown", "id": "3279bed2", "metadata": {}, "source": [ "## Authentication\n", "\n", "`earthaccess` creates and leverages Earthdata Login tokens to authenticate with NASA systems. Earthdata Login tokens expire after a month. To retrieve a token from Earthdata Login, you can either enter your username and password each time you use `earthaccess`, or use a `.netrc` file. A `.netrc` file is a configuration file that is commonly used to store login credentials for remote systems. If you don't have a `.netrc` or don't know if you have one or not, you can use the `persist` argument with the `login` function below to create or update an existing one, then use it for authentication.\n", "\n", "If you do not have an Earthdata Account, you can create one [here](https://urs.earthdata.nasa.gov/home)." ] }, { "cell_type": "code", "execution_count": null, "id": "e74977ca-726f-4809-8360-fe4baa5f5f86", "metadata": {}, "outputs": [], "source": [ "auth = earthaccess.login(persist=True)" ] }, { "cell_type": "markdown", "id": "8cb1f61f", "metadata": {}, "source": [ "## Searching for Collections\n", "\n", "The PACE mission produces several collections or datasets available via the NASA Earthdata cloud archive. \n", "\n", "To view what's available, we can use the `search_datasets` function and with the `keyword` argument.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0f3ad80a", "metadata": {}, "outputs": [], "source": [ "# Retrieve Collections\n", "collections = earthaccess.search_datasets(keyword='PACE Surface Reflectance')\n", "# Print Quantity of Results\n", "print(f'Collections found: {len(collections)}')" ] }, { "cell_type": "markdown", "id": "d7f0579f", "metadata": {}, "source": [ "If you print the `collections` object you can explore all of the json metadata." ] }, { "cell_type": "code", "execution_count": null, "id": "30c60d4b", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# # Print collections\n", "# collections" ] }, { "cell_type": "markdown", "id": "fd3be9e1", "metadata": {}, "source": [ "We can also create a list of the `short-name`, `concept-id`, and `version` of each result collection using list comprehension. These fields are important for specifying and searching for data within collections. " ] }, { "cell_type": "code", "execution_count": null, "id": "1f87b5ef", "metadata": {}, "outputs": [], "source": [ "collections_info = [\n", " {\n", " 'short_name': c['umm']['ShortName'],\n", " 'collection_concept_id': c['meta']['concept-id'],\n", " 'version': c['umm']['Version'],\n", " 'entry_title': c['umm']['EntryTitle']\n", " }\n", " for c in collections\n", "]\n", "pd.set_option('display.max_colwidth', 150)\n", "collections_info = pd.DataFrame(collections_info)\n", "collections_info" ] }, { "cell_type": "markdown", "id": "dc4830bc", "metadata": {}, "source": [ "The collection `concept-id` is the best way to search for data within a collection, as this is unique to each collection. The `short-name` can be used as well, however the `version` should be passed as well as there can be multiple versions available with the same short name. After finding the collection you want to search, you can use the `concept-id` to search for granules within that collection." ] }, { "cell_type": "markdown", "id": "9fbdf386", "metadata": {}, "source": [ "## Searching for Granules\n", "\n", "A `granule` can be thought of as a unique spatiotemporal grouping within a collection. To search for `granules`, we can use the `search_data` function from `earthaccess` and provide the arguments for our search. Its possible to specify search products using several criteria shown in the table below:\n", "\n", "|dataset origin and location|spatio temporal parameters|dataset metadata parameters|\n", "|:---|:---|:---|\n", "|archive_center|bounding_box|concept_id\n", "|data_center|temporal|entry_title\n", "|daac|point|granule_name\n", "|provider|polygon|version\n", "|cloud_hosted|line|short_name\n", "\n", "For this example we will search for the **PACE OCI Level-2 Regional Surface Reflectance Data** using a bounding box and temporal parameters, and add a `cloud_cover` parameter. Note that not all datasets have cloud cover information, so this parameter may not work for all datasets." ] }, { "cell_type": "code", "execution_count": null, "id": "093dfcd2-231c-4340-b703-f985baddb641", "metadata": {}, "outputs": [], "source": [ "tspan = (\"2025-03-01\", \"2025-03-20\")\n", "bbox = (-86.512, 31.70, -86.512, 34.938)\n", " \n", "results = earthaccess.search_data(\n", " short_name=\"PACE_OCI_L2_SFREFL\",\n", " version=\"3.1\",\n", " temporal=tspan,\n", " bounding_box=bbox,\n", " cloud_cover=(0, 50)\n", ")\n", "print(f'Granules found: {len(results)}')" ] }, { "cell_type": "markdown", "id": "36b71f40", "metadata": {}, "source": [ "Results is a `list`, so we can use an index to view a single result." ] }, { "cell_type": "code", "execution_count": null, "id": "2b912222", "metadata": {}, "outputs": [], "source": [ "results[0]" ] }, { "cell_type": "markdown", "id": "76409ff2", "metadata": {}, "source": [ "We can also retrieve specific metadata for a result using `.keys()` since this object also acts as a dictionary." ] }, { "cell_type": "code", "execution_count": null, "id": "47f8b03d", "metadata": {}, "outputs": [], "source": [ "results[0].keys()\n", "# results[0]['meta']" ] }, { "cell_type": "markdown", "id": "d7d93e8c", "metadata": {}, "source": [ "Feel free to explore the metadata of the results, to check if there are fields relevant for your use case.\n", "\n", "From here, we can also retrieve the links for data using `earthaccess`, which we can either download or stream." ] }, { "cell_type": "code", "execution_count": null, "id": "117f1b47", "metadata": {}, "outputs": [], "source": [ "links = [result.data_links()[0] for result in results]\n", "# Display the first 5 links\n", "links[:5]" ] }, { "cell_type": "markdown", "id": "c43bc921", "metadata": {}, "source": [ "To access data from NASA, you'll need to provide your Earthdata Login credentials. When streaming this can best be done using the token or cookies set up by the `earthaccess` library. Since we've already logged in, we can start an `fsspec` session to manage our connection to a remote file, including sending credentials. This allows other libraries to work with a URL as if it is a local file.\n", "\n", "In addition to `fsspec` for this example, we'll also use `BytesIO` to read the full dataset into memory. This is not always the best approach, but for working with a single PACE scene which is ~750 mb, it speeds up the process due to the data structure and some caveats regarding `fsspec`." ] }, { "cell_type": "code", "execution_count": null, "id": "3c049707-5abc-4b2f-8881-cac27ec99bbd", "metadata": {}, "outputs": [], "source": [ "fs = earthaccess.get_fsspec_https_session()" ] }, { "cell_type": "markdown", "id": "032069b8", "metadata": {}, "source": [ "Now, open the connection to the remote file and read into memory." ] }, { "cell_type": "code", "execution_count": null, "id": "f25442a1", "metadata": {}, "outputs": [], "source": [ "with fs.open(links[0]) as f:\n", " file_content = BytesIO(f.read())" ] }, { "cell_type": "markdown", "id": "251c22dd-72f2-4790-97af-4c30e79d3950", "metadata": {}, "source": [ "## Reading dataset" ] }, { "cell_type": "markdown", "id": "99c18aec", "metadata": {}, "source": [ "The PACE_OCI_L2_SFREFL files are heirarchichal netCDF files, so we will use the `xarray.datatree` function to read in and view all of the groups. Using the `datatree` representation object, we can explore the data structure, which is useful for understanding the groups, variables, dimension order, and coordinates. " ] }, { "cell_type": "code", "execution_count": null, "id": "f37c48ba-6c62-4c7b-a5c3-78fd4b7fd681", "metadata": { "scrolled": true }, "outputs": [], "source": [ "datatree = xr.open_datatree(file_content, decode_timedelta=False)\n", "datatree" ] }, { "cell_type": "markdown", "id": "d809f756", "metadata": {}, "source": [ "It will be easier to work with the data if we create a single dataset, which we can do by merging the data and coordinate variables we need from different groups. Longitude and latitude appear as data variables, in the group `navigation_data`, they need to be explicitly set as coordinates.\n", "Note the group `geophysical_data` and its data variables: `rhos` and `l2_flags`. The `rhos` variable are surface reflectances, and the `l2_flags` are quality flags as defined by the Ocean Biology Processing Group (OBPG)." ] }, { "cell_type": "code", "execution_count": null, "id": "cdc767f7", "metadata": { "scrolled": true }, "outputs": [], "source": [ "ds = xr.merge(\n", " (\n", " datatree.ds,\n", " datatree[\"geophysical_data\"].ds[[\"rhos\", \"l2_flags\"]],\n", " datatree[\"sensor_band_parameters\"].coords,\n", " datatree[\"navigation_data\"].ds.set_coords((\"longitude\", \"latitude\")).coords,\n", " )\n", ")\n", "ds" ] }, { "cell_type": "markdown", "id": "8efd2a24-b27b-441a-be72-35d62349c3d4", "metadata": {}, "source": [ "We can also see which wavelengths we have surface reflectance measurements at by accessing the `wavelength_3d` coordinate:" ] }, { "cell_type": "code", "execution_count": null, "id": "4833d7a3-1e3a-431a-9c03-9d43cbfd149b", "metadata": {}, "outputs": [], "source": [ "# Check the wavelengths available for PACE \n", "ds[\"wavelength_3d\"]" ] }, { "cell_type": "markdown", "id": "9925a0b6-7f9f-4530-bd79-91f565d47358", "metadata": {}, "source": [ "Note that “wavelength_3d” is an indexed coordinate, which allows us to subset the dataset by slicing or choosing individual wavelengths. The method=\"nearest\" argument lets us select one wavelength without knowning its exact value.\n", "\n", "Hence, we can select one wavelength even if it does not represent an exact value in the “wavelength_3d” array and plot it using the method=\"nearest\" argument. " ] }, { "cell_type": "code", "execution_count": null, "id": "bca93c87-e4f6-4c29-a54f-b0dd7d5a3d99", "metadata": {}, "outputs": [], "source": [ "rhos_860 = ds[\"rhos\"].sel({\"wavelength_3d\": 860}, method=\"nearest\")\n", "\n", "fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={\"projection\": cartopy.crs.PlateCarree()})\n", "ax.gridlines(draw_labels={\"left\": \"y\", \"bottom\": \"x\"}, linewidth=0.25)\n", "ax.coastlines(linewidth=0.5)\n", "ax.add_feature(cartopy.feature.OCEAN, edgecolor=\"w\", linewidth=0.01)\n", "ax.add_feature(cartopy.feature.LAND, edgecolor=\"w\", linewidth=0.01)\n", "rhos_860.plot(x=\"longitude\", y=\"latitude\", cmap=\"Greys_r\", vmin=0, vmax=1.0)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "a1397c6b-6454-438d-9308-f15ea421d0f6", "metadata": { "scrolled": true }, "outputs": [], "source": [ "rhos_860" ] }, { "cell_type": "markdown", "id": "fb7efe77-943e-4bf1-9da1-b2d15df33fbd", "metadata": {}, "source": [ "## Masking clouds" ] }, { "cell_type": "code", "execution_count": null, "id": "34dc9615-835d-4fc4-a91c-05da71b3334f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "#Now flags \n", "ds[\"l2_flags\"]" ] }, { "cell_type": "markdown", "id": "562a09d7", "metadata": {}, "source": [ "We will use the `cf_xarray` library to create a mask including only land and excluding clouds and ice." ] }, { "cell_type": "code", "execution_count": null, "id": "65c17194-9500-40e1-a6c2-f89c694dc098", "metadata": {}, "outputs": [], "source": [ "# Create and apply Mask\n", "ds[\"l2_flags\"].cf.is_flag_variable\n", "cldwater_mask = (\n", " (ds[\"l2_flags\"].cf == \"LAND\")\n", "& ~(ds[\"l2_flags\"].cf == \"CLDICE\")\n", ")\n", "# Apply mask, creating new dataset\n", "rhos = ds[\"rhos\"].where(cldwater_mask)" ] }, { "cell_type": "code", "execution_count": null, "id": "27e034e2-7a59-4499-9cc0-874d7c44feea", "metadata": {}, "outputs": [], "source": [ "# Plot a band of masked data\n", "rhos_860 = rhos.sel({\"wavelength_3d\": 860}, method=\"nearest\")\n", "\n", "fig, ax = plt.subplots(figsize=(9, 5), subplot_kw={\"projection\": cartopy.crs.PlateCarree()})\n", "ax.coastlines(linewidth=0.25)\n", "ax.gridlines(draw_labels={\"left\": \"y\", \"bottom\": \"x\"}, linewidth=0.25)\n", "ax.add_feature(cartopy.feature.OCEAN, edgecolor=\"w\", linewidth=0.01)\n", "ax.add_feature(cartopy.feature.LAND, edgecolor=\"w\", linewidth=0.01)\n", "ax.add_feature(cartopy.feature.LAKES, edgecolor=\"k\", linewidth=0.1)\n", "rhos_860.plot(x=\"longitude\", y=\"latitude\", cmap=\"Greys_r\", vmin=0, vmax=1.0)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5e92a63f-9c36-4be5-bd88-b56b3cdff6aa", "metadata": {}, "source": [ "## Projecting data " ] }, { "cell_type": "markdown", "id": "d1cda943", "metadata": {}, "source": [ "From viewing the data structure, we can see that the data is in swath format based on the dimensions. This means that we will need to grid the data to work with it in a typical manner. \n", "\n", "There are a couple ways to do this transformation, but for this example, we'll opt for a simple approach using `rioxarray`. Note that the process shown here will create a grid unique to this scene, so it won't easily merge or stack with others created in the same way for something like a time-series analysis. \n", "\n", "We encourage you to check the following resources for additional information on data transformation for PACE datasets: [SeaDAS] and [Python Jupyter notebook] \n", "\n", "To grid the data with `rioxarray`, we need to change the order of the dimensions in the array. Hyperspectral data is often stored (y,x,band) since most processing is conducted along the band dimension; however, `rioxarray` expects the dimensions to be in the order of (band, y, x), so we must transpose the array to grid the data using `rioxarray`. Note that we'll also want to ensure we're using nearest neighbor resampling to preserve the spectral information during the gridding process. If we use another method we could be creating artificial spectra which could lead to inaccurate unmixing or classification results.\n", "\n", "[SeaDAS]: https://seadas.gsfc.nasa.gov/client_server/\n", "[Python Jupyter notebook]: https://github.com/skyecaplan/pace_tutorials/blob/main/oci_project_and_format.ipynb" ] }, { "cell_type": "code", "execution_count": null, "id": "4ee1d846-c1d5-48af-a05a-aa83cda20a5f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Transpose data\n", "sr_src = rhos.transpose(\"wavelength_3d\", ...)\n", "sr_src" ] }, { "cell_type": "code", "execution_count": null, "id": "94108984-122e-4b54-8b87-d50bc97d643a", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Set spatial dimensions and swath CRS\n", "sr_src = sr_src.rio.set_spatial_dims(\"pixels_per_line\", \"number_of_lines\").rio.write_crs(\"epsg:4326\")\n", "sr_src" ] }, { "cell_type": "code", "execution_count": null, "id": "e34ea4f9-b1b0-4a63-8124-1f18b4ed11f4", "metadata": {}, "outputs": [], "source": [ "# Grid the data using rioxarray and rename dimensions\n", "sr_dst = sr_src.rio.reproject(\n", " dst_crs=sr_src.rio.crs,\n", " src_geoloc_array=(\n", " sr_src.coords[\"longitude\"],\n", " sr_src.coords[\"latitude\"],\n", " ),\n", " nodata=np.nan,\n", " resampling=rio.warp.Resampling.nearest,\n", ").rename({'y': 'latitude', 'x': 'longitude'})" ] }, { "cell_type": "code", "execution_count": null, "id": "165db082-a320-406d-b0ff-85c13ec67ff9", "metadata": { "scrolled": true }, "outputs": [], "source": [ "sr_dst" ] }, { "cell_type": "markdown", "id": "42bb2161", "metadata": {}, "source": [ "From this point, we can either subset the data spatially to a smaller region of interest, or extract a series of location points to work with.\n", "\n", "## Interactive Spatial and Spectral Plots\n", "\n", "Combining the Spatial and Spectral information into a single visualization can be a powerful tool for exploring and inspecting imaging spectroscopy data. Using the `streams` module from `Holoviews` we can link a spatial map to a plot of spectra.\n", "\n", "We could plot a single band image as we previously have, but using a multiband image, like an RGB may help infer what targets we're examining. Build function to select and adjust the brightness of bands to create a nice RGB image.\n", "\n", "Select bands to represent red (650 nm), green (560 nm), and blue (470 nm) by finding the nearest to a wavelength chosen to represent that color." ] }, { "cell_type": "code", "execution_count": null, "id": "7f09c1bd-d230-45d3-99d5-b446076ee71e", "metadata": { "scrolled": true }, "outputs": [], "source": [ "rgb_da = sr_dst.sel(wavelength_3d=[650, 560, 470], method='nearest')\n", "rgb_da" ] }, { "cell_type": "markdown", "id": "5971cd8e", "metadata": {}, "source": [ " To get a better rgb image, we can balance the brightness of the selected wavelengths to make a prettier map. **This will not affect the data, just the visuals.** To do this we will use the function below. We can change the bright argument to increase or decrease the brightness of the scene as a whole." ] }, { "cell_type": "code", "execution_count": null, "id": "719bf726", "metadata": {}, "outputs": [], "source": [ "def gamma_adjust(rgb_da, bright=0.3):\n", " \"\"\"\n", " Adjust gamma across all bands in the RGB dataset.\n", " \"\"\"\n", " array = rgb_da.data\n", " # Mask nan and negative values\n", " invalid = np.isnan(array) | (array < 0)\n", " valid = ~invalid\n", "\n", " # Calculate gamma based on the mean of valid values\n", " mean_valid = np.nanmean(array[valid])\n", " gamma = math.log(bright) / math.log(mean_valid)\n", " # Apply scaling and clip\n", " scaled = np.full_like(array, np.nan)\n", " scaled[valid] = np.power(array[valid], gamma)\n", " rgb_da.data = np.clip(scaled, 0, 1)\n", " return rgb_da" ] }, { "cell_type": "code", "execution_count": null, "id": "fa11a9b1", "metadata": {}, "outputs": [], "source": [ "rgb_da = gamma_adjust(rgb_da, bright=0.3)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9bf1d936-caa0-42be-9908-f72785424825", "metadata": {}, "outputs": [], "source": [ "rgb_da.hvplot.rgb(x='longitude', y='latitude', bands='wavelength_3d', frame_height=400, geo=True, crs='EPSG:4326', tiles='OSM', title=\"Stretched RGB PACE Image\")" ] }, { "cell_type": "markdown", "id": "8b0b8729", "metadata": {}, "source": [ "Now that we have an RGB dataset, we can use that to create a spatial plot, and data selected by clicking on our rgb map, which can be inputs for a function to return values from the full dataset at the selected latitude and longitude location using the cell below. To visualize the spectral and spatial data side-by-side, we use the `Point Draw` tool from the `holoviews` library.\n", "\n", "Click in the RGB image to add spectra to the plot. You can also click and hold the mouse button then drag previously placed points. To remove a point click and hold the mouse button down, then press the backspace key." ] }, { "cell_type": "code", "execution_count": null, "id": "47e7f1f7", "metadata": {}, "outputs": [], "source": [ "# Starting Point from center pixel\n", "x_start = -85.3879\n", "y_start = 31.3446\n", "\n", "# Set Point Limit\n", "POINT_LIMIT = 10\n", "\n", "# Set up Color Cycling\n", "color_cycle = hv.Cycle('Category20')\n", "\n", "first_point = ([x_start], [y_start], [0])\n", "\n", "points = gv.Points(first_point, kdims=['longitude','latitude'], vdims='id', crs=cartopy.crs.PlateCarree())\n", "\n", "points_stream = hv.streams.PointDraw(\n", " data=points.columns(),\n", " source=points,\n", " drag=True,\n", " num_objects=POINT_LIMIT,\n", " styles={'fill_color': color_cycle.values[:POINT_LIMIT], 'line_color': 'gray'}\n", ")\n", "\n", "# RGB Plot without Basemap\n", "\n", "rgb_map = rgb_da.hvplot.rgb(x='longitude', y='latitude', bands='wavelength_3d',\n", " frame_height=480, frame_width=480,\n", " crs=cartopy.crs.PlateCarree(),\n", " title=\"Stretched RGB PACE Image\")\n", "\n", "# Coastlines\n", "coastlines = gv.feature.coastline(projection=cartopy.crs.PlateCarree()).opts(line_color='black', line_width=1)\n", "\n", "posxy = hv.streams.PointerXY(source=rgb_map, x=x_start, y=y_start)\n", "clickxy = hv.streams.Tap(source=rgb_map, x=x_start, y=y_start)\n", "\n", "# Function to build spectral plot of clicked location to show on hover stream plot\n", "def click_spectra(data):\n", " coordinates = [c for c in zip(data['longitude'], data['latitude'])]\n", " \n", " plots = {}\n", " for i, coords in enumerate(coordinates):\n", " x, y = coords\n", " selected = sr_dst.sel(longitude=x, latitude=y, method=\"nearest\")\n", " plots[i] = (\n", " selected.hvplot.line(\n", " y=\"rhos\",\n", " x=\"wavelength_3d\",\n", " xlim=(340, 895), #you can comment or change the range here to see the multispectral bands too\n", " label=f\"{i}\"\n", " )\n", " )\n", " points_stream.data[\"id\"][i] = i\n", " return hv.NdOverlay(plots).opts(hv.opts.Curve(color=color_cycle))\n", "\n", "def hover_spectra(x,y):\n", " return sr_dst.sel(longitude=x,latitude=y,method='nearest').hvplot.line(y='rhos',x='wavelength_3d',\n", " color='black', frame_width=480)\n", "# Define the Dynamic Maps\n", "click_dmap = hv.DynamicMap(click_spectra, streams=[points_stream])\n", "hover_dmap = hv.DynamicMap(hover_spectra, streams=[posxy])\n", "\n", "# Plot the Map and Dynamic Map side by side\n", "hv.Layout(hover_dmap*click_dmap + rgb_map * coastlines* points).cols(2).opts(\n", " hv.opts.Points(active_tools=['point_draw'], size=10, tools=['hover'], color='white', line_color='gray'),\n", " hv.opts.Overlay(show_legend=False, show_title=False, fontscale=1.5, frame_height=480)\n", ")" ] }, { "cell_type": "markdown", "id": "74d91b1d", "metadata": {}, "source": [ "We can take these selected points and the corresponding reflectance spectra and save them as a .csv for later use." ] }, { "cell_type": "code", "execution_count": null, "id": "f5489ea3", "metadata": {}, "outputs": [], "source": [ "data = points_stream.data\n", "wavelengths = sr_dst.wavelength_3d.values\n", "\n", "rows = [[\"id\", \"longitude\", \"latitude\"] + [str(i) for i in wavelengths]]\n", " \n", "for p in zip(data['longitude'], data['latitude'], data['id']):\n", " x, y, i = p\n", " spectra = sr_dst.sel(longitude=x, latitude=y, method=\"nearest\").values\n", " row = [i, x, y] + list(spectra)\n", " rows.append(row)" ] }, { "cell_type": "code", "execution_count": null, "id": "abbee4f7", "metadata": {}, "outputs": [], "source": [ "with open('../data/pace_interactive_plot_data.csv', 'w', newline='') as f:\n", " writer = csv.writer(f)\n", " writer.writerows(rows)" ] }, { "cell_type": "code", "execution_count": null, "id": "b5983c82-b903-440c-9afd-3dee59e5774a", "metadata": {}, "outputs": [], "source": [ "# In case you cannot save the file, check your directory\n", "# import os\n", "# print(os.getcwd())" ] }, { "cell_type": "markdown", "id": "f12772ba-19d7-4f00-839a-179fc890c4a7", "metadata": {}, "source": [ "## Calculate unique indices \n", "### Hyperspectral-enabled Vegetation Indices (VI)\n", "We’ll take the Chlorophyll Index Red Edge (CIRE) as an example of a hyperspectral-enabled VI. CIRE uses bands from the red edge and the NIR to get at relative canopy chlorophyll content.\n", "\n", "CIRE: (ρ 800 /ρ 705)−1\n", "\n", "\n", "Carotenoid Content Index (Car):\n", "[(1/ρ495)−(1/ρ705)] * ρ800\n", "\n", "\n", "Information on Hyperspectral enabled indices by PACE [here]. \n", "\n", "[here]: https://www.earthdata.nasa.gov/apt/documents/landvi/v1.0#mathematical_theory " ] }, { "cell_type": "code", "execution_count": null, "id": "db74d5f6-7a36-47b3-856a-3e4e80325c4c", "metadata": {}, "outputs": [], "source": [ "## Chlorophyll Index Red Edge (CIRE)\n", "# Select bands\n", "sr_800 = sr_dst.sel({\"wavelength_3d\": 800}, method=\"nearest\")\n", "sr_705 = sr_dst.sel({\"wavelength_3d\": 705}, method=\"nearest\")\n", "#Calculate\n", "cire = (sr_800 / sr_705) - 1\n", "cire.attrs[\"long_name\"] = \"Chlorophyll Index Red Edge (CIRE)\"" ] }, { "cell_type": "code", "execution_count": null, "id": "0d0622de", "metadata": { "scrolled": true }, "outputs": [], "source": [ "cire" ] }, { "cell_type": "code", "execution_count": null, "id": "78b7bdaa-4c86-410d-a7f6-29b8023b60f9", "metadata": {}, "outputs": [], "source": [ "#Plot\n", "map1 = cire.hvplot.image(x='longitude', y='latitude', frame_height=480, geo=True, cmap=\"viridis\", tiles='CartoDark', title=\"CIRE from PACE OCI\")\n", "map1" ] }, { "cell_type": "code", "execution_count": null, "id": "c7b2a6e5-c161-4a07-90c9-aa27f1774871", "metadata": {}, "outputs": [], "source": [ "## Carotenoid Content Index (Car) \n", "# Select additional bands\n", "sr_495 = sr_dst.sel({\"wavelength_3d\": 495}, method=\"nearest\")\n", "\n", "#Calculate\n", "car = ((1 / sr_495)- (1 / sr_705)) * sr_800\n", "car.attrs[\"long_name\"] = \"Caretenoid Content Index (Car)\"" ] }, { "cell_type": "code", "execution_count": null, "id": "78257966-3ca3-4704-87ad-80a7ec6a89e3", "metadata": {}, "outputs": [], "source": [ "map2 = car.hvplot.image(x='longitude', y='latitude', frame_height=480, geo=True, cmap=\"plasma\", tiles='CartoDark', title=\"CAR from PACE OCI\")\n", "map2" ] }, { "cell_type": "code", "execution_count": null, "id": "04fa0963", "metadata": {}, "outputs": [], "source": [ "# Compare CIRE and CAR side to side\n", "map1.opts(frame_height=400, frame_width=480)+map2.opts(frame_height=400, frame_width=480)" ] }, { "cell_type": "markdown", "id": "1f307b5e-d31b-4cb4-85b1-aa236c776006", "metadata": {}, "source": [ "# Enjoy PACE data " ] }, { "cell_type": "markdown", "id": "b813f119-b9bb-4a10-8404-5ac2b2e81ef2", "metadata": {}, "source": [ "### Calculate your own vegetation index\n", "1. Select relevant bands\n", "2. Calculate index\n", "3. Plot data " ] }, { "cell_type": "markdown", "id": "a1e7cf49", "metadata": {}, "source": [ "## Contact Info: \n", "\n", "Email: LPDAAC@usgs.gov \n", "Voice: +1-866-573-3222 \n", "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹ \n", "Website: \n", "\n", "¹Work performed under USGS contract 140G0126D0001 for NASA contract NNG14HH33I. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.13.14" } }, "nbformat": 4, "nbformat_minor": 5 }