# Import Packages
import os
import earthaccess
import numpy as np
import rasterio as rio
import rioxarray as rxr
import xarray as xr
import hvplot.xarray
import hvplot.pandas
import geopandas as gp
from shapely.geometry import box
import pandas as pd
import panel as pnExploring ECOSTRESS L2T LSTE
Summary
This notebook will show how to access ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) data programmatically using the earthaccess python library leaveraging NASA’s Common Metadata Repository (CMR) and enabling authentication, searching, downloading, and streaming of data with minimal coding. It also shows how to work with ECOSTRESS Tiled Land Surface Temperature and Emissivity (ECOSTRESS_L2T_LSTE) product hosted in the cloud and managed by the Land Processes Distributed Active Archive Center (LP DAAC).
Learning Objectives
- How to search ECOSTRESS data using
earthaccess - How to stream or download ECOSTRESS data
- How to clip ECOSTRESS data to a Region of Interest (ROI)
- How to quality filter ECOSTRESS data
- How to export the processed ECOSTRESS data
Requirements
- NASA Earthdata Login account. If you do not have an Earthdata Account, you can create one here.
Setup
Import the required libraries.
Authentication
Log into Earthdata using the Auth and login functions from the earthaccess library. The persist=True argument will create a local .netrc file if it doesn’t exist, or add your login info to an existing .netrc file. If no Earthdata Login credentials are found in the .netrc, you’ll be prompted for them.
auth = earthaccess.login(persist = True)
# are we authenticated?
print(auth.authenticated)Search for ECOSTRESS Data
In this example, we will use the cloud-hosted ECOSTRESS_L2T_LSTE product but the searching process can be used with other EMIT or ECOSTRESS products, other collections, or different data providers, as well as across multiple catalogs with some modification. The Land Surface Temperature and Emissivity values from ECOSTRESS Level 2 Tiled Land Surface Temperature (ECO_L2T_LSTE) are derived from five thermal infrared (TIR) bands using a physics-based Temperature and Emissivity Separation (TES) algorithm. This tiled data product uses a modified version of the Military Grid Reference System (MGRS) which divides Universal Transverse Mercator (UTM) zones into square tiles that are 109.8 km by 109.8 km with a 70 meter (m) spatial resolution.
Define Your Query Parameters
We can search for granules using attributes such as collection short name, collection ID, acquisition time, and spatial footprint.
# Define our study area: Boulder, Colorado
# These coordinates define a bounding box around Boulder
bbox = (-105.301, 39.957, -105.178, 40.094)
print(f"🎯 Study area: Boulder, Colorado")
print(f"📍 Bounding box: {bbox}")
# Create a simple polygon for our study area (for visualization)
bbox_geom = box(*bbox)
study_area = gp.GeoDataFrame([1], geometry=[bbox_geom], crs='EPSG:4326')
study_area['name'] = 'Boulder, CO'
# Let's visualize our study area
study_area.hvplot(
tiles='ESRI', geo=True, fill_alpha=0, line_color='red', line_width=2,
frame_height=400, frame_width=600, title='Study Area: Boulder, Colorado'
)Below, the parameters including provider, short_name, version, bounding_box, temporal and count are used for our query.
Next, get the downloadable links for LST and quality layers using data_links() method from earthaccess.
# Search for data using this bbox
results = earthaccess.search_data(
provider='LPCLOUD',
short_name='ECO_L2T_LSTE',
version='002',
bounding_box=bbox,
temporal=('2023-07-01','2023-08-01'),
count=100
)Next, get the downloadable links for LSTE, quality, and cloud layers using data_links() method from earthaccess.
# Get links
lst_links = [l for dl in results for l in dl.data_links() if 'LST.tif' in l]
# Show first 5
lst_links[:5]qc_links = [l for dl in results for l in dl.data_links() if 'QC.tif' in l]
qc_links[:5]cl_links = [l for dl in results for l in dl.data_links() if 'cloud.tif' in l]
cl_links[:5]Let’s take a look at the ECOSTRESS tiled data file name:
Filename: **ECOv002_L2T_LSTE_28527_009_13TDE_20230718T081442_0710_01_LST.tif**
ECO : Sensor
v002 : Product Version
L2T : Processing Level and Type (T = Tile)
LSTE : Geophysical Parameter
28527 : Orbit Number
009 : Scene ID
13TDE : Military Grid Reference System (MGRS) Tile ID
20230718 : Date of Acquisition (YYYYMMDD)
T081442 : Time of Acquisition (HHMMSS) (in UTC)
0710 : Build ID of software that generated product, Major+Minor (2+2 digits)
01 : Product Iteration Number
LST : Layer/band Name (each layer is a separate file)
.tif : Data Format for Tile Looking at Military Grid Reference System (MGRS) Tile ID of the outputs, they are all all in UTM Zone 13.
Accessing ECOSTRESS L2T Land Surface Temperature
ECOSTRESS data is stored in NASA’s Earthdata Cloud and can be accessed in different ways.
Downloading – This has been a supported option since the inception of NASA’s DAACs. Users can use the data link(s) to download files to their local working environment. This method works in both cloud and non-cloud environments.
Streaming – Streaming enables on-the-fly reading of remote files (i.e., files not saved locally). However, the accessed data must fit into the workspace’s memory. Streaming works in both cloud and non-cloud environments. Streaming data stored in the cloud without downloading is called in-place access or direct S3 access, this is only available when working in a cloud environment deployed in AWS us-west-2.
The Python libraries used to access COG files in Earthdata Cloud leverage GDAL’s virtual file systems. Whether you are running this code in the Cloud or in a local workspace, GDAL configurations must be set in order to successfully access the ECOSTRESS COG files. For this exercise, we are going to open up a context manager for the notebook using the rasterio.env module to store these configurations. The context manager sends this information, including an authentication token or cookie when connecting to a file and can also customize how the file is handled locally. A list of all available config options can be found in the GDAL config options documentation.
While the context manager is open (env.enter()) we will be able to run the open or get data commands that would typically be executed within a “with” statement. Entering the context manager for multiple cells of the notebook allows us to more freely interact with the data. We’ll close the context manager (env.exit()) when we have all of the data loaded into memory.
In this example, we will show how to stream the data. For that, the gdal configuration is set and the one of our LSTE files is read into the workspace using open_rasterio from the rioxarray library. Since the file consists of only 1 layer, we can squeeze it, removing the band dimension.
rio_env = rio.Env(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))
rio_env.__enter__()eco_lst_ds = rxr.open_rasterio(lst_links[11]).squeeze('band', drop=True)
eco_lst_dsAs mentioned, the ECOSTRESS product we are using here is tiled and the CRS is dependent on UTM zone. For this tile, we can look at the spatial_ref variable through the interactive object above to see details such as the well-known-text (WKT) representation of the CRS and other attributes. we are using hvplot for visualization here. For detailed information on available open-source Python tools and libraries for data visualization see https://pyviz.org/.
Now let’s plot the data using hvplot. reproject function is applied only for the visualization. Make sure to specify the CRS argument within the hvplot.image function so the ESRI imagery RBG background tile aligns with our scene.
size_opts = dict(frame_height=405, frame_width=720, fontscale=2)
eco_lst_ds.rio.reproject('EPSG:4326').hvplot.image(x='x', y='y', **size_opts,
cmap='inferno', tiles='ESRI', xlabel='Longitude',
ylabel='Latitude', title='ECOSTRESS LST (K)',
crs='EPSG:4326')*study_area.hvplot(line_color='blue', fill_alpha=0,
crs='EPSG:4326')Cropping ECOSTRESS Data
clip function from rasterio is used to mask data outside of our region of interest. Before clipping, we need to reproject our ROI to the projection of our dataset which is UTM zone 13N.
study_areaeco_lst_ds.rio.crspolygon_reproj = study_area.to_crs(eco_lst_ds.rio.crs)polygon_reproj.crseco_lst_roi = eco_lst_ds.rio.clip(polygon_reproj.geometry.values, polygon_reproj.crs, all_touched=True)eco_lst_roi.hvplot.image(
geo=True,cmap='inferno',**size_opts, tiles='ESRI',alpha=0.8,
title='Cropped ECOSTRESS LST (K)', xlabel='Longitude',ylabel='Latitude',
crs='EPSG:32613')Next, we will repeat the same process for the associated quality layer.
eco_qc_ds = rxr.open_rasterio(qc_links[14]).squeeze('band', drop=True)
eco_qc_roi = eco_qc_ds.rio.clip(polygon_reproj.geometry.values, polygon_reproj.crs, all_touched=True) # assign a different value to fill value
eco_qc_roi.hvplot.image(
geo=True, cmap='inferno',**size_opts, tiles='ESRI',alpha=0.8,
title='Cropped ECOSTRESS LST (K)', xlabel='Longitude',ylabel='Latitude',
crs='EPSG:32613')Quality Filtering
The quality values are 16 digits bit values with bits 0 and 1 being the mandatory QA flag that will be interpreted as:
00 = Pixel produced, best quality
01 = Pixel produced, nominal quality. Either one or more of the following conditions are met:
1. emissivity in both bands 4 and 5 < 0.95, i.e. possible cloud contamination
2. low transmissivity due to high water vapor loading (<0.4), check PWV values and error estimates
3. Pixel falls on missing scan line in bands 1&5, and filled using spatial neural net. Check error estimates
Recommend more detailed analysis of other QC information
10 = Pixel produced, but cloud detected
11 = Pixel not produced due to missing/bad data, user should check Data quality flag bits The detailed quality information is provided in Table 3-5 in ECOSTRESS Product Specification Document.
Below, the unique quality values are extracted from the clipped data and only the values showing good quality is kept.
quality_vals = np.unique(eco_qc_roi.values).tolist()
good_q = [q for q in quality_vals if np.binary_repr(q, width=16)[-2:] == '00']
good_q.where method is used to filter the quality and keep only the LSTE values generated with the best quality.
eco_lst_roi_mask = eco_lst_roi.where(eco_qc_roi.isin(good_q))
eco_lst_roi_maskImportant: There is a known data issue related to quality of ECOSTRESS Tile LSTE version 2 data:
All users of ECOSTRESS L2 v002 products (ECO_L2T_LSTE, ECO_L2_LSTE, ECO_L2G_LSTE) should be aware that the cloud mask information previously available in the Quality Control (QC) layer in v001, is not available in the v002 QC layer. Instead, users should be using the ‘cloud_mask’ layer in the L2 LSTE product, or the cloud information in the standard cloud mask products (ECO_L2_CLOUD, ECO_L2T_CLOUD, ECO_L2G_CLOUD) to assess if a pixel is clear or cloudy (see section 3 of the User Guide).
Learn more at the product CMR page
To identify clear pixels, you must now apply an extra filtering step.
In the cloud mask layer, value zero means absence of cloud, and one means presence of cloud.
When processing, ensure you retain only pixels with a value of 0 in the cloud mask layer to exclude cloudy pixels.
eco_cl_ds = rxr.open_rasterio(cl_links[14]).squeeze('band', drop=True)
eco_cl_roi = eco_cl_ds.rio.clip(polygon_reproj.geometry.values, polygon_reproj.crs, all_touched=True)
cloud_vals = np.unique(eco_cl_roi.values).tolist()
cloud_valseco_lst_roi_mask_cloud = eco_lst_roi_mask.where(eco_cl_roi.isin([0]))eco_lst_roi_mask_cloud.hvplot.image(
geo=True,cmap='inferno',**size_opts, tiles='ESRI',alpha=0.9,
title='Quality Masked ECOSTRESS LST (K)', xlabel='Longitude',ylabel='Latitude',
crs='EPSG:4326')Interactive Map Exploration of Land Surface Temperatures
def create_simple_thermal_explorer():
"""Create a simplified thermal explorer with C/F temperature display"""
# Temperature unit selector
temp_unit = pn.widgets.RadioButtonGroup(
name="Temperature Unit",
options=['Celsius', 'Fahrenheit', 'Kelvin'],
value='Celsius',
button_type='primary'
)
# Base data - already quality filtered with your approach
base_data = eco_lst_roi_mask # Already quality filtered
# Convert LST data based on selected unit
def get_converted_data(unit):
if unit == 'Celsius':
converted_data = base_data - 273.15 # K to C
unit_label = '°C'
cmap = 'inferno' # Keeping your original colormap
elif unit == 'Fahrenheit':
converted_data = (base_data - 273.15) * 9/5 + 32 # K to F
unit_label = '°F'
cmap = 'inferno'
else: # Kelvin (original)
converted_data = base_data
unit_label = 'K'
cmap = 'inferno'
return converted_data, unit_label, cmap
# Create the plot
def create_plot(unit):
converted_data, unit_label, cmap = get_converted_data(unit)
# Plot settings - using your size_opts
size_opts = dict(frame_height=500, frame_width=800, fontscale=1.2)
# Create the main plot - using geo=True for proper projection
plot = converted_data.hvplot.image(
geo=True, # Important for geographic display
cmap=cmap,
title=f'ECOSTRESS Land Surface Temperature ({unit_label})',
xlabel='Longitude',
ylabel='Latitude',
colorbar=True,
tiles='ESRI', # Keep background tiles
alpha=0.9, # Semi-transparency
crs='EPSG:4326', # Ensure correct CRS
**size_opts
)
return plot
# Create statistics panel
def create_simple_stats(unit):
converted_data, unit_label, _ = get_converted_data(unit)
# Calculate stats on valid (non-NaN) data only
valid_data = converted_data.values[~np.isnan(converted_data.values)]
if len(valid_data) > 0:
stats = {
'Minimum': f"{np.min(valid_data):.1f} {unit_label}",
'Maximum': f"{np.max(valid_data):.1f} {unit_label}",
'Mean': f"{np.mean(valid_data):.1f} {unit_label}",
'Standard Deviation': f"{np.std(valid_data):.1f} {unit_label}",
'Valid Pixels': f"{len(valid_data)} pixels",
'Coverage': f"{(len(valid_data) / (converted_data.shape[0] * converted_data.shape[1]) * 100):.1f}%"
}
else:
stats = {'Note': 'No valid data found after filtering'}
stats_html = "<div style='font-family: Arial; padding: 10px;'>"
stats_html += f"<h3>Temperature Statistics ({unit_label}):</h3>"
for key, value in stats.items():
stats_html += f"<p><strong>{key}:</strong> {value}</p>"
stats_html += "</div>"
return pn.pane.HTML(stats_html, width=300)
# Interactive functions
@pn.depends(temp_unit.param.value)
def interactive_plot(unit_value):
return create_plot(unit_value)
@pn.depends(temp_unit.param.value)
def interactive_stats(unit_value):
return create_simple_stats(unit_value)
# Layout
explorer = pn.Column(
"# ECOSTRESS Thermal Explorer",
"### Quality-filtered LST with temperature unit conversion",
temp_unit,
pn.Row(
interactive_plot,
pn.Column(
interactive_stats,
pn.pane.Markdown("""
**About the data:**
- Only using highest quality pixels (QA = '00')
- Land Surface Temperature from ECOSTRESS
- Toggle between temperature units as needed
""")
)
)
)
return explorer
# Create and display the explorer
simple_explorer = create_simple_thermal_explorer()
simple_explorer.servable()
simple_explorer2.3.3 Writing Outputs
We now have a ECOSTRESS scene that is clipped to our ROI with only good quality values. Finally, we can save this file locally.
out_name = f"../data/{lst_links[14].split('/')[-1].split('.tif')[0]}_clipped.tif"
eco_lst_roi_mask.rio.to_raster(raster_path=out_name, driver='COG')Contact Info:
Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://www.earthdata.nasa.gov/centers/lp-daac
¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.