Skip to article frontmatterSkip to article content

Lessons learned working with the NSIDC dataset.
Dataset: SnowEx 2021; Senator Beck Basin and Grand Mesa
Tutorial Author: Brent Wilder

Computing environment

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

from spectral import *
import numpy as np
import matplotlib.pyplot as plt
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 from spectral import *
      2 import numpy as np
      3 import matplotlib.pyplot as plt

ModuleNotFoundError: No module named 'spectral'

SnowEx21 Spectral Reflectance Dataset

The data were collected using an airborne imaging spectrometer, AVIRIS-NG can be downloaded from here, https://nsidc.org/data/snex21_ssr/versions/1.

  • Reflectance is provided at 5 nm spectral resolution with a range of 380-2500 nm

  • For this dataset, the pixel resolution is 4 m

  • Data span from 19 March 2021 to 29 April 2021, and were collected in two snow-covered environments in Colorado: Senator Beck Basin and Grand Mesa

  • Each file will have a “.img” and “.hdr”. You need to have both of these in the same directory to open data.

Downloading necessary terrain and illumination data

The NSIDC repository does not contain the terrain/illumination information.

However, you can obtain it for the matching flightline (by its timestamp) at the following URL, https://search.earthdata.nasa.gov/ ,

and searching for “AVIRIS-NG L1B Calibrated Radiance, Facility Instrument Collection, V1”

  • You only need to download the “obs_ort” files for the flight of interest. Please note these are different than “obs” files (ort means orthorectified).

  • In the Granule ID search, you can use wildcars “*” on either end of “obs_ort” to reduce your search.

  • You may also want to use this bounding box to reduce your search:

    • SW: 37.55725,-108.58887

    • NE: 39.78206,-106.16309

Using python package, spectral, to open data

# INSERT YOUR PATHS HERE
path_to_aviris = '/data/Albedo/AVIRIS/ang20210429t191025_rfl_v2z1'
path_to_aviris_hdr = '/data/Albedo/AVIRIS/ang20210429t191025_rfl_v2z1.hdr'
path_to_terrain = '/data/Albedo/AVIRIS/ang20210429t191025_rfl_v2z1_obs_ort'
path_to_terrain_hdr = '/data/Albedo/AVIRIS/ang20210429t191025_rfl_v2z1_obs_ort.hdr'
# Open a test image
aviris = envi.open(path_to_aviris_hdr)

# Save to an array in memory
rfl_array = aviris.open_memmap(writeable=True)

# print shape. You can see here we have 425 spectral bands for a grid of 1848x699 pixels
rfl_array.shape
# You can create an array of the bands centers like this
bands = np.array(aviris.bands.centers)
bands
# A simple data visalization by selecting random indices
i = 900
j = 300
pixel = rfl_array[i,j,:]

fig, ax = plt.subplots(1, 1, figsize=(10,5))
plt.rcParams.update({'font.size': 18})
ax.scatter(bands, pixel, color='blue', s=20)
ax.set_xlabel('Wavelength [nm]')
ax.set_ylabel('Reflectance')
plt.show()

Lastly, a very important note!

Please notice that convention for aspect follows π-\pi to π\pi.

# Terrain bands:
# 0 - Path length (m)
# 1 - To sensor azimuth
# 2 - To sensor zenith
# 3 - To sun azimuth
# 4 - To sun zenith
# 5 - Solar phase
# 6 - Slope
# 7 - Aspect
# 8 - cosine(i) (local solar illumination angle)
# 9 - UTC Time
# 10 - Earth-sun distance (AU)

# open envi object
terrain = envi.open(path_to_terrain_hdr)

# Save to an array in memory
terrain_array = terrain.open_memmap(writeable=True)

# Grab just aspect and flatten (remove nan)
aspects = terrain_array[:,:,7].flatten()
aspects = aspects[aspects>-9999]


# Plot a histogram to show aspect range
fig, ax = plt.subplots(1, 1, figsize=(10,5))
plt.rcParams.update({'font.size': 18})
ax.hist(aspects, color='black', bins=50)
ax.set_xlabel('Aspect [degrees]')
ax.set_ylabel('Count')
plt.show()