2. Accessing cloud-hosted ITS_LIVE data

2. Accessing cloud-hosted ITS_LIVE data#

This notebook will demonstrate how to access cloud-hosted Inter-mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) data from AWS S3 buckets. Here, you will find examples of how to successfully access cloud-hosted data as well as some common errors and issues you may run into along the way, what they mean, and how to resolve them.

Learning goals#

Concepts#

  • Query and access cloud-optimized dataset from cloud object storage

  • Create a vector data object representing the footprint of a raster dataset

  • Preliminary visualization of data extent

Techniques#

Other useful resources#

These are resources that contain additional examples and discussion of the content in this notebook and more.

Note

This tutorial was updated Fall 2023 to reflect changes to ITS_LIVE data urls and various software libraries

Software + Setup#

%load_ext watermark

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 get_ipython().run_line_magic('load_ext', 'watermark')

File ~/micromamba/envs/itslive_tutorial_env/lib/python3.13/site-packages/IPython/core/interactiveshell.py:2480, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2478     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2479 with self.builtin_trap:
-> 2480     result = fn(*args, **kwargs)
   2482 # The code below prevents the output from being displayed
   2483 # when using magics with decorator @output_can_be_silenced
   2484 # when the last Python token in the expression is a ';'.
   2485 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/micromamba/envs/itslive_tutorial_env/lib/python3.13/site-packages/IPython/core/magics/extension.py:33, in ExtensionMagics.load_ext(self, module_str)
     31 if not module_str:
     32     raise UsageError('Missing module name.')
---> 33 res = self.shell.extension_manager.load_extension(module_str)
     35 if res == 'already loaded':
     36     print("The %s extension is already loaded. To reload it, use:" % module_str)

File ~/micromamba/envs/itslive_tutorial_env/lib/python3.13/site-packages/IPython/core/extensions.py:62, in ExtensionManager.load_extension(self, module_str)
     55 """Load an IPython extension by its module name.
     56 
     57 Returns the string "already loaded" if the extension is already loaded,
     58 "no load function" if the module doesn't have a load_ipython_extension
     59 function, or None if it succeeded.
     60 """
     61 try:
---> 62     return self._load_extension(module_str)
     63 except ModuleNotFoundError:
     64     if module_str in BUILTINS_EXTS:

File ~/micromamba/envs/itslive_tutorial_env/lib/python3.13/site-packages/IPython/core/extensions.py:77, in ExtensionManager._load_extension(self, module_str)
     75 with self.shell.builtin_trap:
     76     if module_str not in sys.modules:
---> 77         mod = import_module(module_str)
     78     mod = sys.modules[module_str]
     79     if self._call_load_ipython_extension(mod):

File ~/micromamba/envs/itslive_tutorial_env/lib/python3.13/importlib/__init__.py:88, in import_module(name, package)
     86             break
     87         level += 1
---> 88 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1387, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1360, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1324, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'watermark'

%xmode minimal

Exception reporting mode: Minimal

import geopandas as gpd
import hvplot.pandas
import hvplot.xarray
import json
import numpy as np
import os
import rioxarray as rxr
import s3fs
from shapely.geometry import Point, Polygon
import xarray as xr

%config InlineBackend.figure_format='retina'

%watermark

UsageError: Line magic function `%watermark` not found.

%watermark --iversions

UsageError: Line magic function `%watermark` not found.

ITS_LIVE data cube catalog#

The ITS_LIVE project details a number of data access options on their website. Here, we will be accessing ITS_LIVE data in the form of zarr data cubes that are stored in s3 buckets hosted by Amazon Web Services (AWS).

Let’s begin by looking at the GeoJSON data cubes catalog. Click this link to download the file. This catalog contains spatial information and properties of ITS_LIVE data cubes as well as the URL used to access each cube. Let’s take a look at the entry for a single data cube and the information that it contains:

itslive_info

The top portion of the picture shows the spatial extent of the data cube in lat/lon units. Below that, we have properties such as the epsg code of the coordinate reference system, the spatial footprint in projected units and the url of the zarr object.

Let’s take a look at the url more in-depth:

itslive_url

From this link we can see that we are looking at its_live data located in an s3 bucket hosted by amazon AWS. We cans see that we’re looking in the data cube directory and what seems to be version 2. The next bit gives us information about the global location of the cube (N40E080). The actual file name ITS_LIVE_vel_EPSG32645_G0120_X250000_Y4750000.zarr tells us that we are looking at ice velocity data (its_live also has elevation data), in the CRS associated with EPSG 32645 (this code indicates UTM zone 45N). X250000_Y4750000 tells us more about the spatial footprint of the datacube within the UTM zone.

Accessing ITS_LIVE data from python#

We’ve found the url associated with the tile we want to access, let’s try to open the data cube using xarray:

url1 = 'http://its-live-data.s3.amazonaws.com/datacubes/v2/N30E090/ITS_LIVE_vel_EPSG32646_G0120_X750000_Y3350000.zarr'

dc1 = xr.open_dataset(url1)

As you can see, this doesn’t quite work. Passing the url to xr.open_dataset() without specifying a backend, xarray will expect a netcdf file. Because we’re trying to open a zarr file we need to add an additional argument to xr.open_dataset(), shown in the next code cell. You can find more information here. In the following cell, the argument chunks="auto" is passed, which introduces dask into our workflow.

dc1 = xr.open_dataset(url1, engine= 'zarr', chunks="auto")
                                   # storage_options = {'anon':True}) <-- as of Fall 2023 this no longer needed

dc1

<xarray.Dataset> Size: 771GB
Dimensions:                     (mid_date: 25243, y: 833, x: 833)
Coordinates:
  * mid_date                    (mid_date) datetime64[ns] 202kB 2022-06-07T04...
  * x                           (x) float64 7kB 7.001e+05 7.003e+05 ... 8e+05
  * y                           (y) float64 7kB 3.4e+06 3.4e+06 ... 3.3e+06
Data variables: (12/60)
    M11                         (mid_date, y, x) float32 70GB dask.array<chunksize=(25243, 30, 30), meta=np.ndarray>
    M11_dr_to_vr_factor         (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    M12                         (mid_date, y, x) float32 70GB dask.array<chunksize=(25243, 30, 30), meta=np.ndarray>
    M12_dr_to_vr_factor         (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    acquisition_date_img1       (mid_date) datetime64[ns] 202kB dask.array<chunksize=(25243,), meta=np.ndarray>
    acquisition_date_img2       (mid_date) datetime64[ns] 202kB dask.array<chunksize=(25243,), meta=np.ndarray>
    ...                          ...
    vy_error_modeled            (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    vy_error_slow               (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    vy_error_stationary         (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    vy_stable_shift             (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    vy_stable_shift_slow        (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
    vy_stable_shift_stationary  (mid_date) float32 101kB dask.array<chunksize=(25243,), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                CF-1.8
    GDAL_AREA_OR_POINT:         Area
    author:                     ITS_LIVE, a NASA MEaSUREs project (its-live.j...
    autoRIFT_parameter_file:    http://its-live-data.s3.amazonaws.com/autorif...
    datacube_software_version:  1.0
    date_created:               25-Sep-2023 22:00:23
    ...                         ...
    s3:                         s3://its-live-data/datacubes/v2/N30E090/ITS_L...
    skipped_granules:           s3://its-live-data/datacubes/v2/N30E090/ITS_L...
    time_standard_img1:         UTC
    time_standard_img2:         UTC
    title:                      ITS_LIVE datacube of image pair velocities
    url:                        https://its-live-data.s3.amazonaws.com/datacu...

This one worked! Let’s stop here and define a function that we can use for a quick view of this data.

def get_bounds_polygon(input_xr):
    
    xmin = input_xr.coords['x'].data.min()
    xmax = input_xr.coords['x'].data.max()

    ymin = input_xr.coords['y'].data.min()
    ymax = input_xr.coords['y'].data.max()

    pts_ls = [(xmin, ymin), (xmax, ymin),(xmax, ymax), (xmin, ymax), (xmin, ymin)]

    crs = f"epsg:{input_xr.mapping.spatial_epsg}"

    polygon_geom = Polygon(pts_ls)
    polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom]) 
    
    return polygon

Let’s also write a quick function for reading in s3 objects from http urls. This will come in handy when we’re trying to test multiple urls

def read_in_s3(http_url, chunks = 'auto'):
    #s3_url = http_url.replace('http','s3')   <-- as of Fall 2023, can pass http urls directly to xr.open_dataset()
    #s3_url = s3_url.replace('.s3.amazonaws.com','')

    datacube = xr.open_dataset(http_url, engine = 'zarr',
                                #storage_options={'anon':True},
                                chunks = chunks)

    return datacube

Now let’s take a look at the cube we’ve already read in:

bbox = get_bounds_polygon(dc1)

get_bounds_polygon() returns a geopandas.GeoDataFrame object in the same projection as the velocity data object (local UTM). Re-project to latitude/longitude to view the object more easily on a map:

bbox = bbox.to_crs('EPSG:4326')

bbox.hvplot()

DataError: Supplied data does not contain specified dimensions, the following dimensions were not found: ['x', 'y']

PandasInterface expects tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html

poly = bbox.hvplot(legend=True,alpha=0.3, tiles='ESRI', color='red', geo=True)

DataError: Supplied data does not contain specified dimensions, the following dimensions were not found: ['Longitude', 'Latitude']

PandasInterface expects tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html

poly

NameError: name 'poly' is not defined

Now we can see where this granule lies. We are using pandas.hvplot here, which is a great tool for interactive data visualization. You can read more about it here. Later in the tutorial we will demonstrate other situations where these tools are useful.

Searching ITS_LIVE catalog#

Let’s take a look at how we could search the ITS_LIVE data cube catalog for the data that we’re interested in. There are many ways to do this, this is just one example.

First, we will read in the catalog geojson file:

itslive_catalog = gpd.read_file('https://its-live-data.s3.amazonaws.com/datacubes/catalog_v02.json')

Here we’ll show two options for filtering the catalog:

  1. Selecting granules that contain a specific point, and

  2. Selecting all granules within a single UTM zone (specified by epsg code). This will let us take stock of the spatial coverage of data cubes located at working urls within a certain UTM zone.

You could easily tweak these functions (or write your own!) to select granules based on different properties. Play around with the itslive_catalog object to become more familiar with the data it contains and different options for indexing.

import s3fs

fs  = s3fs.S3FileSystem(anon=True)
fs

<s3fs.core.S3FileSystem at 0x7f9e2988a900>

Selecting granules by a single point#

def find_granule_by_point(input_point):
    '''returns url for the granule (zarr datacube) containing a specified point. point must be passed in epsg:4326
    '''
    catalog = gpd.read_file('https://its-live-data.s3.amazonaws.com/datacubes/catalog_v02.json')

    #make shapely point of input point
    p = gpd.GeoSeries([Point(input_point[0], input_point[1])],crs='EPSG:4326')
    #make gdf of point
    gdf = gdf = gpd.GeoDataFrame({'label': 'point', 
                                  'geometry':p})
    #find row of granule 
    granule = catalog.sjoin(gdf, how='inner')

    url = granule['zarr_url'].values[0]
    return url

url = find_granule_by_point([95.180191, 30.645973])

Great, this function returned a single url corresponding to the data cube covering the point we supplied

url

'http://its-live-data.s3.amazonaws.com/datacubes/v2-updated-october2024/N30E090/ITS_LIVE_vel_EPSG32646_G0120_X750000_Y3350000.zarr'

Let’s use the read_in_s3 function we defined to open the datacube as an xarray.Dataset

datacube = read_in_s3(url)

and then the get_bounds_polyon function to take a look at the footprint using hvplot(). Hvplot is an API for data exploration and visualization that you can read more about here.

bbox_dc = get_bounds_polygon(datacube)

poly =  bbox_dc.to_crs('EPSG:4326').hvplot(legend=True,alpha=0.3, tiles='ESRI', color='red', geo=True)

DataError: Supplied data does not contain specified dimensions, the following dimensions were not found: ['Longitude', 'Latitude']

PandasInterface expects tabular data, for more information on supported datatypes see http://holoviews.org/user_guide/Tabular_Datasets.html

poly

NameError: name 'poly' is not defined

Great, now we know how to access the its_live data cubes for a given region and at a specific point.

Let’s take a quick first look at this data cube. The next notebooks will give more examples of inspecting and working with this data. To make it easier for now, we will only plot a few hundred time steps instead of the full time series, and we will look at a spatial subset of the dataset. We use Xarray’s index-based selection methods for this.

datacube_sub = datacube.isel(mid_date = slice(0,50), x= slice(300,500), y=slice(400,600))
#datacube_sub = datacube.isel(mid_date = slice(0,50))

datacube_sub

<xarray.Dataset> Size: 89MB
Dimensions:                     (mid_date: 50, y: 200, x: 200)
Coordinates:
  * mid_date                    (mid_date) datetime64[ns] 400B 2022-06-07T04:...
  * x                           (x) float64 2kB 7.361e+05 7.363e+05 ... 7.6e+05
  * y                           (y) float64 2kB 3.352e+06 ... 3.328e+06
Data variables: (12/60)
    M11                         (mid_date, y, x) float32 8MB dask.array<chunksize=(50, 20, 20), meta=np.ndarray>
    M11_dr_to_vr_factor         (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    M12                         (mid_date, y, x) float32 8MB dask.array<chunksize=(50, 20, 20), meta=np.ndarray>
    M12_dr_to_vr_factor         (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    acquisition_date_img1       (mid_date) datetime64[ns] 400B dask.array<chunksize=(50,), meta=np.ndarray>
    acquisition_date_img2       (mid_date) datetime64[ns] 400B dask.array<chunksize=(50,), meta=np.ndarray>
    ...                          ...
    vy_error_modeled            (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    vy_error_slow               (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    vy_error_stationary         (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    vy_stable_shift             (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    vy_stable_shift_slow        (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    vy_stable_shift_stationary  (mid_date) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
Attributes: (12/19)
    Conventions:                CF-1.8
    GDAL_AREA_OR_POINT:         Area
    author:                     ITS_LIVE, a NASA MEaSUREs project (its-live.j...
    autoRIFT_parameter_file:    http://its-live-data.s3.amazonaws.com/autorif...
    datacube_software_version:  1.0
    date_created:               25-Sep-2023 22:00:23
    ...                         ...
    s3:                         s3://its-live-data/datacubes/v2/N30E090/ITS_L...
    skipped_granules:           s3://its-live-data/datacubes/v2/N30E090/ITS_L...
    time_standard_img1:         UTC
    time_standard_img2:         UTC
    title:                      ITS_LIVE datacube of image pair velocities
    url:                        https://its-live-data.s3.amazonaws.com/datacu...

Use the vmax argument to set a maximum value for the plot:

datacube_sub.v.mean(dim='mid_date').plot();
_images/10eff121cc5e093dc75a88890dcd29d4fcf8554f24ca6588ce74c0aeaba49696.png

The above plot shows the mean of the magnitude of velocity variable along the time dimension. Note that you could plot the entire spatial extent of the data but it would take quite a bit of time. The rest of the notebooks in this chapter will demonstrate how to subset this spatially-large data object to specific areas of interest.

Conclusion#

  • This notebook demonstrated how to query and access zarr data cubes from the ITS_LIVE dataset. We also looked at reading data into python using xarray, introducing dask to our workflows, and some preliminary data visualization.

  • The next notebook in the tutorial will take a look at the dataset in more detail and demonstrate ways to organize and manipulate data using xarray.