Access and Visualize HyTES Data

import earthaccess
import xarray as xr
import numpy as np
import pandas as pd
import geopandas as gpd
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import panel as pn
from os import path
hv.extension('bokeh')
import warnings
warnings.filterwarnings('ignore')

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

auth = earthaccess.login()

datasets = earthaccess.search_datasets(instrument="HyTES")
print(f"HyTES datasets found: {len(datasets)}\n")
for ds in datasets:
    summary = ds.summary()
    print(f"  Short name : {summary.get('short-name', 'N/A')}")

# Load the study area boundary
roi_gdf = gpd.read_file(path.join("data", "langebaan.json"))
roi_gdf.explore(color='red', tiles=xyz, attr=attr)

# BioSCape campaign dates: October–November 2023
date_range = ("2023-10-01", "2023-11-30")

# Extract the bounding box of the ROI for earthaccess queries
# bounding_box = (min_lon, min_lat, max_lon, max_lat)
bounds = roi_gdf.total_bounds  # [minx, miny, maxx, maxy]

# Search for L1 Radiance granules over the study area
# We will reuse l1_granules in later sections, so this search covers both purposes
l1_granules = earthaccess.search_data(
    doi="10.3334/ORNLDAAC/2462",
    bounding_box=tuple(bounds),
    temporal=date_range,
)
print(f"L1 Radiance granules found: {len(l1_granules)}")
for g in l1_granules:
    print(f"  {g['meta']['native-id']}")

# Plot CMR footprints over the study area
flight_paths_gdf = gpd.GeoDataFrame(l1_granules, geometry=gpd.GeoSeries(l1_granules, crs=4326))
m = flight_paths_gdf.geometry.explore(
    tiles=xyz, attr=attr, 
    color="yellow",
    style_kwds={"weight": 2, "fillOpacity": 0.1},
    name="L1 Radiance Footprints",
)
roi_gdf.explore(m=m,color="red", name="ROI")

# csv file paths
srf_f = 'https://data.ornldaac.earthdata.nasa.gov/public/hytes/HyTES_L1_Radiance/comp/HyTES_SRF_v6.csv'
wmx_f = 'https://data.ornldaac.earthdata.nasa.gov/public/hytes/HyTES_L1_Radiance/comp/HyTES_WMX_v6.csv'
# 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()

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')

# Search for co-located L2 Emissivity and LST granules
l2_granules = earthaccess.search_data(
    doi="10.3334/ORNLDAAC/2460",
    bounding_box=tuple(bounds),
    temporal=date_range,
)
print(f"L2 Emissivity/LST granules found: {len(l2_granules)}")
for g in l2_granules:
    print(f"  {g['meta']['native-id']}")

# Open the first L1 and L2 granule from the cloud
l1_files = earthaccess.open([l1_granules[4]], provider="ORNL_CLOUD")
l2_files = earthaccess.open([l2_granules[4]], provider="ORNL_CLOUD")

print(f"L1 file: {l1_files[0].path}")
print(f"L2 file: {l2_files[0].path}")

# Open L1 HDF5 file via xarray (streaming from cloud)
ds_l1_raw = xr.open_datatree(l1_files[0], engine="h5netcdf", chunks="auto", phony_dims="sort")

# Flatten to a single dataset and rename phony dims to meaningful names
ds_l1 = ds_l1_raw.to_dataset().rename({
    "phony_dim_0": "line",
    "phony_dim_1": "sample",
    "phony_dim_2": "band",
})

# # Attach lat/lon as coordinates for easy spatial indexing
# ds_l1.coords["latitude"]  = ds_l1["latitude"]
# ds_l1.coords["longitude"] = ds_l1["longitude"]

print(f"Lines: {ds_l1.sizes['line']}, Samples: {ds_l1.sizes['sample']}, Bands: {ds_l1.sizes['band']}")
ds_l1

ds_l1["radiance_data"].isel(band=179).hvplot.image(
    data_aspect=1,
    aspect="equal",
    cmap="Greys",
    clim=(7, 12),
    frame_width=150,
    title="L1 Radiance — Band 179 (~10.8 µm)",
).opts(invert_yaxis=True)

ds_rgb = ds_l1["radiance_data"].isel(band=[150, 100, 58])
# Assign integer coords so hvplot can index the band dimension
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=150,
    title="False-color composite (R:10.1µm, G:9.2µm, B:8.5µm)",
)

# Three reference locations (lon, lat)
points_df = pd.DataFrame({
    "Type": ["Inshore", "Land", "Water"],
    "lat":  [-33.15014, -33.21875, -33.04839],
    "lon":  [ 18.06706,  18.11642,  18.00005],
})
points_gdf = gpd.GeoDataFrame(
    points_df,
    geometry=gpd.points_from_xy(points_df.lon, points_df.lat),
    crs="EPSG:4326",
)
points_gdf.explore("Type", marker_kwds={"radius": 6}, tiles=xyz, attr=attr)

def get_nearest_pixel(ds, lat, lon):
    """Return the pixel in ds whose lat/lon is closest to the given coordinates."""
    r2 = (ds["latitude"] - lat) ** 2 + (ds["longitude"] - lon) ** 2
    idx = np.where(r2.values == r2.values.min())
    return ds.sel(line=idx[0][0], sample=idx[1][0])

# assign coordinates to the xarray dataset
ds_l1.coords['longitude'] = ds_l1.longitude
ds_l1.coords['latitude'] = ds_l1.latitude

# Retrieve spectra for all three points
spectra = []
for row in points_gdf.itertuples():
    px = get_nearest_pixel(ds_l1, row.lat, row.lon)
    spectra.append(px)

pts = xr.concat(spectra, pd.Index(points_df["Type"].tolist(), name="Surface Type"))

pts["radiance_data"].hvplot.scatter(
    x="band",
    by="Surface Type",
    xlabel="Band index",
    ylabel="Radiance (W m⁻² µm⁻¹ sr⁻¹)",
    title="HyTES L1 Radiance Spectra by Surface Type",
    frame_width=500,
)

# Open L2 HDF5 file via xarray (streaming from cloud)
ds_l2_raw = xr.open_datatree(l2_files[0], engine="h5netcdf", chunks="auto", phony_dims="sort")

# The L2 file has a different dim ordering: band, line, sample, bands_emis
ds_l2 = ds_l2_raw.to_dataset().rename({
    "phony_dim_0": "band",
    "phony_dim_1": "line",
    "phony_dim_2": "sample",
    "phony_dim_3": "bands_emis",
})

print(f"L2 variables: {list(ds_l2.data_vars)}")
ds_l2

ds_l2["L2_LST"].hvplot.image(
    width=300, height=600,
    cmap="seismic",
    clim=(250, 350),
    title="L2 Land Surface Temperature (K)",
).opts(invert_yaxis=True)

# Select the inshore pixel in the L2 dataset (using fixed indices from the L1 lookup)
# L2 dimensions correspond to the same spatial grid as L1
ds_px = ds_l2.sel(sample=218, line=2164)

plot_tes = ds_px.hvplot.scatter(
    y="L2_Emissivity",
    x="L2_Emissivity_Wavelengths",
    color="red",
    label="L2_Emissivity (TES, 164 bands)",
    xlabel="Wavelength (µm)",
    ylabel="Emissivity",
)

plot_pc = ds_px.hvplot.scatter(
    y="L2_Emissivity_PC",
    x="L2_Emissivity_PC_Wavelengths",
    color="green",
    label="L2_Emissivity_PC (PC regression, 256 bands)",
)

(plot_tes * plot_pc).opts(
    title="HyTES L2 Emissivity — Inshore pixel (Langebaan)",
    frame_width=500,
    legend_position="bottom_right",
)