HyTES data - Access and Visualization

HyTES data - Access and Visualization#

Overview#

The ORNL DAAC is the DAAC of Record for many SBG-relevant airborne Facility Instruments, including the Airborne Visibile/Infred Imagging Spectrometer (AVIRIS) and Hyperspectral Thermal Emission Spectrometer (HyTES) Facility Instruments. HyTES data will soon be archived and published to NASA Earthdata through the ORNL DAAC.

In this tutorial, we will access and visualize the airborne Hyperspectral Thermal Emission Spectrometer (HyTES) data for a flight line. The data were ordered from the NASA JPL HyTES website. The HyTES data will be available through NASA Earthdata data, including Earthdata Cloud, in future.

The HyTES is an airborne imaging spectrometer with 256 spectral channels between 7.5 and 12 micrometers in the thermal infrared (TIR) part of the electromagnetic spectrum and 512 pixels cross-track. HyTES provides high spatial and high spectral resolution data on surface temperature and emissivity. HyTES acquires data in the thermal infrared (TIR). TIR data are used to measure land surface temperature (LST), which informs models of water flux from land surface through processes such as evapotranspiration

# import python modules
import xarray as xr
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import panel as pn
from os import path
import pandas as pd
import geopandas as gpd
from glob import glob
import panel.widgets as pnw
import numpy as np
from shapely.geometry import MultiPoint

# esri background basemap for maps
xyz = "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
attr = "ESRI"

HyTES Flights#

A GeoJSON of the HyTES flight lines is available at the shared directory. Let’s plot the Hytes flight lines here.

# HyTES directory
base_dir = "/home/jovyan/shared-public/hytes"
# read flight line geojson
hytes_gdf = gpd.read_file(path.join(base_dir, "HyTES.json"))
hytes_gdf.head()

# plot flight line, colored by locations
hytes_gdf.explore('Location', legend=False, style_kwds={'weight': 3}, 
                  tiles=xyz, attr=attr)

Study Area#

We will be using the boundary of Langebaan Lagoon Marine Protected Area to find the intersecting HyTES flights. Langebaan is a inshore conservation area in South Africa and is declared as a Ramsar Site.

roi_gdf = gpd.read_file(path.join(base_dir, "langebaan.json"))
roi_gdf.explore(color='red', tiles=xyz, attr=attr)

Let’s find the HyTES flight lines that intersects the above region of interest.

# intersect the region of interest with HyTES flight lines
intersect_gdf = gpd.sjoin(hytes_gdf, roi_gdf, predicate='intersects')
intersect_gdf[['Location', 'Acquisition_dt']]

# mapping the flight lines
intersect_gdf.explore('Acquisition_dt', color='red', tiles=xyz, attr=attr)

HyTES flights are available for download from the JPL HyTES Website. We placed an order and downloaded the HyTES data for a flight line at Langebaan Box9ZA for 2023-10-16. They are available to download from here https://hytes.jpl.nasa.gov/orders_complete/HyTES-1747164185996a.zip (4.7GB), and also at the shared folder as shown below:

hytes_orderid = "HyTES-1747164185996a"
hytes_dir = path.join(base_dir, hytes_orderid)
hytes_dir

# list all the files from HyTES order
hytes_f = sorted(glob(path.join(hytes_dir, "*")))
hytes_f

SRF and WMX files#

Also notice there are *.csv files that are named as HyTES_SRF_* and HyTES_WMX_*. These files provides information on Signal Response Function (SRF) and Wave Matrix Data (WMX). Let’s read these CSV files using pandas.

# csv file paths
srf_f = sorted(glob(path.join(hytes_dir, "HyTES_SRF*")))[-1]
wmx_f = sorted(glob(path.join(hytes_dir, "HyTES_WMX*")))[-1]
# open srf files
srf_df = pd.read_csv(srf_f, skiprows=1, index_col='BND')
srf_df.head()

# open wmx file
wmx_df = pd.read_csv(wmx_f, skiprows=1, index_col='BND')
wmx_df.head()

We can now plot to show how the SRF is distributed for each HyTES bands. Use the slider in the plot to scroll across the bands.

def srf_plot(n):
    """returns wmx srf values for a band"""
    return pd.concat([wmx_df.loc[n], srf_df.loc[n]], 
                     keys = ['WMX','SRF'], axis=1)

# band slider
nbands = pn.widgets.IntSlider(name='Band', start=1, end=256)
# interactive SRF plot 
temp_df = hvplot.bind(srf_plot, nbands).interactive(width=600)
temp_df.hvplot.line(x='WMX', y='SRF')

Flight Path#

The flight path along with the altitude of the aircraft is recorded in the *.kml Google KML file. Let’s open the KML file and plot the flight path with flight altitude.

# kml file
kml_f = sorted(glob(path.join(hytes_dir, "*_L1_*.kml")))[0]
# read the flight line kml file to geopandas
gdf = gpd.read_file(kml_f, driver='KML')
gdflist = []
# loop through the rows and extract line coordinates (x, y, z)
for i, r in gdf.iterrows():
    gdflist.append(gpd.GeoDataFrame(crs="EPSG:4326", 
                                    geometry=[MultiPoint(r.geometry.coords)]).explode())
# create geopandas as MultiPoint
kml_gdf = gpd.GeoDataFrame(pd.concat(gdflist, ignore_index=True) )
# create a separate column for altitude
kml_gdf['altitude'] = kml_gdf.geometry.map(lambda val: round(val.z, 4))
# plot, colored by flight altitude
kml_gdf.explore('altitude', cmap = "plasma", tiles=xyz, attr=attr)

HyTES Standard Products#

As we see above, there are two levels of HyTES products currently available for download.

  • HyTES Level 1 (L1) product: Radiance and Locational Information

  • HyTES Level 2 (L2) product: Emissivity and Land Surface Temperature.

  • HyTES Level 3 (L3) product: Multi-species gas products (not available for BioSCape)

# find all hytes L1/L2 HDF5 files
hytes_h5 = sorted(glob(path.join(hytes_dir, "*.hdf5")))
# print
for f in hytes_h5:
    print(path.basename(f))

HyTES L1 product#

L1 data file in *.hdf5 format provides calibrated HyTES data in radiance units of W/m^2/µm/sr. The data is recorded with a band interleaved by pixel (BIP) format. The product also contain locational metadata from the instruments NGDCS, and per-pixel geolocation information, namely latitude, longitude, height, and number of steps taken during ray-casting. It has three dimensions: 1) lines (n = variable), 2) sample (n = 512), and 3) bands (n = 256).

Let’s open the file and look into the variables.

# L1 hdf5 file
hytes_l1 = [f for f in hytes_h5 if "_L1" in f][0]
# open dataset using xarray
ds_l1 = xr.open_datatree(hytes_l1, engine="h5netcdf", chunks='auto', 
                         phony_dims='sort')
# plot xarray datatree
ds_l1

Plotting L1 Radiance#

Let’s plot one of the radiance bands.

# convert to xarray dataset and rename dimensions to sensible names
ds_l1 = ds_l1.to_dataset().rename({'phony_dim_0': 'line',
                                   'phony_dim_1': 'sample', 
                                   'phony_dim_2': 'band'})
# printing xarray dataset
ds_l1

# retrieve radiance and plot a single band
ds_l1.radiance_data.isel(band=179).hvplot.image(data_aspect=1, aspect='equal', 
                                               cmap='Greys', clim = (7, 12), 
                                               frame_width=120).opts(invert_yaxis=True)

Let’s create an RGB radiance composite using channels: 150 (10.1 µm), 100 (9.2 µm), and 58 (8.5 µm).

# create RGB plot
ds_rgb = ds_l1.radiance_data.isel(band=[150, 100, 58]) 
ds_rgb.coords['sample'] = np.arange(ds_rgb.sizes['sample']) 
ds_rgb.coords['line'] = np.arange(ds_rgb.sizes['line']) 
ds_rgb.coords['band'] = np.arange(ds_rgb.sizes['band'])
ds_rgb.hvplot.rgb(x='sample', y='line', bands='band', rasterize=True, flip_yaxis=True, 
                  robust=True, data_aspect=1, aspect='equal', frame_width=120)

Plot pixel radiance#

Let’s plot radiances from three locations: 1) land, 2) inshore, and 3) water.

# 3 locations
df = pd.DataFrame( {"Type": ["Inshore", "Land", "Water"],
                    "y": [-33.15014, -33.21875, -33.04839],
                    "x": [18.06706, 18.11642, 18.00005]})
# create geopandas
gdf = gpd.GeoDataFrame(df, 
                       geometry=gpd.points_from_xy(df.x, df.y), crs="EPSG:4326")
# plot 
gdf.explore('Type',marker_kwds={'radius':5}, tiles=xyz, attr=attr)

Now, lets find the pixel closest to the above coordinates, and retrieve and plot radiance values for all bands.

def get_nearest_pixel(ds, lat, lon):
    """returns a single pixel dataset"""
    y = ds.latitude - lat
    x = ds.longitude - lon
    r2 = y**2 + x**2
    i = np.where(r2 == np.min(r2))
    return ds.sel(line=i[0], sample=i[1])
# assign coordinates to the xarray dataset
ds_l1.coords['longitude'] = ds_l1.longitude
ds_l1.coords['latitude'] = ds_l1.latitude

# loop through the 3 points
rad_ds = []
types = []
for r in gdf.itertuples():
    types.append(r.Type)
    # retrieve single pixel for the point
    rad_ds.append(get_nearest_pixel(ds_l1, r.geometry.y, r.geometry.x))
# create a concatenated single xarray dataset
pts = xr.concat(rad_ds, pd.Index(types, name="Types"))
# plot radiance
pts.radiance_data.hvplot.scatter(x='band', by='Types')

HyTES L2 product#

The L2 product file in *.hdf5 format provides 3 key datasets:

  1. L2_Emissivity: Emissivity spectral data from 8.2-11.5 micrometers, retrieved on 164 Temperature Emissivity Separation (TES) bands instead of 256.

  2. L2_Emissivity_PC: Principal Component (PC) eigenvector regression emissivity data from 7.4-12 micrometers for all 256 channels.

  3. L2_LST: Land surface temperature in Kelvin. Derived from atmospherically corrected level-1 radiance data using the TES algorithm.

Note some HyTES campaigns are flown on “ER2” and others campaigns are “Twin Otter” flights. The L2 products are slightly different between these two flights. Please refer to the file description and the L2 ATBD for more information.

Let’s open the L2 PC Regression Emissivity data file.

hytes_l2 = [f for f in hytes_h5 if "_L2" in f][0]
# open dataset
ds_l2 = xr.open_datatree(hytes_l2, engine="h5netcdf", chunks='auto', phony_dims='sort')
ds_l2 = ds_l2.to_dataset().rename({'phony_dim_0': 'band','phony_dim_1': 'line', 
                      'phony_dim_2': 'sample', 'phony_dim_3': 'bands_emis'})
ds_l2

Emissivity of a single pixel#

Let’s plot L2_Emissivity and L2_Emissivity_PC values from the above inshore location in a same plot.

ds_s = ds_l2.sel(sample=218, line=2164)
((ds_s.hvplot.scatter(y='L2_Emissivity_PC', x='L2_Emissivity_PC_Wavelengths', color='green', label='L2_Emissivity_PC', xlabel='Wavelengths', ylabel='Emissivity')) 
* (ds_s).hvplot.scatter(y='L2_Emissivity', x='L2_Emissivity_Wavelengths', color='red', label='L2_Emissivity'))

L2 Land Surface Temperature#

Let’s plot the variables Land Surface Temperature (LST). The data is in units of Kelvin.

ds_l2.L2_LST.hvplot.image(width=300, height=600,cmap='seismic',
             clim = (250, 350)).opts(invert_yaxis=True)
Contents