Interactive Visualizion with Open Altimetry & Google Earth Engine#

By Philipp Arndt
Scripps Institution of Oceanography, University of California San Diego
Github: @fliphilipp
Contact: parndt@ucsd.edu

Learning Objectives#

Discover some interesting ICESat-2 data by browsing OpenAltimetry
Get the data into a Jupyter Notebook, plot it and overlay the ground track on a map
Use Google Earth Engine to find the closest-in-time Sentinel-2 image that is cloud-free along ICESat-2’s gound track
Make a plot from ICESat-2 data and contemporaneous imagery in Jupyter/matplotlib
Bonus Code: Use Google Earth Engine to extract Sentinel-2 data along the ICESat-2 ground track

Computing environment#

We’ll be using the following Python libraries in this notebook:

%matplotlib widget 
import os
import ee
import geemap
import requests
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
from datetime import datetime
from datetime import timedelta
from ipywidgets import Layout
import rasterio as rio
from rasterio import plot as rioplot
from rasterio import warp

The import below is a class that I wrote myself. It helps us read and store data from the OpenAltimetry API.
If you are interested in how this works, you can find the code in utils/oa.py.

from utils.oa import dataCollector

Google Earth Engine Authentication and Initialization#

GEE requires you to authenticate your access, so if ee.Initialize() does not work you first need to run ee.Authenticate(). This gives you a link at which you can use your google account that is associated with GEE to get an authorization code. Copy the authorization code into the input field and hit enter to complete authentication.

try:
    ee.Initialize()
except: 
    ee.Authenticate()
    ee.Initialize()

# ee.Authenticate()

Downloading data from the OpenAltimetry API#

Let’s say we have found some data that looks weird to us, and we don’t know what’s going on.

image showing OpenAltimetry screenshot

A short explanation of how I got the data:#

I went to openaltimetry.org and selected BROWSE ICESAT-2 DATA. Then I selected ATL 06 (Land Ice Height) on the top right, and switched the projection🌎 to Arctic. Then I selected August 22, 2021 in the calendar📅 on the bottom left, and toggled the cloud☁️ button to show MODIS imagery of that date. I then zoomed in on my region of interest.

To find out what ICESat-2 has to offer here, I clicked on SELECT A REGION on the top left, and drew a rectangle around that mysterious cloud. When right-clicking on that rectangle, I could select View Elevation profile. This opened a new tab, and showed me ATL06 and ATL08 elevations.

It looks like ATL06 can’t decide where the surface is, and ATL08 tells me that there’s a forest canopy on the Greenland Ice Sheet? To get to the bottom of this, I scrolled all the way down and selected 🛈Get API URL. The website confirms that the “API endpoint was copied to clipboard.” Now I can use this to access the data myself.

Note: Instead of trying to find this region yourself, you can access the OpenAltimetry page shown above by going to this annotation🏷️. Just left-click on the red box and select “View Elevation Profile”.

Getting the OpenAltimetry info into python#

All we need to do is to paste the API URL that we copied from the webpage into a string. We also need to specify which beam we would like to look at. The GT1R ground track looks funky, so let’s look at that one!

# paste the API URL from OpenAltimetry below, and specify the beam you are interested in
oa_api_url = 'http://openaltimetry.org/data/api/icesat2/atl03?date=2021-08-22&minx=-24.03772653601395&miny=77.53753839308544&maxx=-23.86438315078483&maxy=77.57604313591165&trackId=909&beamName=gt1r&beamName=gt1l&outputFormat=json'
gtx = 'gt1r'

We can now initialize a dataCollector object, using the copy-pasted OpenAltimetry API URL, and the beam we would like to look at.
(Again, I defined this class in utils/oa.py to do some work for us in the background.)

is2data = dataCollector(oaurl=oa_api_url,beam=gtx, verbose=True)

OpenAltimetry API URL: http://openaltimetry.org/data/api/icesat2/atlXX?date=2021-08-22&minx=-24.03772653601395&miny=77.53753839308544&maxx=-23.86438315078483&maxy=77.57604313591165&trackId=909&outputFormat=json&beamName=gt1r&client=jupyter
Date: 2021-08-22
Track: 909
Beam: gt1r
Latitude limits: [77.53753839308544, 77.57604313591165]
Longitude limits: [-24.03772653601395, -23.86438315078483]

Alternatively, we could use a date, track number, beam, and lat/lon bounding box as input to the dataCollector.

date = '2021-08-22'
rgt = 909
beam = 'gt1r'
latlims = [77.5326, 77.5722]
lonlims = [-23.9891, -23.9503]
is2data = dataCollector(date=date, latlims=latlims, lonlims=lonlims, track=rgt, beam=gtx, verbose=True)

OpenAltimetry API URL: https://openaltimetry.org/data/api/icesat2/atlXX?date=2021-08-22&minx=-23.9891&miny=77.5326&maxx=-23.9503&maxy=77.5722&trackId=909&beamName=gt1r&outputFormat=json&client=jupyter
Date: 2021-08-22
Track: 909
Beam: gt1r
Latitude limits: [77.5326, 77.5722]
Longitude limits: [-23.9891, -23.9503]

Note that this also constructs the API url for us.

Requesting the data from the OpenAltimetry API#

Here we use the requestData() function of the dataCollector class, which is defined in utils/oa.py. It downloads all data products that are available on OpenAltimetry based on the inputs with which we initialized our dataCollector, and writes them to pandas dataframes.

is2data.requestData(verbose=True)

---> requesting ATL03 data... 6314 data points.
---> requesting ATL06 data... 198 data points.
---> requesting ATL07 data... No data.
---> requesting ATL08 data... 26 data points.
---> requesting ATL10 data... No data.
---> requesting ATL12 data... No data.
---> requesting ATL13 data... No data.

The data are now stored as data frames in our dataCollector object. To verify this, we can run the cell below.

vars(is2data)

{'url': 'https://openaltimetry.org/data/api/icesat2/atlXX?date=2021-08-22&minx=-23.9891&miny=77.5326&maxx=-23.9503&maxy=77.5722&trackId=909&beamName=gt1r&outputFormat=json&client=jupyter',
 'date': '2021-08-22',
 'track': 909,
 'beam': 'gt1r',
 'latlims': [77.5326, 77.5722],
 'lonlims': [-23.9891, -23.9503],
 'atl03':             lat        lon          h   conf
 0     77.572273 -23.953863  748.45350  Noise
 1     77.572246 -23.953913  847.64276  Noise
 2     77.572247 -23.953895  785.39874  Noise
 3     77.572232 -23.953940  906.48040  Noise
 4     77.572223 -23.953903  739.40704  Noise
 ...         ...        ...        ...    ...
 6309  77.532606 -23.986983  788.30840   High
 6310  77.532600 -23.986988  787.72420   High
 6311  77.532594 -23.986993  788.89220   High
 6312  77.532594 -23.986993  786.84937   High
 6313  77.532594 -23.986993  787.69586   High
 
 [6314 rows x 4 columns],
 'atl06':            lat        lon          h
 0    77.572104 -23.954022  810.27400
 1    77.571928 -23.954170  809.49520
 2    77.571751 -23.954318  808.62430
 3    77.571575 -23.954465  808.51855
 4    77.571398 -23.954612  808.20930
 ..         ...        ...        ...
 193  77.533473 -23.986256  785.09015
 194  77.533297 -23.986406  785.85700
 195  77.533120 -23.986554  786.59470
 196  77.532944 -23.986702  787.29694
 197  77.532768 -23.986849  787.82263
 
 [198 rows x 3 columns],
 'atl07': Empty DataFrame
 Columns: [lat, lon, h]
 Index: [],
 'atl08':           lat        lon          h     canopy
 0   77.571838 -23.954245  809.23700   3.455017
 1   77.570961 -23.954979  807.26886        NaN
 2   77.570076 -23.955711  804.86220        NaN
 3   77.569191 -23.956451  801.61490        NaN
 4   77.568314 -23.957195  797.30896        NaN
 5   77.567429 -23.957932  792.91420   2.570374
 6   77.566544 -23.958658  788.39480   1.922424
 7   77.565666 -23.959396  784.11694   2.658264
 8   77.564781 -23.960142  779.79500   3.664246
 9   77.563904 -23.960869  776.76560   4.249573
 10  77.563019 -23.961609  767.44324   6.668396
 11  77.562141 -23.962358  762.40380  10.477356
 12  77.544502 -23.977083  763.87840   8.553955
 13  77.543617 -23.977800  764.07556   7.998474
 14  77.542732 -23.978540  765.19090   6.761414
 15  77.541855 -23.979282  768.75446   4.206970
 16  77.540970 -23.980009  771.93854        NaN
 17  77.540085 -23.980747  773.26587        NaN
 18  77.539207 -23.981489  774.43463        NaN
 19  77.538322 -23.982227  776.89930        NaN
 20  77.537445 -23.982954  777.96080        NaN
 21  77.536560 -23.983694  779.56340        NaN
 22  77.535675 -23.984425  769.83990  11.356079
 23  77.534798 -23.985142  770.10900  11.810181
 24  77.533913 -23.985880  772.83750  11.275330
 25  77.533035 -23.986628  775.42786  12.219727,
 'atl10': Empty DataFrame
 Columns: [lat, lon, h]
 Index: [],
 'atl12': Empty DataFrame
 Columns: [lat, lon, h]
 Index: [],
 'atl13': Empty DataFrame
 Columns: [lat, lon, h]
 Index: []}

Plotting the ICESat-2 data#

Now let’s plot this data. Here, we are just creating an empty figure fig with axes ax.

# create the figure and axis
fig, ax = plt.subplots(figsize=[9,5])

# plot the data products
atl06, = ax.plot(is2data.atl06.lat, is2data.atl06.h, c='C0', linestyle='-', label='ATL06')
atl08, = ax.plot(is2data.atl08.lat, is2data.atl08.h, c='C2', linestyle=':', label='ATL08')
if np.sum(~np.isnan(is2data.atl08.canopy))>0:
    atl08canopy = ax.scatter(is2data.atl08.lat, is2data.atl08.h+is2data.atl08.canopy, s=2, c='C2', label='ATL08 canopy')

# add labels, title and legend
ax.set_xlabel('latitude')
ax.set_ylabel('elevation in meters')
ax.set_title('Some ICESat-2 data I found on OpenAltimetry!')
ax.legend(loc='upper left')

# add some text to provide info on what is plotted
info = 'ICESat-2 track {track:d}-{beam:s} on {date:s}\n({lon:.4f}E, {lat:.4f}N)'.format(track=is2data.track, 
                                                                                        beam=is2data.beam.upper(), 
                                                                                        date=is2data.date, 
                                                                                        lon=np.mean(is2data.lonlims), 
                                                                                        lat=np.mean(is2data.latlims))
infotext = ax.text(0.01, -0.08, info,
                   horizontalalignment='left', 
                   verticalalignment='top', 
                   transform=ax.transAxes,
                   fontsize=7,
                   bbox=dict(edgecolor=None, facecolor='white', alpha=0.9, linewidth=0))

# set the axis limits
ax.set_xlim((is2data.atl03.lat.min(), is2data.atl03.lat.max()))
ax.set_ylim((741, 818));

Let’s add the ATL03 photons to better understand what might be going on here.

atl03 = ax.scatter(is2data.atl03.lat, is2data.atl03.h, s=1, color='black', label='ATL03', zorder=-1)
ax.legend(loc='upper left')
fig.tight_layout()

Saving the plot to a file#

fig.savefig('my-plot.jpg', dpi=300)

To make plots easier to produce, the dataCollector class also has a method to plot the data that we downloaded.

fig = is2data.plotData();
fig

../../_images/OpenAltimetry_Earth_Engine_27_0.png

Ground Track Stats#

So far we have only seen the data in elevation vs. latitude space. It’s probably good to know what the scale on the x-axis is here in units that we’re familiar with.

def dist_latlon2meters(lat1, lon1, lat2, lon2):
    # returns the distance between two coordinate points - (lon1, lat1) and (lon2, lat2) along the earth's surface in meters.
    R = 6371000
    def deg2rad(deg):
        return deg * (np.pi/180)
    dlat = deg2rad(lat2-lat1)
    dlon = deg2rad(lon2-lon1)
    a = np.sin(dlat/2) * np.sin(dlat/2) + np.cos(deg2rad(lat1)) * np.cos(deg2rad(lat2)) * np.sin(dlon/2) * np.sin(dlon/2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
    return R * c

lat1, lat2 = is2data.atl08.lat[0], is2data.atl08.lat.iloc[-1]
lon1, lon2 = is2data.atl08.lon[0], is2data.atl08.lon.iloc[-1]

ground_track_length = dist_latlon2meters(lat1, lon1, lat2, lon2)
print('The ground track is about %.1f kilometers long.' % (ground_track_length/1e3))

The ground track is about 4.4 kilometers long.

Putting it all together#

The code above is found more concisely in one more method:

dataCollector.visualize_sentinel2()

This one has some default parameters, which you can change when running it:

- max_cloud_prob = 20
- days_buffer = 10
- gamma_value = 1.8
- title = 'ICESat-2 data'
- imagery_filename = 'my-satellite-image.tif'
- plot_filename = 'my-plot.jpg'

We can now do everything we did in this tutorial in just three lines!

url = 'http://openaltimetry.org/data/api/icesat2/atl03?date=2021-08-22&minx=-23.989069717620517&miny=77.53261687779506&maxx=-23.950348700496775&maxy=77.57222464044443&trackId=909&beamName=gt1r&beamName=gt1l&outputFormat=json'
is2data = dataCollector(oaurl=url,beam='gt1r')
fig = is2data.visualize_sentinel2(max_cloud_prob=5, title='some ICESat-2 data with a Sentinel-2 image')

--> Getting data from OpenAltimetry.
---> requesting ATL03 data... 6314 data points.
---> requesting ATL06 data... 198 data points.
---> requesting ATL07 data... No data.
---> requesting ATL08 data... 26 data points.
---> requesting ATL10 data... No data.
---> requesting ATL12 data... No data.
---> requesting ATL13 data... No data.
The ground track is 4.4 km long.
Looking for Sentinel-2 images from 2021-08-12T12:00:00 to 2021-09-01T12:00:00 --> there are 4 cloud-free images.
Looking for Sentinel-2 images from 2021-08-02T12:00:00 to 2021-09-11T12:00:00 --> there are 8 cloud-free images.
Looking for Sentinel-2 images from 2021-07-23T12:00:00 to 2021-09-21T12:00:00 --> there are 24 cloud-free images.
--> Closest cloud-free Sentinel-2 image to ICESat:
    - product_id: S2B_MSIL2A_20210820T144749_N0301_R082_T26XNM_20210820T173503
    - time difference: 2 days before ICESat-2
    - mean along-track cloud probability: 1.3
--> Downloaded the 8-bit RGB image as my-satellite-image.tif.
--> Saved plot as my-plot.jpg.

And now we can also easily run this on some totally different data!

url = 'http://openaltimetry.org/data/api/icesat2/atl12?date=2020-12-15&minx=-77.858681&miny=25.728091&maxx=-77.831461&maxy=25.832559&trackId=1254&beamName=gt1r&beamName=gt1l&outputFormat=json'
mydata = dataCollector(oaurl=url, beam='gt1r')
fig = mydata.visualize_sentinel2(max_cloud_prob=20,
                                 title='Nearshore Bathymetry ICESat-2 Data', 
                                 imagery_filename='my-other-satellite-image.tif',
                                 plot_filename='nearshore-bathymetry.jpg')

--> Getting data from OpenAltimetry.
---> requesting ATL03 data... 34237 data points.
---> requesting ATL06 data... No data.
---> requesting ATL07 data... No data.
---> requesting ATL08 data... 114 data points.
---> requesting ATL10 data... No data.
---> requesting ATL12 data... 2 data points.
---> requesting ATL13 data... 317 data points.
The ground track is 11.6 km long.
Looking for Sentinel-2 images from 2020-12-05T12:00:00 to 2020-12-25T12:00:00 --> there are 4 cloud-free images.
Looking for Sentinel-2 images from 2020-11-25T12:00:00 to 2021-01-04T12:00:00 --> there are 10 cloud-free images.
--> Closest cloud-free Sentinel-2 image to ICESat:
    - product_id: S2A_MSIL2A_20201213T155521_N0214_R011_T17RRJ_20201213T195958
    - time difference: 2 days before ICESat-2
    - mean along-track cloud probability: 5.7
--> Downloaded the 8-bit RGB image as my-other-satellite-image.tif.
--> Saved plot as nearshore-bathymetry.jpg.

Exercise#

Find some data from OpenAltimetry, and plot it with a cloud-free satellite image.

Look for small-scale features - say a few hundred meters to 20 kilometers along-track. (Hint: OpenAltimetry has a scale bar.)
Bonus points for features where concurrent imagery and the ATL03 photon cloud may give us some information that we would not get from higher-level products alone!
Note: There is no Sentinel-2 coverage at higher latitudes than ~84°, and there will be no concurrent imagery during polar night.

If you don’t know where to start with OpenAltimetry, you can look here.

##### YOUR CODE GOES HERE - uncomment and edit it! 
# url = 'http://???.org/??'
# gtx = 'gt??'
# is2data = dataCollector(oaurl=url, beam=gtx)
# fig = is2data.visualize_sentinel2(title='<your figure title goes here>', 
#                                   imagery_filename='your-satellite-image.tif',
#                                   plot_filename='your-plot.jpg')

Summary#

🎉 Congratulations! You’ve completed this tutorial and have seen how we can put ICESat-2 photon-level data into context using Google Earth Engine and the OpenAltimetry API.

You can explore a few more use cases for this code (and possible solutions for Exercise 2) in this notebook.

References#

To further explore the topics of this tutorial see the following detailed documentation:

Bonus Material: Code for Extracting Along-Track Surface Reflectance Values#

Based on javascript code in this example.

# a function for buffering ICESat-2 along-track depth measurement locations by a footprint radius
def bufferPoints(radius, bounds=False):
    def buffer(pt):
        pt = ee.Feature(pt)
        return pt.buffer(radius).bounds() if bounds else pt.buffer(radius)
    return buffer

# a function for extracting Sentinel-2 band data for ICESat-2 along-track depth measurement locations
def zonalStats(ic, fc, reducer=ee.Reducer.mean(), scale=None, crs=None, bands=None, bandsRename=None,
               imgProps=None, imgPropsRename=None, datetimeName='datetime', datetimeFormat='YYYY-MM-dd HH:mm:ss'):

    # Set default parameters based on an image representative.
    imgRep = ic.first()
    nonSystemImgProps = ee.Feature(None).copyProperties(imgRep).propertyNames()
    if not bands:
        bands = imgRep.bandNames()
    if not bandsRename:
        bandsRename = bands
    if not imgProps:
        imgProps = nonSystemImgProps
    if not imgPropsRename:
        imgPropsRename = imgProps

    # Map the reduceRegions function over the image collection.
    results = ic.map(lambda img: 
        img.select(bands, bandsRename)
        .set(datetimeName, img.date().format(datetimeFormat))
        .set('timestamp', img.get('system:time_start'))
        .reduceRegions(collection=fc.filterBounds(img.geometry()),reducer=reducer,scale=scale,crs=crs)
        .map(lambda f: f.set(img.toDictionary(imgProps).rename(imgProps,imgPropsRename)))
    ).flatten().filter(ee.Filter.notNull(bandsRename))

    return results

ground_track_buffer = 7.5 # radius in meters, for a 15-m diameter footprint
ground_track_coordinates = list(zip(is2data.gt.lon, is2data.gt.lat))
gtx_feature = ee.Geometry.LineString(coords=ground_track_coordinates, proj='EPSG:4326', geodesic=True)
aoi = gtx_feature.buffer(ground_track_buffer)

# clip to the ground track for along-track querying
img_gt = cloudfree_collection.first().clip(aoi)

# get the footprint-buffered query points from the lake object's depth_data dataframe
pts =  ee.FeatureCollection([
    ee.Feature(ee.Geometry.Point([r.lon, r.lat]), {'plot_id': i}) for i,r in is2data.gt.iterrows()
])
ptsS2 = pts.map(bufferPoints(ground_track_buffer))

# query the Sentinel-2 bands at the ground track locations where lake depths are posted
thiscoll = ee.ImageCollection([img_gt])
bandNames = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B11', 'B12', 'probability']
bandIdxs = list(range(len(bandNames)-1)) + [23]
ptsData = zonalStats(thiscoll, ptsS2, bands=bandIdxs, bandsRename=bandNames, imgProps=['PRODUCT_ID'], scale=5)

# get the data to the client side and merge into the lake depth data frame
features_select = ['plot_id'] + bandNames
results = ptsData.select(features_select, retainGeometry=False).getInfo()['features']
data = [x["properties"] for x in results]
dfS2 = pd.DataFrame(data).set_index('plot_id')
df_bands = is2data.gt.join(dfS2, how='left')

# calculate the Normalized Difference Water Index (NDWI)
df_bands['ndwi'] = (df_bands.B2-df_bands.B4) / (df_bands.B2+df_bands.B4)

df_bands

	index	lat	lon	xatc	B1	B11	B12	B2	B3	B4	B5	B6	B7	B8	B8A	B9	probability	ndwi
0	499	77.572104	-23.954022	0.000000	7542.079621	164.685412	112.618040	7679.957684	7304.628062	6419.345212	6360.185412	5687.488864	4959.208797	4554.002227	4227.100223	4271.315702	0.000000	0.089410
1	498	77.572025	-23.954088	8.906222	7519.000000	164.866406	112.482392	7673.486864	7283.530464	6385.437675	6358.749581	5691.365567	4957.435998	4589.643376	4232.547233	4267.000000	0.000000	0.091618
2	497	77.571946	-23.954154	17.812443	7519.000000	164.214807	110.955132	7667.717330	7348.897364	6504.787437	6403.295008	5750.389792	4981.712283	4688.150308	4231.484016	4267.000000	0.000000	0.082055
3	496	77.571867	-23.954220	26.718665	7519.000000	163.100390	110.193530	7614.340212	7219.821528	6464.990519	6421.064696	5779.469604	4995.268823	4696.385945	4231.175683	4267.000000	0.000000	0.081634
4	495	77.571789	-23.954286	35.624887	7519.000000	163.584213	104.806559	7571.846581	7152.758199	6384.584769	6316.176209	5657.226237	4909.948860	4619.942190	4142.191217	4267.000000	0.000000	0.085069
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
495	4	77.533083	-23.986586	4408.579718	8041.000000	201.959686	124.979843	8068.022396	7621.547592	6661.511758	6541.088466	5891.798432	5152.294513	4894.360582	4448.604703	4480.000000	0.386338	0.095489
496	3	77.533004	-23.986653	4417.485940	8041.000000	200.616331	124.308166	8119.082774	7650.724832	6645.447427	6577.359060	5885.081655	5128.785794	4903.615213	4468.755034	4480.000000	0.084452	0.099809
497	2	77.532925	-23.986717	4426.392162	8040.609131	199.819599	123.969933	8081.496659	7619.690423	6637.523385	6592.406459	5883.623608	5117.819599	4823.741648	4475.865256	4477.143653	0.152004	0.098103
498	1	77.532846	-23.986783	4435.298383	8039.849750	189.611018	127.637173	8041.181970	7589.139677	6624.518642	6584.066778	5910.378408	5111.065109	4767.759599	4448.170840	4471.594324	0.152476	0.096597
499	0	77.532768	-23.986849	4444.204605	8038.890128	180.460680	131.649191	8112.696040	7589.586168	6630.025655	6579.622978	5937.141104	5107.816509	4779.498048	4421.994980	4464.581707	0.485778	0.100570

500 rows × 18 columns

fig, axs = plt.subplots(figsize=[9,5], nrows=3, sharex=True)
ax = axs[0]
ax.scatter(is2data.atl03.lat, is2data.atl03.h, s=1, c='k', label='ATL03')
ax.set_ylim((741, 818))
ax.legend()
ax = axs[1]
ax.plot(df_bands.lat,df_bands.ndwi, c='b', label='NDWI')
ax.set_ylim((0, 1))
ax.legend()
ax = axs[2]
scale = 1e-4
ax.plot(df_bands.lat,df_bands.B4*scale, c='r', label='B4')
ax.plot(df_bands.lat,df_bands.B3*scale, c='g', label='B3')
ax.plot(df_bands.lat,df_bands.B2*scale, c='b', label='B2')
ax.set_ylim((0, 1))
ax.legend()
ax.set_xlim((is2data.gt.lat.min(), is2data.gt.lat.max()))
ax.set_xlabel('latitude')

Text(0.5, 0, 'latitude')

ICESat-2 Hackweek 2023

Interactive Visualizion with Open Altimetry & Google Earth Engine

Contents

Interactive Visualizion with Open Altimetry & Google Earth Engine#

Learning Objectives#

Computing environment#

Google Earth Engine Authentication and Initialization#

Downloading data from the OpenAltimetry API#

A short explanation of how I got the data:#

Getting the OpenAltimetry info into python#

Requesting the data from the OpenAltimetry API#

Plotting the ICESat-2 data#

Saving the plot to a file#

Ground Track Stats#

Google Earth Engine#

Show a ground track on a map#

Query for Sentinel-2 images#

Access an image collection#

Filter an image collection by location and time#

Get image collection info#

Add a Sentinel-2 image to the map#

Calculate the along-track cloud probability#

Filter cloudy images#

Sort the collection by time difference from the ICESat-2 overpass#

Show the final image and ground track on the map#

Download images from Earth Engine#

Open a GeoTIFF in rasterio#

Plot a GeoTIFF in Matplotlib#

transform the ground track into the image CRS#

Putting it all together#

Exercise#

Summary#

References#

Bonus Material: Code for Extracting Along-Track Surface Reflectance Values#