Access and Visualization of MASTER L2 Land Surface Temperature and Emissivity

import earthaccess
import xarray
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
import warnings
warnings.filterwarnings('ignore')

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

auth = earthaccess.login()

# Search for the MASTER L2 Facility Instrument Collection
collections = earthaccess.search_datasets(
    doi='10.3334/ORNLDAAC/2489'
)

print(f"Collections found: {len(collections)}")
for c in collections:
    # file dist. info
    fdi = c.get_umm('ArchiveAndDistributionInformation')['FileDistributionInformation'][0]
    print(f"  Title    : {c.get_umm('EntryTitle')}")
    print(f"  ShortName: {c.get_umm('ShortName')}")
    print(f"  ConceptID: {c.concept_id()}")
    print(f"  Version  : {c.get_umm('Version')}")
    print(f"  File Size  : {fdi['TotalCollectionFileSize']} {fdi['TotalCollectionFileSizeUnit'] }")
    print(f"  File format  : {fdi['Format']}")
    print(f"  Citation : {c.citation(format='apa',language="en-US")}")

# dataset short name
short_name='MASTER_L2_Emiss_LST_2489'
# project name
project = 'GEMx'

# earthaccess query
query = earthaccess.DataGranules().short_name(short_name)
query.params['campaign'] = project
granules = query.get_all()

print(f"Granules found for {project} project: {len(granules)}")

granules[0]

gdf = gpd.GeoDataFrame(pd.json_normalize(granules),
                                geometry=gpd.GeoSeries(granules, crs=4326))
gdf

gdf['dt'] = pd.to_datetime(gdf['umm.TemporalExtent.RangeDateTime.BeginningDateTime'])
# montly bins
gdf['month_year'] = gdf['dt'].dt.strftime('%Y-%m')
# Group and count occurrences per interval
gdf.groupby('month_year').size().plot(kind='bar',
                                     edgecolor='black')
plt.xlabel('Flight Dates')
plt.ylabel('Number of granules')
plt.title((
    f"""MASTER Level 2 Granules for {project} project 
    between {gdf['dt'].min().date()} and {gdf['dt'].max().date()}"""
    ))
plt.tight_layout()
plt.show()

gdf['year'] = gdf['dt'].dt.strftime('%Y')
gdf[['meta.native-id', 'year', 'geometry']].explore(
    tiles=xyz, attr=attr, column="year",
    style_kwds={"weight": 2, "fillOpacity": 0.1},
)

gdf.to_file("granule.geojson", driver="GeoJSON")

# Open the first granule as a cloud-streamed file object
file_objs = earthaccess.open([granules[0]])
granule0 = granules[0].get('umm',[]).get('GranuleUR',[])
print(f"Opened granule: {granule0}")

# xarray can read from a file-like object directly
ds = xarray.open_dataset(file_objs[0], chunks={})
ds

ds.info()

ds = ds.assign_coords(lon=ds.Lon, lat=ds.Lat)
ds

ds.LST

print(f"LST range (K):  {float(ds.LST.min()):.2f} – {float(ds.LST.max()):.2f}")

# OpenStreet baselayers
request = cimgt.OSM()

fig = plt.figure(figsize=(10, 8))
ax = plt.axes(projection=request.crs)

# baselayer zoom
ax.add_image(request, 9)

# K to C
lst_celsius = ds.LST - 273.15 
# setting min max
vmin, vmax = np.nanpercentile(lst_celsius, [2, 98])
im = lst_celsius.plot(
        ax=ax, x='lon', y='lat', 
        vmin=vmin, vmax=vmax,
        transform=ccrs.PlateCarree(),
        cmap='inferno', add_colorbar=False 
)
ax.set_title(f"MASTER L2 Land Surface Temperature {granule0}", fontsize=13)
fig.colorbar(im, ax=ax, label=f"LST (Celsius)", pad=0.05, shrink=0.7)
plt.show()

ds.Emissivity

# MASTER TIR emissivity band center wavelengths (µm) for the 6 retrieved bands
wavelengths = ds.Emissivity.attrs.get("wavelength",[])

for i, wl in enumerate(wavelengths):
    emiss = ds.Emissivity.isel({'phony_dim_2': i})
    print(f"Emissivity range for Band {wl} µm:  {float(emiss.min()):.2f} – {float(emiss.max()):.2f}")

# OpenStreet baselayers
request = cimgt.OSM()

fig, axes = plt.subplots(
    nrows=2, ncols=3, figsize=(16, 10),
    subplot_kw={'projection': request.crs})
axes = axes.flatten()

for i, wl in enumerate(wavelengths):
    emiss = ds.Emissivity.isel({'phony_dim_2': i})
    vmin, vmax = np.nanpercentile(emiss, [2, 98])
    axes[i].add_image(request, 9)

     # plot one band
    im = emiss.plot(
        ax=axes[i], x='lon', y='lat', transform=ccrs.PlateCarree(),
        cmap='RdYlGn', vmin=vmin, vmax=vmax,add_colorbar=False 
    )

    axes[i].set_title(f"Emissivity - {wl} µm", fontsize=10)   
    fig.colorbar(im, ax=axes[i], label="Emissivity", pad=0.05, shrink=0.7)
    

plt.suptitle(f"MASTER L2 Emissivity {granule0}", fontsize=13)
plt.tight_layout()
plt.show()

def get_idx(flat_idx, ds_shape, ds_dims):
    """converts flat index into multidimensional index."""
    grid_indices = np.unravel_index(flat_idx, ds_shape)
    return {dim: idx for dim, idx in zip(ds_dims, grid_indices)}

# lst layer
lst = ds.LST.compute()

# flattened index of min value
min_idx = get_idx(lst.argmin(dim=None).item(),
                  lst.shape, lst.dims)
max_idx = get_idx(lst.argmax(dim=None).item(),
                  lst.shape, lst.dims)
                   
plt.figure(figsize=(10, 5))

# emissivity at min max LST plots
ds.Emissivity.isel(min_idx).plot(
    label=f"Cold Pixel, LST {lst[min_idx]:.1f} K", 
    ls="--", marker="o", color="steelblue"
)
ds.Emissivity.isel(max_idx).plot(
    label=f"Hot Pixel, LST {lst[max_idx]:.1f} K", 
    marker="s", color="firebrick"
)
plt.xticks(ticks=np.arange(len(wavelengths)), labels=wavelengths)
plt.xlabel("Wavelength µm")
plt.title(f"MASTER L2 Emissivity {granule0}")
plt.grid(True, linestyle=":", alpha=0.6)
plt.legend()
plt.show()