{ "cells": [ { "cell_type": "markdown", "id": "99125908-912f-4d59-b218-2d7be593818b", "metadata": { "tags": [] }, "source": [ "# 4 Exploring Canopy Water Content (CWC) and Land Surface Temperature (LST) across Vegetation Types in the Jack and Laura Dangermond Preserve\n", "\n", "> This notebook is part of a numbered series (01-05) that focuses specifically on combining EMIT VSWIR and ECOSTRESS TIR data for integrated land surface analysis.\n", "\n", "**Authors** \n", "Christiana Ade1 and Marie Johnson1,2 \n", "1. Jet Propulsion Laboratory, California Institute of Technology \n", "2. University of Montana\n", "\n", "**Summary** \n", "In this notebook, we will examine how Canopy Water Content (CWC) derived from EMIT data and Land Surface Temperature (LST) derived from ECOSTRESS change across three different vegetation types and dates in the Jack and Laura Dangermond preserve.\n", "\n", ">Note in previous notebooks, we referenced the CWC product as Equivalent Water Thickness (EWT) because water signal can still be retrieved from non-vegetated surfaces. But, here we will be investigating vegetated surfaces and will reference it as CWC.*\n", "\n", "**Learning Objectives** \n", "- Use product time series related to a real life scenario concerning invasive species in nature preserves. \n", "\n", "**Requirements** \n", " - [NASA Earthdata Account](https://urs.earthdata.nasa.gov/home) \n", " - *No Python setup requirements if connected to the workshop cloud instance!* \n", " - **Local Only** Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder to set up a local compatible Python environment \n", " - Downloaded necessary files. This is done at the end of the [01_Finding_Concurrent_Data](01_Finding_Concurrent_Data.ipynb) notebook.\n", "\n", "**Background** \n", "\n", "In 2017, The Nature Conservancy (TNC) acquired the preserve through a generous donation from Jack and Laura Dangermond. The preserve is 24,460 acres of a former private ranch situated at Point Conception, California. Point conception represents the boundary between Northern and Southern California, which separates terrestrial, marine, and coastal ecoregions. This is one of only a few areas globally where the boundary between marine and terrestrial ecosystems exist together. The preserve is one of the last \"wild coast\" areas in Southern California and it has some of the highest biodiversity globally. The Dangermond Preserve offers a unique opportunity to study the impacts of global change; sea level rise, the intensification of wildfire and drought providing a natural laboratory.\n", "\n", "The preserve is naturally composed of oak tree forests and native grasslands; however, in recent years it faces several invasion issues from non-native grasses and iceplant (*Carpobrotus spp.*). The currently mapped invasive plant species at the Preserve include non-native grasslands, iceplant mats, and stands of black mustard (*Brassica nigra*), poison hemlock (*Conium maculatum*), thistles, and fennel (*Foeniculum vulgare*). Iceplant is one of the priority species listed under the preserve's invasive species management plan.\n", "\n", "
\n", " \n", " \n", " \n", "
Figure 1. Iceplants along California's coastline.
\n", "
\n", "Iceplant is an invasive species native to South Africa that was introduced to California in the 1500s. Because of its ability to stabilize soil in the 1950's, the Department of Transportation used iceplant to prevent erosion on roadsides. This practice lasted for 20 years before it was discontinued, however, many homeowners still use iceplant for landscaping today. Although iceplant can help prevent coastal erosion, it threatens coastal biodiversity. Iceplant takes up a large amount of soil moisture, thereby taking moisture away from other species. It can also inhibit the establishment of native grass species as a result of a high nitrate accumulation in iceplant soils. Current management efforts to eradicate iceplant at the Dangermond Preserve include cattle grazing. TNC is interested in utilizing enhanced remote sensing techniques to understand iceplant characteristics and detect iceplant locations. \n", "\n", "
\n", " \n", "
Figure 2. Invasive species locations at the preserve in 2021.
\n", "
\n", " \n", "**Exercise information**\n", "\n", "The cropped EMIT and ECOSTRESS imagery used in this notebook were created using code from previous tutorial notebooks. We have simply cropped this imagery to the Dangermond boundaries and saved them to the `data` directory. As a reference, the tutorial code needed crop EMIT and ECOSTRESS products is in `02_Working_with_EMIT_Reflectance_and_ECOSTRESS_LST.ipynb` in sections *2.2.3 Cropping EMIT data to a Region of Interest*, *2.2.4 Write an output*, *2.3.2 Cropping ECOSTRESS Data*, and *2.3.3 Writing Outputs*.\n", "\n", "**Tutorial references** \n", "01_Finding_Concurrent_Data.ipynb \n", "- In this notebook we learned how to use earthaccess to find concurrent EMIT and ECOSTRESS data\n", "- We also learned how to export a list of files and download them programmatically.\n", "\n", "02_Working_with_EMIT_Reflectance_and_ECOSTRESS_LST.ipynb\n", "- How to open and work with EMIT L2A Reflectance and ECOSTRESS L2T LSTE data\n", "- How to apply a quality mask to EMIT datasets\n", "- How to reproject and regrid data\n", "- How to crop EMIT and ECOSTRESS data\n", "\n", "03_EMIT_CWC_from_Reflectance.ipynb\n", "- Calculate the EWT of a single pixel \n", "- Calculate the EWT of a ROI\n", "\n", "\n", "**Required datasets** \n", "1. Image datasets:\n", " This imagery has already been processed and uploaded to the `data` directory using code from previous notebooks. We will use three dates in 2023 to represent different seasons: April 1st, June 29th, and September 23rd. \n", "\n", "- EMIT Canopy Water Content\n", " - EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond_cwc.tif\n", " - EMIT_L2A_RFL_001_20230629T170449_2318011_002_dangermond_cwc.tif\n", " - EMIT_L2A_RFL_001_20230923T232101_2326615_002_dangermond_cwc.tif\n", "\n", "- ECOSTRESS Land Surface Temperature\n", " - ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST_dangermond.tif\n", " - ECOv002_L2T_LSTE_28238_012_10SGD_20230629T170416_0710_01_LST_dangermond.tif\n", " - ECOv002_L2T_LSTE_29576_005_10SGD_20230923T232104_0710_01_LST_dangermond.tif\n", " \n", "\n", "2. Vector datasets: \n", " - *veg_plots_outline.geojson* : This file includes polygons of delineated vegetation plots of iceplant, non-native grassland, and coastal live oak polygons. This file was created specifically for this tutorial and the polygons were delineated using vegetation maps of Dangermond and field data collected during the [SHIFT campaign.](https://climate.nasa.gov/news/3157/california-field-campaign-is-helping-scientists-protect-diverse-ecosystems/)\n", "\n", "**Tutorial Outline** \n", "\n", "4.1 View the CWC and LST images along with the vegetation class polygons \n", "4.2 Extract the CWC Raster Values \n", "4.3 Make a boxplot of CWC Extracted Values \n", "4.4 Extract and create box plot of LST images \n", "4.5 Investigate and discuss the differences between vegetation classes and dates for both CWC and LST \n", "4.6 Interactive playground! - Create your own polygons for extraction \n", " \n", "**References** \n", "* Bossard, C. C., Randall, J. M., & Hoshovsky, M. C. (Eds.). (2000). Invasive Plants of California's Wildlands. University of California Press.\n", "* Butterfield, H.S., M. Reynolds, M.G. Gleason, M. Merrifield, B.S. Cohen, W.N. Heady, D. Cameron, T. Rick, E. Inlander, M. Katkowski, L. Riege, J. Knapp, S. Gennet, G. Gorga, K. Lin, K. Easterday, B. Leahy and M. Bell. 2019. Jack and Laura Dangermond Preserve Integrated Resources Management Plan. The Nature Conservancy. 112 pages. \n", "* Iceplant images \n", " - https://centralcoastparks.org/ice-plant-the-iconic-but-destructive-piece-of-california-coastal-landscape/\n", " - https://sanctuaries.noaa.gov/news/feb15/invasive-species.html" ] }, { "cell_type": "markdown", "id": "96ee0f1b-fb26-4fac-9d2a-ab730e581dfc", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": null, "id": "d4018f37-aa3d-43ac-b258-0b2c5b7ca7eb", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Import required packages \n", "import re\n", "import os\n", "import glob\n", "import numpy as np\n", "import pandas as pd\n", "from osgeo import gdal\n", "import rasterio as rio\n", "import rioxarray as rxr\n", "from matplotlib import pyplot as plt\n", "import hvplot.xarray\n", "import hvplot.pandas\n", "import holoviews as hv\n", "import geopandas as gp\n", "import panel as pn\n", "from datetime import datetime\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "id": "91b0f7ba-626d-4af4-a7c6-1f4b3406b1cf", "metadata": {}, "source": [ "## 4.1 View the vegetation polygons and the EMIT CWC products" ] }, { "cell_type": "markdown", "id": "4b196ef1-80e8-4204-890a-3aab843b5884", "metadata": {}, "source": [ "Let's make a list called `cwc_list` of all the previously cropped EMIT CWC images by searching for tifs in our image directory `cropped/dangermond/ewt/`. " ] }, { "cell_type": "code", "execution_count": null, "id": "fc56af35-9b7f-4c20-a46f-2a678396ae2e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# data directory - location of all images cropped to dangermond. Includes subfolders for ewt and lst. \n", "data_dir = \"../../data/\"\n", "\n", "# For a later interactive portion, load a default set of polygons\n", "my_poly_gdf = gp.read_file(\"../../data/dangermond_default_polygons.geojson\")\n", "\n", "# list all files in the ewt sub-directory of the data_dir path that end in tif\n", "cwc_list = glob.glob(os.path.join(data_dir, \"cwc\", \"*.tif\"))\n", "cwc_list" ] }, { "cell_type": "markdown", "id": "041d585a-ef87-4be6-a719-2e9f8cacbc94", "metadata": {}, "source": [ "`cwc_list` is a list of nine previously processed EMIT scenes over Dangermond. However, several of those scenes do not capture the full preserve. \n", "In these first parts of the notebook, we fill focus on three dates by making a filtered list of the cwc files called `fil_cwc_list`. We will use these dates to subset\n", "- Spring: 2023-04-01\n", "- Summer: 2023-06-29\n", "- Fall: 2023-09-23" ] }, { "cell_type": "code", "execution_count": null, "id": "d0b18ed4-3593-4b54-8bd8-be787721793b", "metadata": { "tags": [] }, "outputs": [], "source": [ "# the date and time codes strings we want to select from our larger cwc_list\n", "cwc_dates = ['20230401T203751', '20230629T170449', '20230923T232101']\n", "\n", "# Written as a loop so the files appear chronologically \n", "fil_cwc_list = [] # empty list\n", "\n", "for date in cwc_dates:\n", " # filter for appropriate dates and then add back\n", " fil_cwc_list.extend([file for file in cwc_list if date in file])\n", "fil_cwc_list" ] }, { "cell_type": "markdown", "id": "0b49a303-eca9-46fd-bfa5-c249eea25ee1", "metadata": {}, "source": [ "### Let's visualize the three CWC scenes " ] }, { "cell_type": "markdown", "id": "a576e92b-6030-4b7e-9d27-cfc1e113898e", "metadata": {}, "source": [ "Next we will write a loop to visualize the change in CWC through time. \n" ] }, { "cell_type": "code", "execution_count": null, "id": "9f0e5bff-2475-4b24-9a4b-fac5d3d2156e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Initialize an empty list to store the plots\n", "plots = []\n", "\n", "# Iterate over each CWC file and create a plot\n", "for file in fil_cwc_list:\n", " # Open the raster file\n", " ras = rxr.open_rasterio(file).squeeze('band', drop=True)\n", " \n", " # extract the date string from the file file name\n", " ras_date = re.search(r\"\\d{8}\", file).group()\n", " date_object = datetime.strptime(ras_date, \"%Y%m%d\")\n", " date_only = date_object.date()#.strftime(\"%Y-%m-%d\")\n", " \n", " # Create a plot using geoviews\n", " plot = ras.hvplot.image(x='x', y='y', geo=True, cmap='blues',\n", " tiles='ESRI', \n", " title=f\"{date_only} {ras.long_name} ({ras.units})\",\n", " xlabel='Longitude', ylabel='Latitude',\n", " frame_width=400, frame_height= 300,\n", " fontscale=1, alpha=0.7)\n", " \n", " # Add the plot to the list\n", " plots.append(plot)\n", "\n", "# Display all plots in a grid layout using Panel\n", "grid = pn.GridSpec(sizing_mode='stretch_both')\n", "for i, plot in enumerate(plots):\n", " grid[i // 2, i % 2] = plot # Adjust 3 to change the number of columns\n", "\n", "grid.servable()" ] }, { "cell_type": "markdown", "id": "77114c58-bc30-4702-a2ca-c50a3a56894f", "metadata": {}, "source": [ "
The southwest corner of Dangermond seems to have a much higher CWC than other regions. From Figure 2, we can see that this is likely a very large iceplant patch. What other patterns do you see?
" ] }, { "cell_type": "markdown", "id": "f6b9d5ea-0fc7-43f1-af8e-6bf0379272e7", "metadata": {}, "source": [ "### Let's Load the Vegetation Polygons" ] }, { "cell_type": "markdown", "id": "3e42f712-25ec-4cdf-8664-dba93cfef9b3", "metadata": {}, "source": [ "Here we will load previously delineated vegetation plots of iceplant, coastal live oak and non-native grasslands from the `veg_plots_outline.geojson` as the variable `veg_poly`. These plots come from a combination of field data and previous maps created for Dangermond preserve management. \n", "\n", "
\n", " \n", "
Figure 3. Examples of the different vegetation types in the veg_plots_outline.geojson.
\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "id": "92877cbc-4c3c-4081-9820-59860d81d4d4", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Read in veg polygons\n", "veg_poly = gp.read_file(\"../../data/veg_plots_outline.geojson\")" ] }, { "cell_type": "code", "execution_count": null, "id": "4427e834-84c5-4ad4-a909-72cf5e25ee0e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Let's load in the September image for the plot visualization \n", "cw3 = rxr.open_rasterio(fil_cwc_list[2]).squeeze('band',drop=True)" ] }, { "cell_type": "markdown", "id": "c159a6d0-7629-46f4-996e-27f885b5d089", "metadata": {}, "source": [ "Let's visualize this polygon dataset on top of the September CWC scene. On the left side, `plot_1` shows the vegetation plots in red overlayed on the CWC scene. On the right side `plot_2` has the polygons colored by vegetation class. \n" ] }, { "cell_type": "code", "execution_count": null, "id": "bf798c9a-9b24-428c-a200-569160931153", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Plot 1 with CWC September for reference \n", "plot_1 = cw3.hvplot.image(x='x',y='y',geo=True, cmap='blues',tiles='ESRI',\n", " title=f\"{cw3.long_name} ({cw3.units})\",\n", " xlabel='Longitude',ylabel='Latitude', frame_width=300,\n", " frame_height=300, fontscale=1, alpha = 0.7) * veg_poly.to_crs('EPSG:3857').hvplot(c='red',alpha=1)\n", "\n", "\n", "# Plot 2 with vegetation plot types varying by color\n", "plot_2 = veg_poly.hvplot.polygons(geo=True, c='Class', hover_cols='all', \n", " xlabel='Longitude',ylabel='Latitude', frame_width=300,\n", " frame_height=300, fontscale=1, title='Dangermond Vegetation Plots')\n", "\n", "plot_1 + plot_2\n" ] }, { "cell_type": "markdown", "id": "99b417f4-b8f1-4fc0-aeec-1ba966edfc89", "metadata": {}, "source": [ "To zoom in click the `box` or `wheel` zoom icon." ] }, { "cell_type": "markdown", "id": "ece5116d-65d2-45ae-a2a5-931b954bc3d8", "metadata": {}, "source": [ "### 4.2 Extract CWC Pixel Values from Vegetation Polygons" ] }, { "cell_type": "markdown", "id": "1e7e113d-a506-43d2-996d-534ab26ebf9d", "metadata": {}, "source": [ "In this section, we will use a previously written extract function called `extract_raster_values` from the accompanying script `jlpd_ras_funcs.py`. This script includes all the more advanced functions written for this notebook which is why we are using `import *.`. Several of the other will be referenced later in the interactive portion of the notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "d9a82b83-ee6c-4493-b982-5d8276d57ad6", "metadata": { "tags": [] }, "outputs": [], "source": [ "# This is one way to load functions written from another script\n", "from vitals.jldp_ras_funcs import *" ] }, { "cell_type": "markdown", "id": "ae1917ca-c8fb-4692-a8c2-e6b379623bf0", "metadata": {}, "source": [ "Here we will create a loop to extract the raster values from every pixel in each polygon of the `veg_plots_outline.shp` file for the spring, summer, and fall `cwc_dates`. The final output dataframe file will be called `subset_df`." ] }, { "cell_type": "code", "execution_count": null, "id": "71e651db-1568-43cd-a31b-8b4b5634b13e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Create an empty list to store extracted values in\n", "ex_df = [] \n", "\n", "# Loop through each raster in fil_cwc_list\n", "for r in fil_cwc_list:\n", " # Open the raster \n", " ras = rxr.open_rasterio(r).squeeze('band', drop=True)\n", " \n", " # Extract raster values\n", " ex_ras = extract_raster_values(ras, veg_poly)\n", " \n", " # Extract the image date from the raster itself and format correctly\n", " date_object = datetime.strptime(ras.attrs['time_coverage_end'], \"%Y-%m-%dT%H:%M:%S%z\")\n", " date_only = date_object.date()\n", " \n", " # Add date to new column in data frame\n", " ex_ras['rasDate'] = date_only\n", " \n", " # Append to list\n", " ex_df.append(ex_ras)\n", " \n", "# concatenate all list entries into a dataframe\n", "ex_cwc = pd.concat(ex_df).reset_index(drop=True)\n", " # Rename column value to CWC\n", "ex_cwc.rename(columns={'value': 'CWC'}, inplace=True)\n", "# Replace -9999 with NA in the CWC column\n", "ex_cwc['CWC'] = ex_cwc['CWC'].replace(-9999, np.nan) \n", "# Remove rows with NA values in the CWC column\n", "ex_cwc = ex_cwc.dropna(subset=['CWC']) " ] }, { "cell_type": "code", "execution_count": null, "id": "5cbacaa5-a78e-4e24-9606-864160d353b0", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Preview data\n", "ex_cwc.head()" ] }, { "cell_type": "markdown", "id": "8a70e0e8-e8af-4c47-b9a9-32c06701dd60", "metadata": {}, "source": [ "Let's make a quick boxplot of our extracted values across seasons just to check that our extract function worked. " ] }, { "cell_type": "code", "execution_count": null, "id": "fbe5771e-b688-4c0a-8fa5-0894181975ec", "metadata": { "tags": [] }, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7,3))\n", "# add boxplot \n", "b = sns.boxplot(x = 'rasDate', \n", " y = 'CWC', \n", " data = ex_cwc,\n", " ax = ax,\n", " hue = 'rasDate',\n", " legend = False,\n", " width = 0.3,\n", " palette = [\"#56B4E9\", \"#009E73\", \"#CC79A7\"]) # Create color palette\n", "b.set_ylabel(r\"Canopy Water Content g/cm$\\mathregular{^2}$\", fontsize = 12)\n", "b.set_xlabel(\"Date\", fontsize = 12)\n", "b.set_title(\"Seasonal Changes in Canopy Water Content\", fontsize = 16)\n", "\n", "# overlay points for all values\n", "b = sns.stripplot(data = ex_cwc,\n", " x = 'rasDate', \n", " y = 'CWC', \n", " ax = ax,\n", " color = \"black\", # Colours the dots\n", " linewidth = 1, # Dot outline width\n", " alpha = 0.4) # Makes them transparent" ] }, { "cell_type": "markdown", "id": "e8069b11-0927-4e7f-8db9-ed8bde41b6fc", "metadata": { "tags": [] }, "source": [ "## 4.3 Generate Boxplots of CWC Across Vegetation types" ] }, { "cell_type": "markdown", "id": "abc552d1-be76-439b-aba5-e08f4c11dc20", "metadata": {}, "source": [ "Great! We managed to extract all the pixel values, but we are more interested in investigating the differences across **vegetation plot types**. Here, we separate the above plot into three faceted plots, one for each vegetation type. The final plot is called `p_class`. \n" ] }, { "cell_type": "code", "execution_count": null, "id": "9fafd0a2-cbd5-4fd0-b5d4-9093ad1a34bc", "metadata": { "tags": [] }, "outputs": [], "source": [ "p_class = sns.catplot(\n", " x='rasDate', \n", " y='CWC', \n", " col='Class', # This will create a separate plot for each unique value in the 'Class' column\n", " data=ex_cwc, \n", " kind='box',\n", " col_wrap=3, # Adjust this depending on how many plots per row you want\n", " sharex=False, \n", " sharey=False,\n", " height=5, \n", " aspect=1,\n", " width = 0.3,\n", " hue = 'rasDate',\n", " legend = False,\n", " palette = [\"#56B4E9\", \"#009E73\", \"#CC79A7\"]\n", ")\n", "\n", "# set titles and x and y labels\n", "p_class.set_titles(size = 18)\n", "p_class.set_xlabels(\"Date\", fontsize = 18)\n", "p_class.set_ylabels(r\"Canopy Water Content g/cm$\\mathregular{^2}$\", fontsize = 18)" ] }, { "cell_type": "markdown", "id": "7de3edd0-df71-4d7a-ab7f-02bb6c88986a", "metadata": {}, "source": [ "
What changes in CWC do you see across vegetation types and dates?
" ] }, { "cell_type": "markdown", "id": "f353848c-b280-424d-8c43-e976f2c4fcf6", "metadata": {}, "source": [ "## 4.4 Extract LST values and create corresponding boxplot" ] }, { "cell_type": "markdown", "id": "98776ca6-d5f1-4488-a0ac-46825c2fa833", "metadata": {}, "source": [ "We are also interested in looking at changes in LST across time, so let's repeat a similar process to extract LST values by vegetation type. The final extracted data frame is called `final_df_lst`" ] }, { "cell_type": "markdown", "id": "92a0c588-7797-48ae-a577-be8a39be843f", "metadata": {}, "source": [ "Create a list of all LST files called `lst_list`." ] }, { "cell_type": "code", "execution_count": null, "id": "0dbff30d-ebb2-474a-b8ff-271da377dab8", "metadata": { "tags": [] }, "outputs": [], "source": [ "# list all files in the ewt sub-directory of the data_dir path that end in tif\n", "lst_list = glob.glob(os.path.join(data_dir, \"lst\", \"*.tif\"))\n", "lst_list" ] }, { "cell_type": "markdown", "id": "780337ed-de0e-4e89-a47e-c6e7f181769b", "metadata": {}, "source": [ "Create a subset of list of LST files using datetimes that are concurrent with EMIT acquisitions. The list is called `fil_lst_list`." ] }, { "cell_type": "code", "execution_count": null, "id": "4b1b9204-250a-466c-9112-8612940f3a6a", "metadata": { "tags": [] }, "outputs": [], "source": [ "# the date and time codes strings we want to select from our larger lst_list\n", "lst_dates = ['20230401T203733', '20230629T170416', '20230923T232104']\n", "\n", "# #loop so that the files appear chronologically like our lst_date list\n", "fil_lst_list = [] # empty list\n", "for date in lst_dates:\n", " # filter for appropriate dates and then add back\n", " fil_lst_list.extend([file for file in lst_list if date in file])\n", "fil_lst_list" ] }, { "cell_type": "markdown", "id": "f1110129-a0bb-4ce7-869d-4325b8e75330", "metadata": { "tags": [] }, "source": [ "### Let's visualize all of our LST scenes." ] }, { "cell_type": "code", "execution_count": null, "id": "ec68b40d-45e4-4f59-b44e-3c5cd1b4b6dd", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Initialize an empty list to store the plots\n", "plots = []\n", "\n", "# Iterate over each LST file and create a plot\n", "for file in fil_lst_list:\n", " # Open the raster file\n", " ras = rxr.open_rasterio(file).squeeze('band', drop=True)\n", " \n", " # Unlike EMIT rasters, extract date from filename using regex since it's not in metadata.\n", " ras_date = re.search(r\"\\d{8}\", file).group() # search for groups of 8 digits\n", " date_object = datetime.strptime(ras_date, \"%Y%m%d\")\n", " date_only = date_object.date()\n", " \n", " # Convert Kelvin to Celsius\n", " ras -= 273.15 \n", " \n", " # Add attributes for plotting\n", " long_name_value = \"Land Surface Temperature\"\n", " ras.attrs['long_name'] = long_name_value\n", " units_value = \"°C\" # Change the units to Celsius\n", " ras.attrs['units'] = units_value\n", " \n", " # Create a plot using geoviews\n", " plot = ras.hvplot.image(x='x', y='y', geo=True, cmap='reds',\n", " tiles='ESRI', \n", " title=f\"{date_only} {ras.long_name} ({ras.units})\",\n", " xlabel='Longitude', ylabel='Latitude',\n", " frame_width=400, frame_height= 300,\n", " fontscale=1, alpha=0.7)\n", " \n", " # Add the plot to the list\n", " plots.append(plot)\n", "\n", "# Display all plots in a grid layout using Panel\n", "grid = pn.GridSpec(sizing_mode='stretch_both')\n", "for i, plot in enumerate(plots):\n", " grid[i // 2, i % 2] = plot # Adjust 3 to change the number of columns\n", "\n", "grid.servable()" ] }, { "cell_type": "markdown", "id": "9a6aaa23-a6cc-4353-95e1-46956b3a5dde", "metadata": {}, "source": [ "### Similar to before let's extract all our raster values" ] }, { "cell_type": "code", "execution_count": null, "id": "4e61f301-3729-49f0-a19e-c7b12cae50f2", "metadata": { "tags": [] }, "outputs": [], "source": [ "# create an empty list\n", "ex_df_lst = []\n", "\n", "# loop through each raster in fil_lst_list\n", "for r in fil_lst_list:\n", " # extract values \n", " ras = rxr.open_rasterio(r).squeeze('band',drop=True)\n", " \n", " # Convert from K to C\n", " ras -= 273.15 \n", " ex_ras = extract_raster_values(ras, veg_poly)\n", " \n", " # add date to the dataframe\n", " rDate = re.search(r\"\\d{8}\", r).group()\n", " date_object = datetime.strptime(rDate, \"%Y%m%d\")\n", " date_only = date_object.date()#.strftime(\"%Y-%m-%d\")\n", " \n", " # add column \n", " ex_ras['rasDate'] = date_only\n", " # add back to data frame with all values\n", " ex_df_lst.append(ex_ras)\n", "\n", "# concatenate all extracts \n", "final_df_lst = pd.concat(ex_df_lst).reset_index(drop=True)\n", "final_df_lst.rename(columns = {'value':'LST'}, inplace = True) # Rename values column to LST\n", "final_df_lst = final_df_lst[final_df_lst['LST'] > 0 ] # Removes any value less than 0" ] }, { "cell_type": "code", "execution_count": null, "id": "09b32ff9-a77d-493f-98a2-c52a90fcf528", "metadata": { "tags": [] }, "outputs": [], "source": [ "# preview data\n", "#final_df_lst.head()" ] }, { "cell_type": "markdown", "id": "a2074dee-b919-4068-93e7-917ba95e443c", "metadata": {}, "source": [ "Let's take these values and make a box plot similar to before. " ] }, { "cell_type": "code", "execution_count": null, "id": "427d5213-5ad5-4627-b700-37595e2ea0cb", "metadata": { "tags": [] }, "outputs": [], "source": [ "p_class_lst = sns.catplot(\n", " x='rasDate', \n", " y='LST', \n", " col='Class', # This will create a separate plot for each unique value in the 'Class' column\n", " data=final_df_lst, \n", " kind='box',\n", " col_wrap=3, # Adjust this depending on how many plots per row you want\n", " sharex=False, \n", " sharey=False,\n", " height=5, \n", " aspect=1,\n", " width = 0.3,\n", " hue = 'rasDate',\n", " legend = False,\n", " palette = [\"#56B4E9\", \"#009E73\", \"#CC79A7\"]\n", ")\n", "p_class_lst.set_titles(size = 18)\n", "p_class_lst.set_xlabels(\"Date\", fontsize = 18)\n", "p_class_lst.set_ylabels(\"Land Surface Temperature C\", fontsize = 18)" ] }, { "cell_type": "markdown", "id": "c19a20b5-e4ba-45ab-9895-16bc49ca3cb5", "metadata": { "tags": [] }, "source": [ "
What changes in LST do you see across vegetation types and dates?
" ] }, { "cell_type": "markdown", "id": "c8e3771a-2ee1-4a06-8f36-1234a4a4fcc7", "metadata": {}, "source": [ "## 4.5 Look at CWC and LST together" ] }, { "cell_type": "markdown", "id": "27ac32ef-d7db-4412-a46a-35c2cb98700a", "metadata": {}, "source": [ "Often times we want to visualize two variables together, so let's look at LST and CWC side by side." ] }, { "cell_type": "code", "execution_count": null, "id": "5ac1d035-0526-4b32-a0b0-5b2f74b6fdae", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Generate subplots with boxplots for different vegetation classes comparing CWC and LST\n", "n_classes = ex_cwc['Class'].nunique()\n", "fig, axes = plt.subplots(n_classes, 2, figsize=(15, 5 * n_classes), sharex='col')\n", "\n", "# Iterate over each class and plot\n", "for i, class_name in enumerate(ex_cwc['Class'].unique()):\n", " # Filter data for each class\n", " class_subset = ex_cwc[ex_cwc['Class'] == class_name]\n", " class_lst_subset = final_df_lst[final_df_lst['Class'] == class_name]\n", "\n", " # Create boxplot for CWC on the first column\n", " sns.boxplot(\n", " x='rasDate', \n", " y='CWC', \n", " data=class_subset, \n", " ax=axes[i, 0], # Plot on the first column\n", " width=0.3, \n", " hue = 'rasDate',\n", " legend = False,\n", " palette={\"#56B4E9\", \"#009E73\", \"#CC79A7\"}\n", " )\n", "\n", " # Create boxplot for LST on the second column\n", " sns.boxplot(\n", " x='rasDate', \n", " y='LST', \n", " data=class_lst_subset, \n", " ax=axes[i, 1], # Plot on the second column\n", " width=0.3, \n", " hue = 'rasDate',\n", " legend = False,\n", " palette=[\"#56B4E9\", \"#009E73\", \"#CC79A7\"]\n", " )\n", "\n", " # Set titles and labels\n", " axes[i, 0].set_title(f'Class {class_name} - CWC')\n", " axes[i, 1].set_title(f'Class {class_name} - LST')\n", " axes[i, 0].set_xlabel('Date')\n", " axes[i, 0].set_ylabel('CWC')\n", " axes[i, 1].set_xlabel('Date')\n", " axes[i, 1].set_ylabel('LST')\n", "\n", "# Adjust layout\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "37132dd3-f4e7-4074-9fa6-6e85fa78b10f", "metadata": {}, "source": [ "
What do you notice about the changes in LST and CWC across vegetation types and dates? The CWC for iceplant is interesting because it remains very high across all dates indicating that these plants contain a lot of water throughout the year.
" ] }, { "cell_type": "markdown", "id": "2c0d25cd-7ac3-422e-b804-f8e88ccb5a17", "metadata": {}, "source": [ "## 4.6 Interactive playground" ] }, { "cell_type": "markdown", "id": "95c4a0e0-0783-4732-b44c-59f1f3a05f17", "metadata": {}, "source": [ "
\n", " Now let's have some fun with what you have learned and let you draw your own polygons to study vegetation change!\n", "
\n", "\n" ] }, { "cell_type": "markdown", "id": "f49492e2-f425-4860-859c-660065190e70", "metadata": { "tags": [] }, "source": [ "This section relies on several functions from the `jldp_ras_funcs.py` script. The interactive polygon portion is based on https://github.com/auspatious/hyperspectral-notebooks/blob/main/04_EMIT_Interactive_Polygons.ipynb, with added functions to reformat the output polygons into a GeoPandas dataframe for use in the `extract_raster_values` function. " ] }, { "cell_type": "markdown", "id": "413299db-2416-4daa-be6b-12b27ffa4ca1", "metadata": { "tags": [] }, "source": [ "Here we are limiting the number of polygons you can draw using `POLY_LIMIT = 5` and are only displaying the canopy water content." ] }, { "cell_type": "code", "execution_count": null, "id": "d72a9af8-e8f8-4efb-8b41-5e4a8ac8d17d", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Limit the number of drawn polygons\n", "POLY_LIMIT = 5" ] }, { "cell_type": "markdown", "id": "59e6d3be-13cf-4d18-afc2-8d2cb9e2d5e7", "metadata": {}, "source": [ "### Let's display a CWC image and draw some polygons!" ] }, { "cell_type": "markdown", "id": "961af9e0-475d-4e43-b348-5f95d65618c1", "metadata": {}, "source": [ "Here, we will just load the canopy water content from september as our background image. This is because we already loaded it all the way above in section 4.1 when we first visualized the polygons. We have it opened as the variable `cw3`, which we then copy to the `ds` variable for using the visualization code. " ] }, { "cell_type": "code", "execution_count": null, "id": "b4834603-5f14-4e7a-8e86-ae0b55d5f450", "metadata": { "tags": [] }, "outputs": [], "source": [ "ds = cw3" ] }, { "cell_type": "markdown", "id": "10f8c49d-07e3-42be-a98f-d81ce681d87d", "metadata": {}, "source": [ "If you wanted to change the backround image you may do so using the code below." ] }, { "cell_type": "code", "execution_count": null, "id": "a5ee4b5e-d851-4996-8cfe-ea9fb51eab56", "metadata": { "tags": [] }, "outputs": [], "source": [ "# if you wanted to change the background display do it here\n", "# for example, let's set it to the spring image (index zero in our file list)\n", "# ds = rxr.open_rasterio(fil_cwc_list[0]).squeeze('band',drop=True)" ] }, { "cell_type": "markdown", "id": "e94570b3-5d83-4d14-b75a-325e44e41105", "metadata": {}, "source": [ "Set up plotting parameters and change the color of each polygon that we will draw." ] }, { "cell_type": "code", "execution_count": null, "id": "b125ac38-0f48-4455-b12a-abcaa6376992", "metadata": { "tags": [] }, "outputs": [], "source": [ "# create color cycle and list of colors for the number of polygons \n", "color_cycle = hv.Cycle('Category10')\n", "colors = [color_cycle[i] for i in range(POLY_LIMIT)]\n", "\n", "# add map variable similar to before\n", "map = ds.hvplot.image(x='x',y='y',geo=True, cmap='blues',tiles='ESRI',\n", " title=f\"{ds.long_name} ({ds.units})\",\n", " xlabel='Longitude',ylabel='Latitude', frame_width=500,\n", " frame_height=500, fontscale=1.5, alpha = 0.7)\n", "\n", "\n", "# Set up a holoviews points array to enable plotting of the clicked points\n", "xmid = ds.x.values[int(len(ds.x) / 2)]\n", "ymid = ds.y.values[int(len(ds.y) / 2)]\n", "\n", "# create holoview polygons\n", "polygons = hv.Polygons(\n", " [],\n", " kdims=[\"Longitude\", \"latitude\"],\n", ")\n", "# stream the drawn polygons\n", "polygons_stream = hv.streams.PolyDraw(\n", " data=polygons.columns(),\n", " source=polygons,\n", " num_objects=POLY_LIMIT,\n", " styles={'fill_color': color_cycle.values[0:POLY_LIMIT]}\n", ")\n" ] }, { "cell_type": "markdown", "id": "d2873319-7073-4e40-87db-51f2e359ca6a", "metadata": {}, "source": [ "**To draw a polygon activate the Polygon Draw Tool** from the toolbar on the right hand side (fourth icon from the top). Then **double click or click and hold (depending on your version of holoviews)** somewhere on the map to start your polygon and double click to end drawing your polygon. *If the hover window is in your way disable the hover tool by clicking on the last icon in the toolbar*." ] }, { "cell_type": "code", "execution_count": null, "id": "d94411e1-8cdc-493c-882e-806f498089f6", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Plot the Map and Dynamic Map side by side\n", "(map * polygons)\n" ] }, { "cell_type": "markdown", "id": "e5cc881f-7215-4f47-b837-7f3a809ab4e1", "metadata": {}, "source": [ "We will use the `hv_stream_to_rio_geometries` function to create a geojson style list of hand-drawn polygons and then convert them to a list. This list will then be converted to a geopandas dataframe for raster extraction.\n", "\n", "By default, the cell below is commented out so the notebook can run unattended using a preselected set of default polygons. If you'd like to use your own drawn polygons instead, uncomment the cell below." ] }, { "cell_type": "code", "execution_count": null, "id": "5d8dfabe-fbd7-4f6e-84ac-9d556cd442e7", "metadata": { "tags": [] }, "outputs": [], "source": [ "# list out the geometries of the different polygons \n", "# my_geometries = hv_stream_to_rio_geometries(polygons_stream.data)\n", "# geo_contents = list(my_geometries)\n", "# my_poly_gdf = create_geodataframe(geo_contents, transform_needed=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "561f5a3c-26b6-4a48-b46a-690a3eb2448e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# View polygon gdf and geometries\n", "my_poly_gdf" ] }, { "cell_type": "code", "execution_count": null, "id": "31966803-68bc-404b-86ef-a67fc994ff34", "metadata": {}, "outputs": [], "source": [ "# Write your polygon to file\n", "my_poly_gdf.to_file(\"../../data/dangermond_user_polygons.geojson\")" ] }, { "cell_type": "markdown", "id": "33620ab1-af79-47ab-a777-477f3001fd16", "metadata": {}, "source": [ "Now use the `create_geodataframe` function from `jldp_ras_funcs.py` to format your drawn polygons similarly to our `veg_plot_outline.geojson` file. This time we will not have a class column, but will have unique id column called `poly_fid column`. This column will have values that represent the order in which you drew polygons above, starting with index 0." ] }, { "cell_type": "markdown", "id": "c77ec06b-4f5d-4220-b5de-8f8dd5870417", "metadata": { "tags": [] }, "source": [ "### Now that we have our polygons let's make a box plot for our reference raster" ] }, { "cell_type": "markdown", "id": "6e34be9a-f1c8-46ea-9471-be89c3e536b4", "metadata": {}, "source": [ "Our reference raster was variable `ds`. " ] }, { "cell_type": "code", "execution_count": null, "id": "c4a51cb0-fbf2-4211-8fe6-e40a66fdeea4", "metadata": { "tags": [] }, "outputs": [], "source": [ "# extract values with previously loaded function\n", "ex_ras = extract_raster_values(ds,my_poly_gdf)\n", "# remove all values less than 0 (which are -9999, all na values)\n", "ex_ras = ex_ras[ex_ras['value'] > 0 ] \n" ] }, { "cell_type": "markdown", "id": "f428535e-2ba7-44f2-aec7-c1428ff4ad63", "metadata": { "tags": [] }, "source": [ "Using the same plotting code from above. " ] }, { "cell_type": "code", "execution_count": null, "id": "2be43a79-a655-40e1-9836-dc4713d927f0", "metadata": { "tags": [] }, "outputs": [], "source": [ "p_class_selected = sns.catplot(\n", " y='value', \n", " col='poly_fid', # This will create a separate plot for each unique value in the 'poly_fid' column\n", " data=ex_ras, \n", " kind='box',\n", " col_wrap=3, # Adjust this depending on how many plots per row you want\n", " sharex=False, \n", " sharey=False,\n", " height=5, \n", " aspect=1,\n", " width = 0.3,\n", " hue = 'poly_fid',\n", " legend = False,\n", " palette = [\"#56B4E9\", \"#009E73\", \"#CC79A7\"]\n", ")\n", "p_class_lst.set_titles(size = 18)\n", "p_class_lst.set_xlabels(\"Date\", fontsize = 18)\n", "p_class_lst.set_ylabels(\"Canopy Water Content\", fontsize = 18)" ] }, { "cell_type": "markdown", "id": "de7937f2-35e4-4dd8-820e-472e68c95254", "metadata": {}, "source": [ "Great! Our extraction and plotting worked!" ] }, { "cell_type": "markdown", "id": "adf8b7f6-c078-4c65-ac1d-3f0a2cf12aee", "metadata": {}, "source": [ "### Now let's extract data from more CWC dates" ] }, { "cell_type": "markdown", "id": "a0bc981e-d95c-4808-a68d-6c890535961c", "metadata": {}, "source": [ "If you noticed at the beginning there were more than three CWC dates available when we printed `cwc_list`. **Let's add two more dates to our boxplots.**\n", "* 02-19-2023 (this file will say merged because it is comprised of two EMIT scenes taken on the same day) \n", "* 10-14-2023\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e18e019b-364e-4337-a53b-101376de8338", "metadata": { "tags": [] }, "outputs": [], "source": [ "# the date and time codes strings we want to select from our larger cwc_list\n", "# for the first one we are using the merged 0219 dataset\n", "cwc_dates = ['20230219T202939_2305013_002_dangermond_cwc_merged',\n", " '20230401T203751', '20230629T170449', '20230923T232101',\n", " '20231014T224006']\n", "#loop so that the files appear chronologically like our cwc_dates list\n", "fil_cwc_list = [] # empty list\n", "for date in cwc_dates:\n", " # filter for appropriate dates and then add back\n", " fil_cwc_list.extend([file for file in cwc_list if date in file])\n", "fil_cwc_list" ] }, { "cell_type": "markdown", "id": "acadad59-8a76-46ae-8a6d-02b288a8b1fa", "metadata": {}, "source": [ "### Let's visualize all five images before extracting the data" ] }, { "cell_type": "code", "execution_count": null, "id": "26927e13-f200-4af9-bb68-a438b7f0864e", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Initialize an empty list to store the plots\n", "plots = []\n", "\n", "# Iterate over each CWC file and create a plot\n", "for file in fil_cwc_list:\n", " # Open the raster file\n", " ras = rxr.open_rasterio(file).squeeze('band', drop=True)\n", " \n", " # extract the date string from the file name using string matching\n", " ras_date = re.search(r\"\\d{8}\", file).group()\n", " date_object = datetime.strptime(ras_date, \"%Y%m%d\")\n", " date_only = date_object.date()\n", " \n", " # Create a plot using geoviews\n", " plot = ras.hvplot.image(x='x', y='y', geo=True, cmap='blues',\n", " tiles='ESRI', \n", " title=f\"{date_only} {ras.long_name} ({ras.units})\",\n", " xlabel='Longitude', ylabel='Latitude',\n", " frame_width=400, frame_height= 300,\n", " fontscale=1, alpha=0.7)\n", " \n", " # Add the plot to the list\n", " plots.append(plot)\n", "\n", "# Display all plots in a grid layout using Panel\n", "grid = pn.GridSpec(sizing_mode='stretch_both')\n", "for i, plot in enumerate(plots):\n", " grid[i // 2, i % 2] = plot # Adjust 3 to change the number of columns\n", "\n", "grid.servable()" ] }, { "cell_type": "markdown", "id": "302af9e1-bb37-4144-8428-f12532cdc714", "metadata": {}, "source": [ "### Now let's extract the data for the polygons we drew! " ] }, { "cell_type": "markdown", "id": "b8cb117f-3ca8-4a1e-b926-2d031a1277df", "metadata": {}, "source": [ "You should be pretty familiar with this loop by now. We did edit the line `ex_ras = extract_raster_values(ras, my_poly_gdf)` to use our polygon geodataframe `my_poly_gdf` for raster extraction instead of `veg_poly`." ] }, { "cell_type": "code", "execution_count": null, "id": "de3a5b62-3090-4b7f-b3e0-60f370f79f74", "metadata": { "tags": [] }, "outputs": [], "source": [ "# create an empty list for the extracted values\n", "ex_df = []\n", "for r in fil_cwc_list:\n", " ## extract values\n", " ras = rxr.open_rasterio(r).squeeze('band',drop=True)\n", " \n", " ## here we replaced veg_poly with my_poly_gdf\n", " ex_ras = extract_raster_values(ras, my_poly_gdf)\n", " \n", " # add date to the dataframe\n", " date_object = datetime.strptime(ras.attrs['time_coverage_end'], \"%Y-%m-%dT%H:%M:%S%z\")\n", " date_only = date_object.date()\n", " \n", " # add column \n", " ex_ras['rasDate'] = date_only\n", " # add back to data frame with all values\n", " ex_df.append(ex_ras)\n", "\n", "ex_cwc_five = pd.concat(ex_df).reset_index(drop=True)\n", "\n", "# change the name of the extracted data column\n", "ex_cwc_five.rename(columns = {'value':'CWC'}, inplace = True)\n", "# remove all nan values from the raster\n", "ex_cwc_five = ex_cwc_five[ex_cwc_five['CWC'] != -9999 ] " ] }, { "cell_type": "markdown", "id": "74d0c0c0-3423-43ae-8052-022e5eabff0c", "metadata": {}, "source": [ "Print boxplots for all your extracted polygons. You should have 5 dates this time. " ] }, { "cell_type": "code", "execution_count": null, "id": "a3cf1d16-3804-4010-9fb5-5847f9937380", "metadata": { "tags": [] }, "outputs": [], "source": [ "p_class_sel = sns.catplot(\n", " x='rasDate', \n", " y='CWC', \n", " col='poly_fid', # This will create a separate plot for each unique value in the 'poly_fid' column\n", " data=ex_cwc_five, \n", " kind='box',\n", " col_wrap=2, # Adjust this depending on how many plots per row you want\n", " sharex=False, \n", " sharey=False,\n", " height=5, \n", " aspect=1,\n", " width = 0.3,\n", ")\n", "\n", "p_class_sel.set_titles(size = 18)\n", "p_class_sel.set_xlabels(\"Date\", fontsize = 18)\n", "p_class_sel.set_ylabels(r\"Canopy Water Content g/cm$\\mathregular{^2}$\", fontsize = 18)" ] }, { "cell_type": "markdown", "id": "3318d7e8-60e0-4440-ab52-500edb5f4db6", "metadata": {}, "source": [ "
\n", " If you have time, see if you can repeat this process for additional LST dates.\n", "\n" ] } ], "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 }