Notebook

Single Track Demo¶

Process a single ATL03 granule using SlideRule's ATL06-SR algorithm and compare the results to the existing ATL06 data product.

What is demonstrated¶

The icesat2.atl06 API is used to perform a SlideRule processing request of a single ATL03 granule
The icesat2.h5 API is used to read existing ATL06 datasets
The matplotlib package is used to plot the elevation profile of all three tracks in the granule (with the first track overlaid with the expected profile)
The cartopy package is used to produce different plots representing the geolocation of the gridded elevations produced by SlideRule.

Points of interest¶

Most use cases for SlideRule use the higher level icesat2.atl06p API which works on a region of interest; but this notebook shows some of the capabilities of SlideRule for working on individual granules.

In [1]:

import sys

import pandas as pd
import numpy as np
import math

from sliderule import icesat2

In [2]:

# Configure Session #
icesat2.init("icesat2sliderule.org", True)
resource = "_20181019065445_03150111_003_01.h5"

In [3]:

#
# ATL06 Read Request
#
def expread(resource, asset):

    # Read ATL06 Data
    delta_time  = icesat2.h5("/gt1r/land_ice_segments/delta_time",  resource, asset)
    heights     = icesat2.h5("/gt1r/land_ice_segments/h_li",        resource, asset)

    # Build Dataframe of SlideRule Responses
    df = pd.DataFrame(data=list(zip(heights, delta_time)), index=delta_time, columns=["h_mean", "delta_time"])
    df['time'] = pd.to_datetime(np.datetime64(icesat2.ATLAS_SDP_EPOCH) + (delta_time * 1000000.0).astype('timedelta64[us]'))

    # Filter Bad Elevations #
    df = df[df["h_mean"] < 4000]    
    
    # Return DataFrame
    print("Retrieved {} points from ATL06, returning {} points".format(len(heights), len(df.values)))
    return df

In [4]:

#
# SlideRule Processing Request
#
def algoexec(resource):

    # Build ATL06 Request
    parms = {
        "cnf": 4,
        "ats": 20.0,
        "cnt": 10,
        "len": 40.0,
        "res": 20.0,
        "maxi": 1
    }

    # Request ATL06 Data
    gdf = icesat2.atl06(parms, resource)

    # Return DataFrame
    print("Reference Ground Tracks: {} to {}".format(min(gdf["rgt"]), max(gdf["rgt"])))
    print("Cycle: {} to {}".format(min(gdf["cycle"]), max(gdf["cycle"])))
    print("Retrieved {} points from SlideRule".format(len(gdf["h_mean"])))
    return gdf

In [5]:

# Execute SlideRule Algorithm
act = algoexec("ATL03"+resource)

INFO:sliderule.sliderule:atl06 processing initiated on ATL03_20181019065445_03150111_003_01.h5 ...
INFO:sliderule.sliderule:... continuing to read ATL03_20181019065445_03150111_003_01.h5 (after 10 seconds)
INFO:sliderule.sliderule:... continuing to read ATL03_20181019065445_03150111_003_01.h5 (after 20 seconds)
INFO:sliderule.sliderule:... continuing to read ATL03_20181019065445_03150111_003_01.h5 (after 30 seconds)
INFO:sliderule.sliderule:processed 216303 segments in ATL03_20181019065445_03150111_003_01.h5 (after 40 seconds)
INFO:sliderule.sliderule:processing of ATL03_20181019065445_03150111_003_01.h5 complete (730524/39931/0)

Reference Ground Tracks: 315 to 315
Cycle: 1 to 1
Retrieved 622412 points from SlideRule

In [6]:

# Read ATL06 Expected Results
exp = expread("ATL06"+resource, "atlas-s3")

Retrieved 121753 points from ATL06, returning 114153 points

In [7]:

# Build Track (Actual) Datasets
standard1 = exp.sort_values(by=['time'])
track1 = act[act["spot"].isin([1, 2])]
track2 = act[act["spot"].isin([3, 4])]
track3 = act[act["spot"].isin([5, 6])]

Matplotlib Plots¶

In [8]:

# Import MatPlotLib Package
import matplotlib.pyplot as plt

In [9]:

# Create Elevation Plot
fig = plt.figure(num=None, figsize=(12, 6))

# Plot Track 1 Elevations
ax1 = plt.subplot(131)
ax1.set_title("Along Track 1 Elevations")
ax1.plot(track1.index.values, track1["h_mean"].values, linewidth=1.0, color='b')
ax1.plot(standard1["time"].values, standard1["h_mean"].values, linewidth=1.0, color='g')

# Plot Track 2 Elevations
ax2 = plt.subplot(132)
ax2.set_title("Along Track 2 Elevations")
ax2.plot(track2.index.values, track2["h_mean"].values, linewidth=1.0, color='b')

# Plot Track 3 Elevations
ax3 = plt.subplot(133)
ax3.set_title("Along Track 3 Elevations")
ax3.plot(track3.index.values, track3["h_mean"].values, linewidth=1.0, color='b')

# Show Plot
plt.show()

Cartopy Plots¶

In [10]:

import cartopy
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

In [11]:

# Create PlateCarree Plot
fig = plt.figure(num=None, figsize=(24, 12))

################################
# add global plot
################################
ax1 = plt.subplot(121,projection=cartopy.crs.PlateCarree())
ax1.set_title("Ground Tracks")

# add coastlines with filled land and lakes
ax1.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax1.add_feature(cartopy.feature.LAKES)
ax1.set_extent((-180,180,-90,90),crs=cartopy.crs.PlateCarree())

# format grid lines
gl = ax1.gridlines(crs=cartopy.crs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator([-180, -120, -60, 0, 60, 120, 180])
gl.ylocator = mticker.FixedLocator([-90, -60, -30, 0, 30, 60, 90])

# plot ground tracks
ax1.plot(track1.geometry.x, track1.geometry.y, linewidth=1.5, color='r',zorder=2, transform=cartopy.crs.Geodetic())

################################
# add zoomed plot
################################
ax2 = plt.subplot(122,projection=cartopy.crs.PlateCarree())
ax2.set_title("Ground Tracks")

# add coastlines with filled land and lakes
ax2.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax2.add_feature(cartopy.feature.LAKES)
ax2.set_extent((-80,-60,-90,-70),crs=cartopy.crs.PlateCarree())

# format grid lines
gl = ax2.gridlines(crs=cartopy.crs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

# plot ground tracks
ax2.plot(track1.geometry.x, track1.geometry.y, linewidth=1.5, color='r',zorder=2, transform=cartopy.crs.Geodetic())

# show plot
plt.show()

In [12]:

# Select Projection
#p = cartopy.crs.PlateCarree()
#p = cartopy.crs.LambertAzimuthalEqualArea()
p = cartopy.crs.SouthPolarStereo()

# Create SouthPolar Plot
fig = plt.figure(num=None, figsize=(12, 6))

ax1 = plt.axes(projection=p)
ax1.gridlines()

# Plot Ground Tracks
ax1.set_title("Ground Tracks")
ax1.plot(track1.geometry.x, track1.geometry.y,linewidth=1.5, color='r',zorder=2, transform=cartopy.crs.Geodetic())

# Plot Bounding Box
box_lon = [-70, -70, -65, -65, -70]
box_lat = [-80, -82.5, -82.5, -80, -80]
ax1.plot(box_lon, box_lat, linewidth=1.5, color='b', zorder=2, transform=cartopy.crs.Geodetic())

# add coastlines with filled land and lakes
ax1.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
ax1.add_feature(cartopy.feature.LAKES)

# add margin around plot for context
ax1.set_xmargin(1.0)
ax1.set_ymargin(1.0)

# show plot
plt.show()