ECOSTRESS Workshop Fall 2022

Gregory Halverson Jet Propulsion Laboratory, California Institute of Technology

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, and was sponsored by ECOSTRESS and the National Aeronautics and Space Administration (80NM0018D0004).

© 2022. All rights reserved.

Importing Libraries

These are some built-in Python functions we need for this notebook, including functions for handling filenames and dates.

from os.path import join, abspath, basename, splitext
from glob import glob
from datetime import datetime, date, timedelta
from zipfile import ZipFile
import warnings

We’re using the rioxarray package for loading raster data from GeoTIFF files, and we’re importing it as rxr. We’re using the numpy library to handle arrays, and we’re importing it as np. We’re using the rasterstats package for zonal statistics.

!mamba install -q -y rasterstats

import rioxarray as rxr
import numpy as np
from rasterstats import zonal_stats

We’re using the geopandas library to load vector data from GeoJSON files, and we’re importing it as gpd. We’re using the shapely library to handle vector data and the pyproj library to handle projections.

import geopandas as gpd
from shapely.geometry import Point, box
from shapely.ops import transform
from pyproj import Transformer

We’re using the pandas library to handle tables, and we’re importing it as pd.

import pandas as pd

We’re using the seaborn library to produce our graphs, and we’re importing it as sns. We’re using the hvplot library to produce our maps. We’re using the matplotlib library to handle plotting figures, and we’re importing it as plt.

import seaborn as sns
import hvplot.xarray
import hvplot.pandas
import matplotlib.pyplot as plt

Defining Constants

These constants define the dimensions of our figures. Feel free to adjust these to fit your display.

FIG_WIDTH_PX = 1080
FIG_HEIGHT_PX = 720
FIG_WIDTH_IN = 16
FIG_HEIGHT_IN = 9
FIG_ALPHA = 0.7
BASEMAP = "ESRI"
SEABORN_STYLE = "whitegrid"
sns.set_style(SEABORN_STYLE)

This is the location of the example ECOSTRESS Collection 2 product files.

DATA_DIRECTORY = "/home/jovyan/shared/2022-fall-ecostress-workshop/Carpinteria_ECOSTRESS_Collection_2"
print(f"data directory: {DATA_DIRECTORY}")

Loading an ECOSTRESS Collection 2 Granule

First, let’s trying opening a data layer from a product file.

filename = join(DATA_DIRECTORY, "ECOv002_L2T_LSTE_12139_005_11SKU_20200826T191453_0700_01.zip")
print(f"example L2T LSTE file: {filename}")

The granule ID for this granule can be parsed from the filename by dropping the .zip extension.

granule_ID = splitext(basename(filename))[0]
print(f"granule ID: {granule_ID}")

This product bundle, stored in zip format, contains a number of files, including raster data layers in GeoTIFF format as .tif files, and GeoJPEG browse images as .jpeg files. The GeoTIFF files can be loaded into GIS software, such as QGIS and ArcGIS. The GeoJPEG files can be loaded into Google Earth.

with ZipFile(filename) as zip_file:
    for internal_file in zip_file.filelist:
        print(internal_file.filename)

The ECOSTRESS Collection 2 tiled products include metadata in JSON format as a .json text file.

with ZipFile(filename) as zip_file:
    metadata = zip_file.read(f"{granule_ID}/{granule_ID}.json").decode()

print(metadata)

To open the temperature layer of this file, we’ll form a Universion Resource Identifier (URI) with the pattern: zip://{filename}!/{granule_ID}/{granule_ID}_{variable}.tif

URI = f"zip://{abspath(filename)}!/{granule_ID}/{granule_ID}_LST.tif"
print(f"URI: {URI}")

We’re using rioxarray to open the surface temperature layer from the L2T_LSTE product on the 11SKU tile covering the Carpinteria Salt Marsh. We’re passing the URI pointing to the GeoTIFF file contained within this zip file. If you unzip this product bundle or download the GeoTIFF file on its own, you can pass the filename of the GeoTIFF file directly into rioxarray.

ST_K_raster = rxr.open_rasterio(URI).squeeze('band', drop=True)
ST_K_raster

This xarray.DataArray object contains both an array of image values and spatial metadata. The rioxarray package extends xarray with a .rio attribute containing the metadata. Here we’re examining the coordinate reference system (CRS) of this image in the rioxarray metadata.

CRS = ST_K_raster.rio.crs
CRS.to_wkt()

This image is projected in UTM zone 11 north. Distances in this projection system are in meters. All tiled ECOSTRESS products are projected in UTM, but tiles are projected into different UTM zones, depending on where they are in the world.

We can also check the spatial resolution of the grid cells in this image in the .rio metadata.

cell_width, cell_height = ST_K_raster.rio.resolution()
print(f"cell width: {cell_width} cell height: {cell_height}")

This image is projected with 70m square pixels, as are all tiled ECOSTRESS products.

To know the observation time for this granule, we’re parsing the timestamp from the filename. This timestamp is given as Coordinated Universal Time (UTC).

datetime_UTC = datetime.strptime(basename(filename).split("_")[-3], "%Y%m%dT%H%M%S")
print(f"date/time UTC: {datetime_UTC:%Y-%m-%d %H:%M:%S}")

We want to know the centroid coordinate of this tile so that we can adjust the UTC time given to solar apparent time. We’re calculating the centroid as the average of x coordinate values and average of y coordinate values.

centroid_UTM = Point(np.nanmean(ST_K_raster.x), np.nanmean(ST_K_raster.y))
print(f"centroid UTM: {centroid_UTM}")

This centroid coordinate is in meters. We want to convert these projected x and y values in meters to geographic latitude and longitude in degress.

centroid_latlon = transform(Transformer.from_crs(CRS, "EPSG:4326", always_xy=True).transform, centroid_UTM)
print(f"centroid lat/lon: {centroid_latlon.wkt}")

We’re shifting the UTC time to local time according to longitude.

datetime_solar = datetime_UTC + timedelta(hours=(np.radians(centroid_latlon.x) / np.pi * 12))
print(f"date/time solar: {datetime_solar:%Y-%m-%d %H:%M:%S}")

The hvplot package extends xarray to allow us to plot maps. We’re reprojecting the raster geographic projection EPSG 4326 to overlay on the basemap with a latitude and longitude graticule. We’re using the jet color scheme to render temperature with a rainbow of colors with red meaning hot and blue meaning cool. We’re setting the alpha to make the raster semi-transparent on top of the basemap. We’re filtering out values lower than the 2% percentile and higher than the 98% percentile to make the variation in the image more visible.

ST_CMAP = "jet"

ST_K_map = ST_K_raster.rio.reproject("EPSG:4326").hvplot.image(
    geo=True,
    cmap=ST_CMAP,
    tiles=BASEMAP,
    alpha=FIG_ALPHA,
    width=FIG_WIDTH_PX,
    height=FIG_HEIGHT_PX,
    clim=(ST_K_raster.quantile(0.02), ST_K_raster.quantile(0.98)),
    title=f"{granule_ID} Surface Temperature (Kelvin)"
)

ST_K_map = ST_K_map.options(xlabel="Longitude", ylabel="Latitude")
ST_K_map

The temperatures in the L2T_LSTE product are given in Kelvin. To convert them to Celsius, we subtract 273.15.

ST_C_raster = ST_K_raster.copy()
ST_C_raster.data = ST_K_raster.data - 273.15

ST_C_map = ST_C_raster.rio.reproject("EPSG:4326").hvplot.image(
    geo=True,
    cmap=ST_CMAP,
    tiles=BASEMAP,
    alpha=FIG_ALPHA,
    width=FIG_WIDTH_PX,
    height=FIG_HEIGHT_PX,
    clim=(ST_C_raster.quantile(0.02), ST_C_raster.quantile(0.98)),
    title=f"{granule_ID} Surface Temperature (Celsius)"
)

ST_C_map = ST_C_map.options(xlabel="Longitude", ylabel="Latitude")
ST_C_map

Including Vector Data in Analysis

We want to analyze the ECOSTRESS images using the Carpinteria Salt Marsh Habitat Polygons provided by the USGS. This dataset is included here in GeoJSON format, which we’ll load using the geopandas package.

landcover_latlon = gpd.read_file("landcover.geojson")
landcover_latlon

To align this vector dataset with the raster datasets, we need to project it to the UTM projection used for the rasters.

landcover_UTM = landcover_latlon.to_crs(ST_C_raster.rio.crs)
landcover_UTM

This vector dataset contains polygons classifying the surface of the Carpinteria Salt Marsh into channel, salt flat, upland, pan, and marsh.

landcover_colors = {
    "channel": "blue",
    "marsh": "yellow",
    "pan": "green",
    "salt flat": "white",
    "upland": "brown"
}

landcover_map = landcover_latlon.to_crs("EPSG:4326").hvplot.polygons(
    geo=True,
    color=landcover_UTM["type"].apply(lambda type: landcover_colors[type]),
    tiles=BASEMAP,
    alpha=FIG_ALPHA,
    width=FIG_WIDTH_PX,
    height=FIG_HEIGHT_PX,
    title="Carpinteria Salt Marsh Habitat Polygons"
)

landcover_map = landcover_map.options(xlabel="Longitude", ylabel="Latitude")
landcover_map

To compare the raster image to the vector dataset, we want to subset the raster to the bounds of the vector dataset. We’re getting the bounds of our study area in meters from the convex hull of our land-cover polygons with a 100 meter buffer.

xmin, ymin, xmax, ymax = landcover_UTM.unary_union.convex_hull.buffer(100).bounds
xmin, ymin, xmax, ymax

We’re using this bounding box obtained from the extent of the vector dataset to clip the extent of the raster dataset. Now we can look at the ECOSTRESS surface temperature over the Carpinteria Salt Marsh.

ST_C_subset = ST_C_raster.rio.clip([box(xmin, ymin, xmax, ymax)])

ST_C_subset_map = ST_C_subset.rio.reproject("EPSG:4326").hvplot.image(
    geo=True,
    cmap=ST_CMAP,
    tiles=BASEMAP,
    alpha=FIG_ALPHA,
    width=FIG_WIDTH_PX,
    height=FIG_HEIGHT_PX,
    clim=(ST_C_subset.quantile(0.02), ST_C_subset.quantile(0.98)),
    title = "Carpinteria Salt Marsh Surface Temperature (Celsius)"
)

ST_C_subset_map = ST_C_subset_map.options(xlabel="Longitude", ylabel="Latitude")
ST_C_subset_map

We’re comparing the raster dataset to the polygon dataset by calculating zonal statistics with the rasterstats package.

landcover_stats = pd.DataFrame(zonal_stats(
    landcover_UTM,
    ST_C_subset.data,
    affine=ST_C_subset.rio.transform(),
    nodata=np.nan,
    stats=["min", "median", "max"]
))

landcover_stats["type"] = landcover_UTM["type"]
landcover_stats["variable"] = "ST_C"
landcover_stats["datetime_solar"] = datetime_solar
landcover_stats = landcover_stats[["type", "variable", "datetime_solar", "min", "median", "max"]]
landcover_stats

Loading a Time-Series of ECOSTRESS Files in Bulk

We’re using the built-in glob function to search for L2T_LSTE filenames in our downloaded collection of ECOSTRESS files.

L2T_LSTE_filenames = sorted(glob(join(DATA_DIRECTORY, "*_L2T_LSTE_*.zip")))
L2T_LSTE_filenames

Let’s run our zonal statistics analysis on each of the ECOSTRESS granules we have available, and produce a table of these statistics.

salt_marsh_ST = pd.DataFrame({}, columns=["datetime_solar", "type", "ST"])

for i, filename in enumerate(L2T_LSTE_filenames):
    granule_ID = splitext(basename(filename))[0]
    print(f"({i+1}/{len(L2T_LSTE_filenames)}) {granule_ID}")
    URI = f"zip://{filename}!/{granule_ID}/{granule_ID}_LST.tif"
    ST_K_raster = rxr.open_rasterio(URI).squeeze('band', drop=True)
    ST_C_raster = ST_K_raster.copy()
    ST_C_raster.data = ST_K_raster.data - 273.15
    datetime_UTC = datetime.strptime(basename(filename).split("_")[-3], "%Y%m%dT%H%M%S")
    centroid_UTM = Point(np.nanmean(ST_K_raster.x), np.nanmean(ST_K_raster.y))
    centroid_latlon = transform(Transformer.from_crs(CRS, "EPSG:4326", always_xy=True).transform, centroid_UTM)
    datetime_solar = datetime_UTC + timedelta(hours=(np.radians(centroid_latlon.x) / np.pi * 12))
    ST_C_subset = ST_C_raster.rio.clip([box(xmin, ymin, xmax, ymax)])

    landcover_stats = pd.DataFrame(zonal_stats(
        landcover_UTM,
        ST_C_subset.data,
        affine=ST_C_subset.rio.transform(),
        nodata=np.nan,
        stats=["median"]
    ))

    landcover_stats["type"] = landcover_UTM["type"]
    landcover_stats["datetime_solar"] = datetime_solar
    landcover_stats["ST"] = landcover_stats["median"].apply(lambda value: np.nan if value is None else value)
    landcover_stats = landcover_stats[["datetime_solar", "type", "ST"]]
    salt_marsh_ST = pd.concat([salt_marsh_ST, landcover_stats])

salt_marsh_ST = salt_marsh_ST.dropna()
salt_marsh_ST

Visualizing Distributions with Boxplots

Let’s examine the distribution of ECOSTRESS surface temperatures seen in the different land-cover types within the Carpinteria Salt Marsh. To visualize these distributions, we’re plotting them as boxplots. We’re using the seaborn library to generate these boxplots.

We tend to see cooler temperatures where water settles in the channel and warmer temperatures in the upland and salt flat areas, which are above the water and lack vegetation. Moderate temperatures are seen in the vegetated marsh and pan areas in between.

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    salt_marsh_ST_overall_median = salt_marsh_ST.groupby("type").median().reset_index()[["type", "ST"]].sort_values("ST")

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))

sns.boxplot(
    data=salt_marsh_ST,
    x="type",
    y="ST",
    order=salt_marsh_ST_overall_median.type,
    palette=landcover_colors,
    ax=ax
)

ax.set(xlabel="Land Cover Type", ylabel="Median Surface Temperature (Celsius)")
yticks = range(int(np.nanmin(salt_marsh_ST.ST)), int(np.nanmax(salt_marsh_ST.ST)) + 2)
ax.set_yticks(yticks)
yticklabels = [f"{tick}°C" for tick in yticks]
ax.set_yticklabels(yticklabels)
plt.title("Carpinteria Salt Marsh Median ECOSTRESS Surface Temperature Boxplots by Land Cover")
plt.show()
plt.close(fig)

print(salt_marsh_ST_overall_median)

Seasonal Agreggate Time-Series

Let’s examine a time-series of ECOSTRESS surface temperature. To make this time-series smooth, we’ll aggregate our zonal statistics by season.

salt_marsh_ST_seasonal = salt_marsh_ST.copy()
salt_marsh_ST_seasonal["season"] = salt_marsh_ST_seasonal.apply(
    lambda row: date(int(row.datetime_solar.year), int(((row.datetime_solar.month - 1) // 3) * 3 + 1), 1), axis=1)
salt_marsh_ST_seasonal = salt_marsh_ST_seasonal[["season", "type", "ST"]]
salt_marsh_ST_seasonal = salt_marsh_ST_seasonal.dropna().groupby(
    by=["season", "type"]).median().reset_index()
salt_marsh_ST_seasonal

We’re using the seaborn package again to visualize this timeseries as a line-plot. In this line-plot, we can see that temperatures in the Carptineria Salt Marsh tend to drop in the Winter and rise in the Summer. Channel and marsh temperatures tend to be the coolest for much of the year, and upland and salt flat temperatures tend to be the warmest for much of the year.

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))

sns.lineplot(
    data=salt_marsh_ST_seasonal,
    x="season",
    y="ST",
    hue="type",
    ax=ax
)

ax.set(xlabel="Season", ylabel="Seasonal Median Surface Temperature (Celsius)")
xticks = sorted(np.unique(salt_marsh_ST_seasonal.season))
ax.set_xticks(xticks)
seasons = {1: "Winter", 4: "Spring", 7: "Summer", 10: "Fall"}
ax.set_xticklabels([f"{seasons[xtick.month]}\n{xtick.year}" for xtick in xticks])
yticks = range(int(min(salt_marsh_ST_seasonal.ST)), int(max(salt_marsh_ST_seasonal.ST)) + 1)
ax.set_yticks(yticks)
ax.set_yticklabels([f"{ytick:0.0f}°C" for ytick in yticks])

plt.legend(title="Marsh Land Type")
plt.title("Carpinteria Salt Marsh Seasonal Median ECOSTRESS Surface Temperature Timeline by Land Cover")
plt.show()
plt.close(fig)

ECOSTRESS Evapotranspiration

Now let’s repeat our time-series analysis using the ECOSTRESS evapotranspiration product. ECOSTRESS Collection 2 includes the JPL Evapotranspiration Ensemble (JET) product, which runs several models to create a balanced ET estimate. We’re searching through our set of downloaded ECOSTRESS files for L3T_JET granules using glob.

L3T_JET_filenames = sorted(glob(join(DATA_DIRECTORY, "*_L3T_JET_*.zip")))
L3T_JET_filenames

We’re running our zonal statistics time-series analysis again, but now loading daily evapotranspiration (ET) in millimeters per day from the L3T JET product.

salt_marsh_ET = pd.DataFrame({}, columns=["datetime_solar", "type", "ET"])

for i, filename in enumerate(L3T_JET_filenames):
    granule_ID = splitext(basename(filename))[0]
    print(f"({i+1}/{len(L3T_JET_filenames)}) {granule_ID}")
    URI = f"zip://{filename}!/{granule_ID}/{granule_ID}_ETdaily.tif"
    ET_raster = rxr.open_rasterio(URI).squeeze('band', drop=True)
    datetime_UTC = datetime.strptime(
        basename(filename).split("_")[-3], "%Y%m%dT%H%M%S")
    centroid_UTM = Point(np.nanmean(ST_K_raster.x), np.nanmean(ST_K_raster.y))
    centroid_latlon = transform(Transformer.from_crs(
        CRS, "EPSG:4326", always_xy=True).transform, centroid_UTM)
    datetime_solar = datetime_UTC + \
        timedelta(hours=(np.radians(centroid_latlon.x) / np.pi * 12))
    ET_subset = ET_raster.rio.clip([box(xmin, ymin, xmax, ymax)])

    landcover_stats = pd.DataFrame(zonal_stats(
        landcover_UTM,
        ET_subset.data,
        affine=ET_subset.rio.transform(),
        nodata=np.nan,
        stats=["median"]
    ))

    landcover_stats["type"] = landcover_UTM["type"]
    landcover_stats["datetime_solar"] = datetime_solar
    landcover_stats["ET"] = landcover_stats["median"].apply(lambda value: np.nan if value is None or value == 0 else value)
    landcover_stats = landcover_stats[["datetime_solar", "type", "ET"]]
    salt_marsh_ET = pd.concat([salt_marsh_ET, landcover_stats])

salt_marsh_ET = salt_marsh_ET.dropna()
salt_marsh_ET

Now let’s produce seaborn box-plots again to visualize the distribution of ECOSTRESS evapotranspiration. The non-vegetated upland and salt flat areas tend to be the driest, while the channel and marsh areas tend to be the wettest.

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    salt_marsh_ET_overall_median = salt_marsh_ET.groupby("type").median().reset_index()[["type", "ET"]].sort_values("ET")

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))

sns.boxplot(
    data=salt_marsh_ET,
    x="type",
    y="ET",
    order=salt_marsh_ET_overall_median.type,
    palette=landcover_colors,
    ax=ax
)

ax.set(xlabel="Land Cover Type", ylabel="Median Daily Evapotranspiration (mm)")
yticks = np.arange(np.round(np.nanmin(salt_marsh_ET.ET)), np.round(np.nanmax(salt_marsh_ET.ET)) + 0.2, 0.1)
ax.set_yticks(yticks)
yticklabels = [f"{tick:.1f} mm" for tick in yticks]
ax.set_yticklabels(yticklabels)
plt.title("Carpinteria Salt Marsh Median ECOSTRESS Daily Evapotranspiration Boxplots by Land Cover")
plt.show()
plt.close(fig)
print(salt_marsh_ET_overall_median)

Let’s create another seasonal aggregate, now for ECOSTRESS evapotranspiration.

salt_marsh_ET_seasonal = salt_marsh_ET.copy()
salt_marsh_ET_seasonal["season"] = salt_marsh_ET_seasonal.apply(
    lambda row: date(int(row.datetime_solar.year), int(((row.datetime_solar.month - 1) // 3) * 3 + 1), 1), axis=1)
salt_marsh_ET_seasonal = salt_marsh_ET_seasonal[["season", "type", "ET"]]
salt_marsh_ET_seasonal = salt_marsh_ET_seasonal.dropna().groupby(
    by=["season", "type"]).median().reset_index()
salt_marsh_ET_seasonal

We’re using a seaborn line-plot again to visualize the seasonal aggregate of ECOSTRESS evapotranspiration by land-cover type. The channel and pan areas become dry in the heat of the Summer but then become wet again in the Autumn.

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))

sns.lineplot(
    data=salt_marsh_ET_seasonal,
    x="season",
    y="ET",
    hue="type",
    ax=ax
)

ax.set(xlabel="Season", ylabel="Seasonal Median Daily Evapotranspiration (mm)")
xticks = sorted(np.unique(salt_marsh_ET_seasonal.season))
ax.set_xticks(xticks)
seasons = {1: "Winter", 4: "Spring", 7: "Summer", 10: "Fall"}
ax.set_xticklabels([f"{seasons[xtick.month]}\n{xtick.year}" for xtick in xticks])
yticks = np.arange(np.round(np.nanmin(salt_marsh_ET_seasonal.ET), 1), np.round(np.nanmax(salt_marsh_ET_seasonal.ET), 1) + 0.1, 0.1)
ax.set_yticks(yticks)
yticklabels = [f"{tick:.1f} mm" for tick in yticks]
ax.set_yticklabels(yticklabels)
plt.legend(title="Marsh Land Type")
plt.title("Carpinteria Salt Marsh Seasonal Median ECOSTRESS Evapotranspiration Timeline by Land Cover")
plt.show()
plt.close(fig)

Relationship Between Surface Temperature and Evapotranspiration

Now let’s compare the seasonal time-series of surface temperature and evapotranspiration. First we’re merging the surface temperature and evapotranspiration tables we’re produced into a single table using pandas.

salt_marsh_seasonal = pd.merge(salt_marsh_ST_seasonal, salt_marsh_ET_seasonal, how="inner").dropna().groupby(["season", "type"]).median().reset_index()
salt_marsh_seasonal

Let’s compare the seasonal timelines of ECOSTRESS surface temperature and evapotranspiration with one last seaborn line-plot. We’re plotting these lines together with a shared x axis for the season and separate y axes to compare temperature in Celsius to evapotranspiration in millimeters per day. We see a consistent and complementary seasonal cycle between these two variables, with cool temperatures and high evapotranspiration in the Winter and warm temperatures and low evapotranspiration in the Summer.

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))
ax.set(xlabel="Season", ylabel="Seasonal Median Surface Temperature (Celsius)")
ax.grid(True)
yticks = np.arange(np.round(np.nanmin(salt_marsh_seasonal.ST)), np.round(np.nanmax(salt_marsh_seasonal.ST)))
xticks = sorted(np.unique(salt_marsh_seasonal.season))
ax.set_xticks(xticks)
seasons = {1: "Winter", 4: "Spring", 7: "Summer", 10: "Fall"}
ax.set_xticklabels([f"{seasons[xtick.month]}\n{xtick.year}" for xtick in xticks])
ax.set_yticks(yticks)
ax.set_yticklabels([f"{int(tick)} °C" for tick in yticks])

sns.lineplot(
    data=salt_marsh_seasonal,
    x="season",
    y="ST",
    color="red",
    ax=ax,
    label="Surface Temperature"
)

plt.legend(loc="upper left")
ax2 = ax.twinx()
ax2.set(ylabel="Seasonal Median Daily Evapotranspiration (mm)")
ax2.grid(False)
yticks = np.arange(np.round(np.nanmin(salt_marsh_seasonal.ET), 1), np.round(np.nanmax(salt_marsh_seasonal.ET), 1), 0.1)
ax2.set_yticks(yticks)
ax2.set_yticklabels([f"{tick:0.1f} mm" for tick in yticks])

sns.lineplot(
    data=salt_marsh_seasonal,
    x="season",
    y="ET",
    color="blue",
    ax=ax2,
    label="Evapotranspiration"
)

plt.legend(loc="upper right")
plt.title("Carpinteria Salt Marsh Seasonal Median ECOSTRESS Surface Temperature & Evapotranspiration Timeline")
plt.show()
plt.close(fig)

Finally, let’s visualize this inverse relationship between surface temperature and evapotranspiration with a scatter-plot, again using seaborn. The surface temperature is the x-axis and the evapotranspiration is the y-axis. In this scatter of points, we see a roughly linear relationship from cool and wet points in the upper left to warm and dry points in the lower right. We’re using regplot to plot this trend-line.

fig, ax = plt.subplots(figsize=(FIG_WIDTH_IN, FIG_HEIGHT_IN))
xticks = np.arange(np.round(np.nanmin(salt_marsh_seasonal.ST)), np.round(np.nanmax(salt_marsh_seasonal.ST)) + 1)
ax.set_xticks(xticks)
ax.set_xticklabels([f"{int(tick)} °C" for tick in xticks])
yticks = np.arange(np.round(np.nanmin(salt_marsh_seasonal.ET), 1), np.round(np.nanmax(salt_marsh_seasonal.ET), 1) + 0.1, 0.1)
ax.set_yticks(yticks)
ax.set_yticklabels([f"{tick:0.1f} mm" for tick in yticks])

sns.scatterplot(
    data=salt_marsh_seasonal,
    x="ST",
    y="ET",
    hue="type",
    palette=landcover_colors,
    edgecolor="black",
    linewidth=1,
    ax=ax
)

sns.regplot(
    data=salt_marsh_seasonal,
    x="ST",
    y="ET",
    scatter=False,
    color="gray",
    ax=ax
)

ax.set(xlabel="Seasonal Median Surface Temperature (Celsius)", ylabel="Seasonal Median Daily Evapotranspiration (mm)")
plt.legend(title="Marsh Land Type")
plt.title("Carpinteria Salt Marsh Seasonal Surface Temperature and Evapotranspiration Scatterplot")
plt.show()
plt.close(fig)

To quantify this inverse relationship, first we’ll calculate the correlation between these two variables using numpy. We see an inverse relationship between seasonally aggregated ECOSTRESS surface temperature and evapotranspiration in the Carpinteria Salt Marsh with a correlation coefficient of -0.49.

np.round(np.corrcoef(x=salt_marsh_seasonal.ST, y=salt_marsh_seasonal.ET)[0][1], 2)