NEON Single Pixel¶

This notebook is an excercise in running ISOFIT in pixel mode (dynamically from python), in order to get into some configuation details.

Prerequisites:

• Run the neon.ipynb notebook succesfully - we'll use some of the outputs created from that intiial sample run

Running Per-Pixel Retrievals¶

That gave us a baseline run, but its not clear what might be improved. Let's load up some individual retrievals so we can understand the performance a bit better.

To get started with executing ISOFIT for iterable improvements, ISOFIT is comprised of three primary pieces:

1. The ForwardModel object
2. The IO object
3. The Inverse object

There is also an Isofit object that holds each of the other three, and is commonly used for more generic runs (it is what is called by apply_oe).

In [1]:

%matplotlib inline

import os
import logging
from pathlib import Path
from types import SimpleNamespace

import matplotlib.pyplot as plt
import numpy as np
from spectral.io import envi

from isofit.core.isofit import Isofit
from isofit.core.fileio import IO
from isofit.core.forward import ForwardModel
from isofit.data import env
from isofit.inversion.inverse import Inversion
from isofit.inversion.inverse_simple import invert_algebraic
from isofit.configs import configs
from isofit.core.geometry import Geometry

logging.getLogger().setLevel(logging.INFO)

%matplotlib inline import os import logging from pathlib import Path from types import SimpleNamespace import matplotlib.pyplot as plt import numpy as np from spectral.io import envi from isofit.core.isofit import Isofit from isofit.core.fileio import IO from isofit.core.forward import ForwardModel from isofit.data import env from isofit.inversion.inverse import Inversion from isofit.inversion.inverse_simple import invert_algebraic from isofit.configs import configs from isofit.core.geometry import Geometry logging.getLogger().setLevel(logging.INFO)

In [2]:

# Extract the image locations of each point of interest (POI)
# These are defined in the NEON report as pixel locations, so we round here to convert to indices
report = {}
report['173647'] = {     # Upp L Y  | Low R Y  | Upp L X | Low R X
    'WhiteTarp': np.round([2224.9626, 2230.9771, 316.0078, 324.9385,]).astype(int),
    'BlackTarp': np.round([2224.9626, 2231.0032, 328.0086, 333.9731,]).astype(int),
    'Veg'      : np.round([2245.0381, 2258.8103, 343.9006, 346.9423,]).astype(int),
    'RoadEW'   : np.round([2214.9905, 2216.9978, 348.9902, 373.0080,]).astype(int),
    'RoadNS'   : np.round([2205.9580, 2225.9612, 357.9536, 359.9608,]).astype(int)
}
report['174150'] = {     # Upp L Y | Low R Y | Upp L X  | Low R X
    'WhiteTarp': np.round([653.9626, 659.9771, 3143.0078, 3151.9385]).astype(int),
    'BlackTarp': np.round([653.9626, 660.0032, 3155.0086, 3160.9731]).astype(int),
    'Veg'      : np.round([674.0381, 687.8103, 3170.9006, 3173.9423]).astype(int),
    'RoadEW'   : np.round([643.9905, 645.9978, 3175.9902, 3200.0080]).astype(int),
    'RoadNS'   : np.round([634.9580, 654.9612, 3184.9536, 3186.9608]).astype(int)
}

# Which NEON date to process - change this to process a different date
neon_id  = list(report.keys())[0]
neon_str = f"NIS01_20210403_{neon_id}"

# Select the locations from the neon id -- roi == Regions of Interest
roi = report[neon_id]

# Extract the image locations of each point of interest (POI) # These are defined in the NEON report as pixel locations, so we round here to convert to indices report = {} report['173647'] = { # Upp L Y | Low R Y | Upp L X | Low R X 'WhiteTarp': np.round([2224.9626, 2230.9771, 316.0078, 324.9385,]).astype(int), 'BlackTarp': np.round([2224.9626, 2231.0032, 328.0086, 333.9731,]).astype(int), 'Veg' : np.round([2245.0381, 2258.8103, 343.9006, 346.9423,]).astype(int), 'RoadEW' : np.round([2214.9905, 2216.9978, 348.9902, 373.0080,]).astype(int), 'RoadNS' : np.round([2205.9580, 2225.9612, 357.9536, 359.9608,]).astype(int) } report['174150'] = { # Upp L Y | Low R Y | Upp L X | Low R X 'WhiteTarp': np.round([653.9626, 659.9771, 3143.0078, 3151.9385]).astype(int), 'BlackTarp': np.round([653.9626, 660.0032, 3155.0086, 3160.9731]).astype(int), 'Veg' : np.round([674.0381, 687.8103, 3170.9006, 3173.9423]).astype(int), 'RoadEW' : np.round([643.9905, 645.9978, 3175.9902, 3200.0080]).astype(int), 'RoadNS' : np.round([634.9580, 654.9612, 3184.9536, 3186.9608]).astype(int) } # Which NEON date to process - change this to process a different date neon_id = list(report.keys())[0] neon_str = f"NIS01_20210403_{neon_id}" # Select the locations from the neon id -- roi == Regions of Interest roi = report[neon_id]

In [3]:

# Below are the default values for the ISOFIT environment. Change these if your environment differs

env.load('~/.isofit/isofit.ini') # Ini file to load
env.changeSection('DEFAULT')     # Section of the ini to use
# env.changeBase('~./isofit')      # Base path for ISOFIT extras (data, examples, etc)
# env.changePath('srtmnet', '/path/to/sRTMnet_v120.h5') # Overwrite the path to sRTMnet - copy this line for other products such as sixs if in non-default locations

print('Using environment paths:')
for key, path in env.items():
    print(f"- {key} = {path}")

# Below are the default values for the ISOFIT environment. Change these if your environment differs env.load('~/.isofit/isofit.ini') # Ini file to load env.changeSection('DEFAULT') # Section of the ini to use # env.changeBase('~./isofit') # Base path for ISOFIT extras (data, examples, etc) # env.changePath('srtmnet', '/path/to/sRTMnet_v120.h5') # Overwrite the path to sRTMnet - copy this line for other products such as sixs if in non-default locations print('Using environment paths:') for key, path in env.items(): print(f"- {key} = {path}")

Using environment paths:
- data = /home/runner/work/isofit-tutorials/data
- examples = /home/runner/work/isofit-tutorials/isofit-tutorials
- imagecube = /home/runner/work/isofit-tutorials/imagecube
- srtmnet = /home/runner/work/isofit-tutorials/srtmnet
- sixs = /home/runner/work/isofit-tutorials/sixs
- surface = /home/runner/work/isofit-tutorials/surface
- plots = /home/runner/work/isofit-tutorials/plots
- libradtran = /home/runner/work/isofit-tutorials/libradtran
- srtmnet.file = sRTMnet_v120.h5
- srtmnet.aux = sRTMnet_v120_aux.npz
- libradtran.version = libRadtran-2.0.6

In [4]:

# Set the paths for this tutorial
base = Path(env.path('examples', 'NEON'))
raws = base / 'data'
data = base / 'neon_subset'

paths = SimpleNamespace(
    rdn = str(data / f'{neon_str}_rdn_ort'),
    loc = str(data / f'{neon_str}_loc_ort'),
    obs = str(data / f'{neon_str}_obs_ort'),
    insitu = raws,
    output = base / 'output',
    working = base / f'output/NIS01_20210403_{neon_id}',
    surface = str(base / 'output/surface.mat'),
    surface_config = env.path('examples', '20171108_Pasadena', 'configs', 'ang20171108t184227_surface.json')
)

# Set the paths for this tutorial base = Path(env.path('examples', 'NEON')) raws = base / 'data' data = base / 'neon_subset' paths = SimpleNamespace( rdn = str(data / f'{neon_str}_rdn_ort'), loc = str(data / f'{neon_str}_loc_ort'), obs = str(data / f'{neon_str}_obs_ort'), insitu = raws, output = base / 'output', working = base / f'output/NIS01_20210403_{neon_id}', surface = str(base / 'output/surface.mat'), surface_config = env.path('examples', '20171108_Pasadena', 'configs', 'ang20171108t184227_surface.json') )

In [5]:

# Now run actual retrievals
config = configs.create_new_config(paths.working / "config" / f"{neon_str}_isofit.json")
config.input.measured_radiance_file = paths.rdn
config.input.obs_file = paths.obs
config.input.loc_file = paths.loc

fm = ForwardModel(config)
io = IO(config, fm)
inv = Inversion(config, fm)

# Now run actual retrievals config = configs.create_new_config(paths.working / "config" / f"{neon_str}_isofit.json") config.input.measured_radiance_file = paths.rdn config.input.obs_file = paths.obs config.input.loc_file = paths.loc fm = ForwardModel(config) io = IO(config, fm) inv = Inversion(config, fm)

INFO:root:Loading config file: /home/runner/work/isofit-tutorials/isofit-tutorials/NEON/output/NIS01_20210403_173647/config/NIS01_20210403_173647_isofit.json

/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:3824: RuntimeWarning: Mean of empty slice
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)
INFO:/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/isofit/radiative_transfer/radiative_transfer_engine.py:Loading from wavelength_file: /home/runner/work/isofit-tutorials/isofit-tutorials/NEON/output/NIS01_20210403_173647/data/wavelengths.txt

INFO:/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/isofit/radiative_transfer/radiative_transfer_engine.py:Prebuilt LUT provided

INFO:isofit.radiative_transfer.luts:Loading LUT into memory

WARNING:isofit.radiative_transfer.luts:thermal_upwelling is fully NaN, leaving as-is

WARNING:isofit.radiative_transfer.luts:thermal_downwelling is fully NaN, leaving as-is

INFO:/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/isofit/radiative_transfer/radiative_transfer_engine.py:LUT grid loaded from file

Inversions¶

In [6]:

region = "Veg"

# Find the tightened bounding box of all the regions of interest
r = np.stack(list(report[neon_id].values())).T
Y = r[0].min() - 5, r[1].max() + 5
X = r[2].min() - 5, r[3].max() + 5

area = roi[region]
bounds = slice(*area[:2] - Y[0]), slice(*area[2:] - X[0])

read = lambda file: np.mean(envi.open(file).open_memmap(interleave="bip")[bounds].copy(), axis=(0, 1))
rdn = read(f"{config.input.measured_radiance_file}.hdr")
obs = read(f"{config.input.obs_file}.hdr")
loc = read(f"{config.input.loc_file}.hdr")

region = "Veg" # Find the tightened bounding box of all the regions of interest r = np.stack(list(report[neon_id].values())).T Y = r[0].min() - 5, r[1].max() + 5 X = r[2].min() - 5, r[3].max() + 5 area = roi[region] bounds = slice(*area[:2] - Y[0]), slice(*area[2:] - X[0]) read = lambda file: np.mean(envi.open(file).open_memmap(interleave="bip")[bounds].copy(), axis=(0, 1)) rdn = read(f"{config.input.measured_radiance_file}.hdr") obs = read(f"{config.input.obs_file}.hdr") loc = read(f"{config.input.loc_file}.hdr")

In [7]:

geom = Geometry(obs=obs, loc=loc)
states = inv.invert(rdn, geom)

x_surface, x_RT, x_instrument = fm.unpack(states[-1,:])
rfl_est, coeffs = invert_algebraic(
    fm.surface, 
    fm.RT, 
    fm.instrument, 
    x_surface, 
    x_RT, 
    x_instrument,
    rdn, 
    geom,
)

geom = Geometry(obs=obs, loc=loc) states = inv.invert(rdn, geom) x_surface, x_RT, x_instrument = fm.unpack(states[-1,:]) rfl_est, coeffs = invert_algebraic( fm.surface, fm.RT, fm.instrument, x_surface, x_RT, x_instrument, rdn, geom, )

WARNING:root:Earth sun distance not provided. Proceeding without might cause some inaccuracies down the line

/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:3824: RuntimeWarning: Mean of empty slice
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/runner/work/isofit-tutorials/isofit-tutorials/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)

Plotting¶

In [8]:

def closest_wl(mv):
    return np.argmin(np.abs(io.meas_wl-mv))

wl_nan = io.meas_wl.copy()
wl_nan[closest_wl(1360):closest_wl(1410)] = np.nan
wl_nan[closest_wl(1800):closest_wl(1970)] = np.nan

fig = plt.figure(figsize=(14,5))

for n in range(0, states.shape[0], 3):
    plt.plot(wl_nan, states[n, :-2] + 0.04*n,
        label = f"Step {n}" if n else "Algebraic Inversion at Initial Guess"
    )

plt.legend()

def closest_wl(mv): return np.argmin(np.abs(io.meas_wl-mv)) wl_nan = io.meas_wl.copy() wl_nan[closest_wl(1360):closest_wl(1410)] = np.nan wl_nan[closest_wl(1800):closest_wl(1970)] = np.nan fig = plt.figure(figsize=(14,5)) for n in range(0, states.shape[0], 3): plt.plot(wl_nan, states[n, :-2] + 0.04*n, label = f"Step {n}" if n else "Algebraic Inversion at Initial Guess" ) plt.legend()

Out[8]:

<matplotlib.legend.Legend at 0x7fe13ee68050>