Xarray Basics

This notebook demonstrates some of the basic xarray functionality.

Open an interactive notebook:

Sign into the SHIFT SMCE Daskhub and select an instance size in a different tab
Follow this link
Select the “notebook” kernel

Streamlined Xarray Interface to AVIRIS Datasets

We have an experimental capability to access the entire set of AVIRIS gridded mosaics as a single Xarray dataset (that can be subset by space, time, and wavelength). To open the dataset, code like the following (pay careful attention to the open_dataset arguments, position of brackets, quotes, etc.)

[1]:

import sys
sys.path.append('/efs/SHIFT-Python-Utilities/')
from shift_python_utilities.intake_shift import shift_catalog
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython import display
import warnings
import pandas as pd
warnings.filterwarnings("ignore")

# Read in the data using the shift python utilities library
cat = shift_catalog()
dat = cat.aviris_v1_gridded.read_chunked()

# Data can also be oppened using xarray
# ds = xr.open_dataset("reference://", engine="zarr", backend_kwargs={
#     "consolidated": False,
#     "storage_options": {"fo": "s3://dh-shift-curated/aviris/v1/gridded/zarr.json"}
# })
dat

[1]:

<xarray.Dataset>
Dimensions:      (time: 13, y: 12023, wavelength: 425, x: 13739)
Coordinates:
  * time         (time) datetime64[us] 2022-02-24 2022-02-28 ... 2022-05-29
  * wavelength   (wavelength) float32 377.2 382.2 387.2 ... 2.496e+03 2.501e+03
  * x            (x) float64 7.177e+05 7.177e+05 ... 7.864e+05 7.864e+05
  * y            (y) float64 3.866e+06 3.866e+06 ... 3.806e+06 3.806e+06
Data variables:
    reflectance  (time, y, wavelength, x) float32 dask.array<chunksize=(1, 1, 425, 13739), meta=np.ndarray>
    spatial_ref  int64 ...
Attributes: (12/13)
    description:               flight_products/20220224/box_mosaics/box_rfl_p...
    samples:                   13739
    lines:                     12023
    bands:                     425
    header offset:             0
    file type:                 ENVI Standard
    ...                        ...
    interleave:                bil
    byte order:                0
    map info:                  ['UTM', '1', '1', '717720', '3865865', '5', '5...
    coordinate system string:  ['PROJCS["WGS_1984_UTM_Zone_10N"', 'GEOGCS["GC...
    wavelength:                ['377.1956495', '382.20564950000005', '387.215...
    fwhm:                      ['5.57', '5.58', '5.58', '5.58', '5.5900000000...

Select a pixel at a time step and plot the reflectance

[4]:

dsub = dat.sel(x=750_000, y=3_830_000, time="2022-02-24", method="nearest")
dsub.reflectance.plot()

[4]:

[<matplotlib.lines.Line2D at 0x7eff8e4bd730>]

../_images/notebooks_xarray_basics_6_1.png

Extract a time series for a pixel

[3]:

dsub = dat.sel(x=750_000, y=3_830_000, method="nearest").reflectance

fig, ax = plt.subplots()
plt.ylabel("reflectance")
plt.xlabel("wavelength")

x = dsub[0].wavelength
line, = ax.plot(x, dsub[0])
ax.set_title(str(pd.Timestamp(dsub.time[0].item())))

def animate(i):
    line.set_ydata(dsub[i])  # update the data
    ax.set_title(str(pd.Timestamp(dsub.time[i].item())))
    return line,

ani = animation.FuncAnimation(fig, animate, repeat=True,
                                    frames=len(dsub), interval=1000)


plt.close()
video = ani.to_html5_video()
html = display.HTML(video)
display.display(html)

Select multiple adjacent pixels at a particular time

[4]:

dsub = dat.sel(x=750_000, time="2022-02-24", method="nearest").drop("x").sel(y=slice(3_830_020, 3_830_000)).reflectance

fig, ax = plt.subplots()
plt.ylabel("reflectance")
plt.xlabel("wavelength")

x = dsub[0].wavelength
line, = ax.plot(x, dsub[0])
ax.set_title(str(dsub.y[0].item()))

def animate(i):
    line.set_ydata(dsub[i])  # update the data
    ax.set_title("x:750000 y: " + str(dsub.y[i].item()))
    return line,

ani = animation.FuncAnimation(fig, animate, repeat=True,
                                    frames=len(dsub), interval=1000)

plt.close()
video = ani.to_html5_video()
html = display.HTML(video)
display.display(html)

Calculate a vegetation index through time for an area

[5]:

# Slice an area using x, y, and time
aoi = dat.sel(x=slice(730300,731000), y=slice(3819660,3819050), time="2022-03-08")

# Reshape to increase computational speed
aoi = aoi.stack(combined=('y', 'x'))

# Calculate NDVI
red = aoi.sel(wavelength=660, method="nearest").reflectance
nir = aoi.sel(wavelength=880, method="nearest").reflectance
ndvi = (nir - red) / (nir + red)

# Return to the orginal shape
ndvi = ndvi.unstack('combined')

[6]:

# Plot
plt.figure(figsize=(8, 6), dpi=80)
ndvi.plot.pcolormesh("x", "y", robust=True, add_colorbar=True);

../_images/notebooks_xarray_basics_13_0.png