SlideRule ATL03 Subsetting: Interactive Widget¶

SlideRule is an on-demand science data processing service that runs in on Amazon Web Services and responds to REST API calls to process and return science results. SlideRule was designed to enable researchers and other data systems to have low-latency access to custom-generated, high-level, analysis-ready data products using processing parameters supplied at the time of the request.

Documentation for using SlideRule is available from the project website

Background¶

SlideRule can subset ATL03 geolocated photon height data on-demand and calculate photon classifications to suit different needs.

Jupyter and SlideRule¶

Jupyter widgets are used to set parameters for the SlideRule API.

Regions of interest for submitting to SlideRule are drawn on a ipyleaflet map.
The results from SlideRule can be displayed on the interactive ipyleaflet map along with additional contextual layers.

Load necessary packages¶

In [1]:
import warnings
# turn off warnings for demo
warnings.filterwarnings('ignore')# autoreload

from sliderule import icesat2, ipysliderule, sliderule, io, earthdata
import ipywidgets as widgets
import geopandas
import logging

%load_ext autoreload
%autoreload 2

Initiate SlideRule API¶

  • Sets the URL for accessing the SlideRule service
  • Builds a table of servers available for processing data
In [2]:
# set the url for the sliderule service
# set the logging level
icesat2.init("slideruleearth.io", organization="bathy", loglevel=logging.WARNING)

Set options for making science data processing requests to SlideRule¶

SlideRule can provide different sources for photon classification. This is useful for example in cases using the photon returns to estimate ocean or lake bathymetry, vegetated canopy heights, or ground position in the presence of a vegetated canopy.

  • ATL03 photon confidence values, based on algorithm-specific classification types for land, ocean, sea-ice, land-ice, or inland water
  • ATL08 Land and Vegetation Height product photon classification
  • Experimental YAPC (Yet Another Photon Classification) photon-density-based classification
In [3]:
# display widgets for setting SlideRule parameters
SRwidgets = ipysliderule.widgets()
SRwidgets.set_atl03_defaults()
widgets.VBox([
    SRwidgets.start_date,
    SRwidgets.end_date,
    SRwidgets.classification,
    SRwidgets.surface_type,
    SRwidgets.confidence,
    SRwidgets.quality,
    SRwidgets.land_class,
    SRwidgets.yapc_knn,
    SRwidgets.yapc_win_h,
    SRwidgets.yapc_win_x,
    SRwidgets.yapc_min_ph,
])
Out[3]:
VBox(children=(DatePicker(value=datetime.datetime(2018, 10, 13, 0, 0), description='Start Date', step=1, toolt…

Interactive Mapping with Leaflet¶

Interactive maps within the SlideRule python API are built upon ipyleaflet.

Leaflet Basemaps and Layers¶

There are 3 projections available within SlideRule for mapping (Global, North and South). There are also contextual layers available for each projection.

Global (Web Mercator, EPSG:3857) North (Alaska Polar Stereographic, EPSG:5936) South (Antarctic Polar Stereographic, EPSG:3031)
  • USGS 3DEP Hillshade
  • ASTER GDEM Hillshade
  • ESRI Imagery
  • Global Land Ice Measurements from Space (GLIMS)
  • Randolph Glacier Inventory (RGI)
  • ESRI Imagery
  • ArcticDEM
  • Landsat Image Mosaic of Antarctica (LIMA)
  • MODIS Mosaic of Antarctica (MOA)
  • Radarsat Antarctic Mapping Project (RAMP)
  • Reference Elevation Model of Antarctica (REMA)

In addition, most xyzservice providers can be added as contextual layers to the global Web Mercator maps

In [4]:
widgets.VBox([
    SRwidgets.projection,
    SRwidgets.layers,
    SRwidgets.raster_functions
])
Out[4]:
VBox(children=(Dropdown(description='Projection:', options=('Global', 'North', 'South'), tooltip='Projection: …

Select regions of interest for submitting to SlideRule¶

Here, we create polygons or bounding boxes for our regions of interest.
This map is also our viewer for inspecting our SlideRule ICESat-2 data returns.

In [5]:
# create ipyleaflet map in specified projection
m = ipysliderule.leaflet(SRwidgets.projection.value)
# install click handler callback
m.add_selected_callback(SRwidgets.atl03_click_handler)
m.map
Out[5]:
Map(center=[39, -108], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out_t…
In [6]:
m.add_layer(
    layers=SRwidgets.layers.value,
    rendering_rule=SRwidgets.rendering_rule
)

Inspect list of available granules¶

  • SlideRule will query the NASA Common Metadata Repository (CMR) for ATL03 data within our region of interest
In [7]:
%%time
# sliderule asset and data release
asset = SRwidgets.asset.value
release = SRwidgets.release.value
# find granule for each region of interest
granules_list = []
# for each region of interest
for poly in m.regions:
    granules = earthdata.cmr(short_name="ATL03",
        polygon=poly,
        time_start=SRwidgets.time_start,
        time_end=SRwidgets.time_end,
        version=release)
    granules_list.extend(granules)
# inspect granules list
print(f'Available granules: {len(granules_list)}')
Available granules: 43
CPU times: user 27.7 ms, sys: 20 ms, total: 47.7 ms
Wall time: 1 s

Transmit requests to SlideRule¶

  • When using the icesat2 asset, the ICESat-2 ATL03 data are then accessed from the NSIDC AWS s3 bucket in us-west-2
  • The ATL03 granules is spatially subset within SlideRule to our exact region of interest
  • Photon classification parameters can then be extracted or calculated for our ATL03 data
  • The completed data is streamed concurrently back and combined into a geopandas GeoDataFrame within the Python client
In [8]:
# build sliderule parameters using latest values from widget
parms = SRwidgets.build_atl03()

# clear existing geodataframe results
elevations = [sliderule.emptyframe()]

# for each region of interest
for poly in m.regions:
    # add polygon from map to sliderule parameters
    parms["poly"] = poly 
    # make the request to the SlideRule (ATL03-SR) endpoint
    # and pass it the request parameters to request ATL03 Data
    elevations.append(icesat2.atl03sp(parms, resources=granules_list))
    
gdf = geopandas.pd.concat(elevations)

Review GeoDataFrame output¶

Can inspect the columns, number of returns and returns at the top of the GeoDataFrame.

See the ICESat-2 documentation for descriptions of each column

In [9]:
print(f'Returned {gdf.shape[0]} records')
gdf.head()
Returned 3027642 records
Out[9]:
geometry solar_elevation segment_dist segment_id background_rate region sc_orient rgt pair cycle ... height atl03_cnf landcover relief atl08_class quality_ph snowcover x_atc y_atc spot
2018-11-19 20:59:13.509515520 POINT (-108.83775 39.15001) 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 ... 1662.840942 0.0 255.0 0.0 4.0 0.0 255.0 -9.930932 -5537.269043 2.0
2018-11-19 20:59:13.509515520 POINT (-108.83774 39.15002) 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 ... 1609.357056 0.0 255.0 0.0 4.0 0.0 255.0 -9.970532 -5538.292969 2.0
2018-11-19 20:59:13.509515520 POINT (-108.83771 39.15002) 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 ... 1448.442017 2.0 255.0 0.0 2.0 0.0 255.0 -10.089421 -5541.372070 2.0
2018-11-19 20:59:13.509515520 POINT (-108.83768 39.15002) 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 ... 1303.463745 0.0 255.0 0.0 4.0 0.0 255.0 -10.196534 -5544.146973 2.0
2018-11-19 20:59:13.509515520 POINT (-108.83766 39.15002) 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 ... 1219.619141 0.0 255.0 0.0 4.0 0.0 255.0 -10.258735 -5545.751465 2.0

5 rows × 22 columns

Add GeoDataFrame to map¶

For stability of the leaflet map, SlideRule will as a default limit the plot to have up to 10000 points from the GeoDataFrame

GeoDataFrames can be plotted in any available matplotlib colormap

In [10]:
widgets.VBox([
    SRwidgets.variable,
    SRwidgets.cmap,
    SRwidgets.reverse,
])
Out[10]:
VBox(children=(Dropdown(description='Variable:', index=5, options=('atl03_cnf', 'atl08_class', 'cycle', 'delta…
In [11]:
%matplotlib inline
# ATL03 fields for hover tooltip
fields = m.default_atl03_fields()
m.GeoData(gdf, column_name=SRwidgets.variable.value, cmap=SRwidgets.colormap,
    max_plot_points=10000, tooltip=True, colorbar=True, fields=fields)
Deprecated. Will be removed in a future release

Create plots for a single track¶

  • scatter: Will plot data returned by SlideRule for a single RGT, ground track and cycle
In [12]:
widgets.VBox([
    SRwidgets.plot_classification,
    SRwidgets.rgt,
    SRwidgets.ground_track,
    SRwidgets.cycle,
])
Out[12]:
VBox(children=(Dropdown(description='Classification', index=1, options=('atl03', 'atl08', 'yapc', 'none'), too…
In [14]:
%matplotlib widget
SRwidgets.plot(atl03=gdf, kind='scatter', title='Photon Cloud',
    cmap=SRwidgets.colormap, legend=True, legend_frameon=True,
    classification=SRwidgets.plot_classification.value, 
    segments=False)
Figure
No description has been provided for this image

Save GeoDataFrame to output file¶

  • pytables HDF5: easily read back as a Geopandas GeoDataFrame
  • netCDF: interoperable with other programs
In [15]:
display(SRwidgets.filesaver)
Text(value='ATL03-SR_20240404194513_006.h5')
In [16]:
# append sliderule api version to attributes
version = sliderule.get_version()
parms['version'] = version['icesat2']['version']
parms['commit'] = version['icesat2']['commit']
# save to file in format (HDF5 or netCDF)
io.to_file(gdf, SRwidgets.file,
    format=SRwidgets.format,
    driver='pytables',
    parameters=parms,
    regions=m.regions,
    verbose=True)

Read GeoDataFrame from input file¶

  • pytables HDF5
  • netCDF
In [17]:
display(SRwidgets.fileloader)
Text(value='')
In [18]:
# read from file in format (HDF5 or netCDF)
gdf,parms,regions = io.from_file(SRwidgets.file,
    format=SRwidgets.format,
    driver='pytables',
    return_parameters=True,
    return_regions=True)

Review GeoDataFrame input from file¶

In [19]:
gdf.head()
Out[19]:
solar_elevation segment_dist segment_id background_rate region sc_orient rgt pair cycle track ... atl03_cnf landcover relief atl08_class quality_ph snowcover x_atc y_atc spot geometry
2018-11-19 20:59:13.509515520 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 3.0 ... 0.0 255.0 0.0 4.0 0.0 255.0 -9.930932 -5537.269043 2.0 POINT (-108.83775 39.15001)
2018-11-19 20:59:13.509515520 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 3.0 ... 0.0 255.0 0.0 4.0 0.0 255.0 -9.970532 -5538.292969 2.0 POINT (-108.83774 39.15002)
2018-11-19 20:59:13.509515520 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 3.0 ... 2.0 255.0 0.0 2.0 0.0 255.0 -10.089421 -5541.372070 2.0 POINT (-108.83771 39.15002)
2018-11-19 20:59:13.509515520 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 3.0 ... 0.0 255.0 0.0 4.0 0.0 255.0 -10.196534 -5544.146973 2.0 POINT (-108.83768 39.15002)
2018-11-19 20:59:13.509515520 25.088615 4.355700e+06 217195.0 713556.370971 2.0 1.0 798.0 0.0 1.0 3.0 ... 0.0 255.0 0.0 4.0 0.0 255.0 -10.258735 -5545.751465 2.0 POINT (-108.83766 39.15002)

5 rows × 22 columns

Set parameters and add saved regions to map¶

In [20]:
SRwidgets.set_values(parms)
m.add_region(regions)
Out[20]:
<sliderule.ipysliderule.leaflet at 0xffff735cd940>
In [ ]: