diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..610de1e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+imagery/
\ No newline at end of file
diff --git a/emit_nc.py b/emit_nc.py
new file mode 100644
index 0000000..1c4bbb9
--- /dev/null
+++ b/emit_nc.py
@@ -0,0 +1,579 @@
+"""
+Convert NASA EMIT L1B netCDF radiance + observation files to ISOFIT-ready ENVI format.
+
+Unlike NEON, EMIT already stores physically-scaled float32 values (no
+integer/scale-factor reconstruction is needed) and ships radiance and
+observation geometry as two separate granule files:
+
+- ``EMIT_L1B_RAD_*.nc`` — at-sensor radiance (``/radiance``) plus per-pixel
+  geolocation (``/location/lon``, ``lat``, ``elev``)
+- ``EMIT_L1B_OBS_*.nc`` — per-pixel observation geometry (``/obs``), with its
+  own (duplicate) copy of ``/location``
+
+Both files carry native, non-orthorectified pushbroom data on a
+(downtrack, crosstrack, bands) grid, along with a GLT (``glt_x``/``glt_y``)
+that maps that native grid onto a map-projected output grid. The GLT is only
+needed to orthorectify retrieval *outputs* after the fact — apply_oe itself
+runs on the native grid, so this module does not use it and does not write
+ENVI map info for the outputs it produces.
+
+This module writes the same three ENVI binary files that ISOFIT's apply_oe
+expects:
+
+- ``.rad``  — at-sensor radiance, Band-Interleaved-by-Line (BIL)
+- ``.loc``  — per-pixel location (longitude, latitude, elevation), BIL
+- ``.obs``  — per-pixel observation geometry (angles, path length, time), BIP
+
+Call graph:
+----------
+::
+
+    convert_emit_nc                         # main entry point
+    ├── convert_rad                         # copy & write radiance
+    │       ├── _decode                     # bytes/numpy scalar → str
+    │       └── _write_envi_header          # write .hdr text file
+    ├── convert_loc                         # write location cube
+    │       └── _write_envi_header
+    └── convert_obs                         # write geometry cube
+            ├── _decode
+            └── _write_envi_header
+"""
+
+from __future__ import annotations
+
+import re
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+# Band name constants exposed so callers can reference them without magic strings
+LOC_BAND_NAMES = ["longitude", "latitude", "elevation_m"]
+
+NO_DATA = -9999.0
+
+
+# ---------------------------------------------------------------------------
+# Private helpers
+# ---------------------------------------------------------------------------
+
+
+# --------------------------------
+# Called by: convert_obs
+# Calls:     _decode (recursive)
+# --------------------------------
+def _decode(x) -> str:
+    """Return a plain string from an H5 scalar/element that may be bytes or a numpy void."""
+    if isinstance(x, bytes):
+        return x.decode("utf-8", errors="ignore")
+    if hasattr(x, "shape") and x.shape == () and x.dtype.kind in ("S", "O"):
+        return _decode(x[()])
+    return str(x)
+
+
+# --------------------------------
+# Called by: convert_rad, convert_loc, convert_obs
+# Calls:     (none)
+# --------------------------------
+def _write_envi_header(
+    path: str,
+    *,
+    samples: int,
+    lines: int,
+    bands: int,
+    interleave: str,
+    no_data: float = NO_DATA,
+    wavelengths: np.ndarray | None = None,
+    fwhm: np.ndarray | None = None,
+    band_names: list[str] | None = None,
+    map_info: str | None = None,
+    coord_sys: str | None = None,
+    description: str = "EMIT netCDF → ENVI",
+) -> None:
+    """
+    Write a minimal ENVI ``.hdr`` file.
+
+    Parameters
+    ----------
+    path : str
+        Output path for the header (typically ``<data_file>.hdr``).
+    samples, lines, bands : int
+        Spatial and spectral dimensions.
+    interleave : str
+        ``"bil"``, ``"bip"``, or ``"bsq"``.
+    no_data : float
+        Fill value written to the ``data ignore value`` field.
+    wavelengths : np.ndarray, optional
+        Centre wavelengths in nm (written for ``.rad`` files).
+    fwhm : np.ndarray, optional
+        Full-width half-maximum values in nm.
+    band_names : list[str], optional
+        Per-band labels.
+    map_info : str, optional
+        ENVI map info string. Only meaningful for data on a regular map
+        grid — the native EMIT swath grid is not, so ``.rad``/``.loc``/
+        ``.obs`` outputs never set this; only orthorectified outputs do.
+    coord_sys : str, optional
+        Coordinate system (WKT) string, paired with ``map_info``.
+    description : str
+        Short description written into the header.
+    """
+    rows = [
+        "ENVI",
+        f"description = {{{description}}}",
+        f"samples = {samples}",
+        f"lines   = {lines}",
+        f"bands   = {bands}",
+        "header offset = 0",
+        "file type = ENVI Standard",
+        "data type = 4",  # float32
+        "byte order = 0",  # little-endian
+        f"interleave = {interleave.lower()}",
+        f"data ignore value = {no_data}",
+    ]
+    if wavelengths is not None:
+        rows.append("wavelength = {" + ", ".join(f"{w:.6f}" for w in wavelengths) + "}")
+    if fwhm is not None:
+        rows.append("fwhm = {" + ", ".join(f"{v:.6f}" for v in fwhm) + "}")
+    if band_names is not None:
+        rows.append("band names = {" + ", ".join(band_names) + "}")
+    if map_info is not None:
+        rows.append(f"map info = {{{map_info}}}")
+    if coord_sys is not None:
+        rows.append(f"coordinate system string = {{{coord_sys}}}")
+
+    Path(path).write_text("\n".join(rows) + "\n")
+
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+
+# --------------------------------
+# Called by: convert_emit_nc
+# Calls:     _write_envi_header
+# --------------------------------
+def convert_rad(
+    rad_nc_path: str, out_basename: str, chunk_rows: int = 64
+) -> tuple[str, str]:
+    """
+    Convert an EMIT L1B RAD netCDF file to an ENVI BIL ``.rad`` file.
+
+    EMIT radiance is already physically scaled (µW cm⁻² sr⁻¹ nm⁻¹), so this
+    is a straight copy/reshape rather than a reconstruction — unlike NEON,
+    there is no integer + scale-factor decoding step.
+
+    Parameters
+    ----------
+    rad_nc_path : str
+        Path to an ``EMIT_L1B_RAD_*.nc`` file.
+    out_basename : str
+        Output path without extension.  A ``.rad`` data file and a ``.hdr``
+        header are written next to each other.
+    chunk_rows : int, default 64
+        Number of downtrack rows processed at once. EMIT radiance cubes are
+        large (~2 GB); reduce this if memory is limited.
+
+    Returns
+    -------
+    rad_path : str
+    hdr_path : str
+    """
+    rad_path = out_basename + ".rad"
+    hdr_path = out_basename + ".rad.hdr"
+
+    with h5py.File(rad_nc_path, "r") as f:
+        ds = f["radiance"]  # (downtrack, crosstrack, bands)
+        n_lines, n_samples, n_bands = ds.shape
+
+        wavelengths = f["sensor_band_parameters"]["wavelengths"][:]
+        fwhm = f["sensor_band_parameters"]["fwhm"][:]
+
+        with open(rad_path, "wb") as out:
+            row = 0
+            while row < n_lines:
+                r1 = min(row + chunk_rows, n_lines)
+                R = ds[row:r1].astype(np.float32)
+
+                # BIL: each row on disk is (bands, samples) — transpose (samples, bands)
+                for i in range(R.shape[0]):
+                    out.write(R[i].T.tobytes(order="C"))
+                row = r1
+
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bil",
+        wavelengths=wavelengths,
+        fwhm=fwhm,
+        description="EMIT L1B RAD netCDF → ENVI radiance",
+    )
+    return rad_path, hdr_path
+
+
+# --------------------------------
+# Called by: convert_emit_nc
+# Calls:     _write_envi_header
+# --------------------------------
+def convert_loc(rad_nc_path: str, out_basename: str) -> tuple[str, str]:
+    """
+    Export the EMIT ``location`` group (longitude, latitude, elevation) to an
+    ENVI BIL ``.loc`` file.
+
+    The ``location`` group is present (identically) in both the RAD and OBS
+    granule files; either may be passed here.
+
+    Parameters
+    ----------
+    rad_nc_path : str
+        Path to an ``EMIT_L1B_RAD_*.nc`` (or ``EMIT_L1B_OBS_*.nc``) file.
+    out_basename : str
+        Output path without extension.
+
+    Returns
+    -------
+    loc_path : str
+    hdr_path : str
+    """
+    with h5py.File(rad_nc_path, "r") as f:
+        loc_grp = f["location"]
+        lon = loc_grp["lon"][:]
+        lat = loc_grp["lat"][:]
+        elev = loc_grp["elev"][:]
+
+    loc = np.stack([lon, lat, elev], axis=-1).astype(np.float32)  # (lines, samples, 3)
+    n_lines, n_samples, n_bands = loc.shape
+
+    loc_path = out_basename + ".loc"
+    with open(loc_path, "wb") as out:
+        for i in range(n_lines):
+            out.write(loc[i].T.tobytes(order="C"))  # BIL: (bands, samples) per row
+
+    hdr_path = loc_path + ".hdr"
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bil",
+        band_names=LOC_BAND_NAMES,
+        description="EMIT netCDF → ENVI location",
+    )
+    return loc_path, hdr_path
+
+
+# --------------------------------
+# Called by: convert_emit_nc
+# Calls:     _decode, _write_envi_header
+# --------------------------------
+def convert_obs(obs_nc_path: str, out_basename: str) -> tuple[str, str]:
+    """
+    Export the EMIT L1B OBS ``obs`` cube (observation geometry) to an ENVI
+    BIP ``.obs`` file.
+
+    EMIT's 11 observation bands (already documented by name inside the file,
+    at ``/sensor_band_parameters/observation_bands``) are: path length,
+    to-sensor azimuth/zenith, to-sun azimuth/zenith, solar phase, slope,
+    aspect, cosine(i), UTC time (already decimal hours), and Earth-sun
+    distance. Band names are read directly from the file rather than
+    hardcoded, so this stays correct across EMIT product versions.
+
+    ISOFIT expects observation geometry in BIP (Band-Interleaved-by-Pixel)
+    format, which is why this file uses a different interleave from ``.rad``
+    and ``.loc``.
+
+    Parameters
+    ----------
+    obs_nc_path : str
+        Path to an ``EMIT_L1B_OBS_*.nc`` file.
+    out_basename : str
+        Output path without extension.
+
+    Returns
+    -------
+    obs_path : str
+    hdr_path : str
+    """
+    with h5py.File(obs_nc_path, "r") as f:
+        obs = f["obs"][:].astype(np.float32)  # (downtrack, crosstrack, bands)
+        band_names = [_decode(b) for b in f["sensor_band_parameters"]["observation_bands"][:]]
+
+    n_lines, n_samples, n_bands = obs.shape
+
+    obs_path = out_basename + ".obs"
+    with open(obs_path, "wb") as out:
+        for i in range(n_lines):
+            # BIP: each row is (samples, bands) in C order — no transpose needed
+            out.write(obs[i].tobytes(order="C"))
+
+    hdr_path = obs_path + ".hdr"
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bip",
+        band_names=band_names,
+        description="EMIT L1B OBS netCDF → ENVI observation geometry",
+    )
+    return obs_path, hdr_path
+
+
+# --------------------------------
+# Called by: (user / external callers)
+# Calls:     convert_rad, convert_loc, convert_obs
+# --------------------------------
+def convert_emit_nc(
+    rad_nc_path: str, obs_nc_path: str, out_root: str, *, subdir: bool = True
+) -> dict[str, str]:
+    """
+    Convert an EMIT L1B RAD/OBS netCDF pair to the full ISOFIT input triplet
+    (.rad, .loc, .obs).
+
+    Parameters
+    ----------
+    rad_nc_path : str
+        Path to the ``EMIT_L1B_RAD_*.nc`` granule file.
+    obs_nc_path : str
+        Path to the matching ``EMIT_L1B_OBS_*.nc`` granule file.
+    out_root : str
+        Root output directory.  With ``subdir=True`` (default), files are
+        placed in ``<out_root>/<rad_stem>/``; with ``subdir=False`` directly
+        in ``<out_root>/``.
+    subdir : bool, default True
+        Whether to create a per-file subdirectory named after the RAD file's
+        stem.
+
+    Returns
+    -------
+    dict[str, str]
+        Mapping of ``{"rad", "rad_hdr", "loc", "loc_hdr", "obs", "obs_hdr"}``
+        to the written file paths.
+
+    Examples
+    --------
+    >>> paths = convert_emit_nc("EMIT_L1B_RAD_...nc", "EMIT_L1B_OBS_...nc", "./envi_output")
+    >>> paths["rad"]
+    './envi_output/EMIT_L1B_RAD_.../EMIT_L1B_RAD_....rad'
+    """
+    rad_nc_path = Path(rad_nc_path)
+    out_dir = Path(out_root) / rad_nc_path.stem if subdir else Path(out_root)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    base = str(out_dir / rad_nc_path.stem)
+
+    rad, rad_hdr = convert_rad(str(rad_nc_path), base)
+    loc, loc_hdr = convert_loc(str(rad_nc_path), base)
+    obs, obs_hdr = convert_obs(str(obs_nc_path), base)
+
+    return {
+        "rad": rad,
+        "rad_hdr": rad_hdr,
+        "loc": loc,
+        "loc_hdr": loc_hdr,
+        "obs": obs,
+        "obs_hdr": obs_hdr,
+    }
+
+
+# ---------------------------------------------------------------------------
+# Orthorectification
+#
+# apply_oe runs on the native (non-orthorectified) swath grid produced above
+# and its reflectance output stays on that same grid. Getting a map-projected
+# raster out of that output is a separate resampling step, done with the GLT
+# (glt_x/glt_y) carried in the RAD/OBS files' ``location`` group: for each
+# cell of the map-projected ("ortho") grid, the GLT gives the 1-based
+# (row, col) of the native-grid pixel to sample — 0 means no native pixel
+# maps there. This section reads an arbitrary ENVI file (e.g. ISOFIT's
+# reflectance output), remaps it through the GLT, and writes the result as a
+# new ENVI file with proper map info this time, since it now genuinely is on
+# a regular map grid.
+# ---------------------------------------------------------------------------
+
+_ENVI_DTYPES = {
+    1: np.uint8,
+    2: np.int16,
+    3: np.int32,
+    4: np.float32,
+    5: np.float64,
+    12: np.uint16,
+}
+
+
+# --------------------------------
+# Called by: _read_envi
+# Calls:     (none)
+# --------------------------------
+def _parse_envi_header(hdr_path: str) -> dict[str, str]:
+    """Parse an ENVI ``.hdr`` file into a ``{field: value}`` dict of raw strings."""
+    text = Path(hdr_path).read_text()
+    fields: dict[str, str] = {}
+    for m in re.finditer(r"([a-zA-Z_ ]+?)\s*=\s*(\{[^}]*\}|[^\n]+)", text, re.DOTALL):
+        key = m.group(1).strip().lower()
+        val = m.group(2).strip()
+        if val.startswith("{") and val.endswith("}"):
+            val = val[1:-1].strip()
+        fields[key] = val
+    return fields
+
+
+# --------------------------------
+# Called by: orthorectify
+# Calls:     _parse_envi_header
+# --------------------------------
+def _read_envi(hdr_path: str) -> np.ndarray:
+    """
+    Read an ENVI raster into a ``(lines, samples, bands)`` array.
+
+    Supports BIL, BIP, and BSQ interleave and the handful of ENVI data-type
+    codes ISOFIT actually writes (see ``_ENVI_DTYPES``).
+    """
+    fields = _parse_envi_header(hdr_path)
+    samples = int(fields["samples"])
+    lines = int(fields["lines"])
+    bands = int(fields["bands"])
+    interleave = fields["interleave"].lower()
+    dtype = _ENVI_DTYPES[int(fields["data type"])]
+    offset = int(fields.get("header offset", 0))
+
+    data_path = hdr_path[: -len(".hdr")] if hdr_path.endswith(".hdr") else hdr_path
+    raw = np.fromfile(data_path, dtype=dtype, offset=offset)
+
+    if interleave == "bsq":
+        return raw.reshape(bands, lines, samples).transpose(1, 2, 0)
+    if interleave == "bil":
+        return raw.reshape(lines, bands, samples).transpose(0, 2, 1)
+    if interleave == "bip":
+        return raw.reshape(lines, samples, bands)
+    raise ValueError(f"Unsupported interleave: {interleave!r}")
+
+
+# --------------------------------
+# Called by: orthorectify
+# Calls:     _decode
+# --------------------------------
+def read_glt(nc_path: str) -> dict:
+    """
+    Read the GLT (geometric lookup table) and map georeferencing from an
+    EMIT L1B RAD or OBS netCDF file's ``location`` group.
+
+    Either file works — both carry an identical ``location`` group.
+
+    Returns
+    -------
+    dict with keys ``glt_x``, ``glt_y`` (1-based native-pixel indices, 0 =
+    no data), ``geotransform`` (GDAL-style 6-element array), and
+    ``spatial_ref`` (WKT coordinate system string).
+    """
+    with h5py.File(nc_path, "r") as f:
+        glt_x = f["location"]["glt_x"][:]
+        glt_y = f["location"]["glt_y"][:]
+        geotransform = f.attrs["geotransform"][:]
+        spatial_ref = _decode(f.attrs["spatial_ref"])
+    return {
+        "glt_x": glt_x,
+        "glt_y": glt_y,
+        "geotransform": geotransform,
+        "spatial_ref": spatial_ref,
+    }
+
+
+# --------------------------------
+# Called by: orthorectify
+# Calls:     (none)
+# --------------------------------
+def apply_glt(
+    cube: np.ndarray, glt_x: np.ndarray, glt_y: np.ndarray, fill_value: float = NO_DATA
+) -> np.ndarray:
+    """
+    Resample a native-grid ``(lines, samples, bands)`` cube onto the ortho
+    grid defined by a GLT, via nearest-neighbor lookup (no interpolation —
+    each ortho cell just copies the one native pixel the GLT points at).
+
+    Parameters
+    ----------
+    cube : np.ndarray
+        Native-grid data, shape ``(lines, samples, bands)``.
+    glt_x, glt_y : np.ndarray
+        Ortho-grid arrays of 1-based native-pixel column/row indices, shape
+        ``(ortho_lines, ortho_samples)``. A value of 0 in either means no
+        native pixel maps to that ortho cell.
+    fill_value : float
+        Value written to ortho cells with no corresponding native pixel.
+
+    Returns
+    -------
+    np.ndarray
+        Shape ``(ortho_lines, ortho_samples, bands)``.
+    """
+    n_ortho_lines, n_ortho_samples = glt_x.shape
+    n_bands = cube.shape[-1]
+
+    valid = (glt_x > 0) & (glt_y > 0)
+    out = np.full((n_ortho_lines, n_ortho_samples, n_bands), fill_value, dtype=cube.dtype)
+    out[valid] = cube[glt_y[valid] - 1, glt_x[valid] - 1]
+    return out
+
+
+# --------------------------------
+# Called by: (user / external callers)
+# Calls:     _read_envi, read_glt, apply_glt, _write_envi_header
+# --------------------------------
+def orthorectify(in_envi_hdr_path: str, nc_path: str, out_basename: str) -> tuple[str, str]:
+    """
+    Orthorectify an ENVI file (e.g. apply_oe's native-grid reflectance
+    output) onto EMIT's map-projected ortho grid, using the GLT carried in
+    an EMIT L1B RAD or OBS netCDF file.
+
+    Note this can be memory-heavy for full scenes: the output array is
+    ``ortho_lines x ortho_samples x bands`` in memory at once (e.g. for a
+    typical EMIT scene with a 285-band reflectance cube, several GB).
+
+    Parameters
+    ----------
+    in_envi_hdr_path : str
+        Path to the ``.hdr`` file of the native-grid ENVI raster to
+        orthorectify (BIL, BIP, or BSQ; the matching data file is assumed to
+        sit alongside it with the ``.hdr`` suffix stripped).
+    nc_path : str
+        Path to the EMIT L1B RAD or OBS netCDF file carrying the GLT for
+        this scene.
+    out_basename : str
+        Output path without extension.  A ``.ortho`` data file (BIP) and a
+        ``.ortho.hdr`` header are written next to each other.
+
+    Returns
+    -------
+    ortho_path : str
+    hdr_path : str
+    """
+    cube = _read_envi(in_envi_hdr_path)
+    glt = read_glt(nc_path)
+
+    ortho = apply_glt(cube, glt["glt_x"], glt["glt_y"])
+    n_ortho_lines, n_ortho_samples, n_bands = ortho.shape
+
+    ortho_path = out_basename + ".ortho"
+    np.ascontiguousarray(ortho).tofile(ortho_path)  # (lines, samples, bands) C-order == BIP
+
+    gt = glt["geotransform"]
+    map_info = (
+        f"Geographic Lat/Lon, 1.0000, 1.0000, {gt[0]:.6f}, {gt[3]:.6f}, "
+        f"{abs(gt[1]):.8f}, {abs(gt[5]):.8f}, WGS-84, units=Degrees"
+    )
+
+    hdr_path = ortho_path + ".hdr"
+    _write_envi_header(
+        hdr_path,
+        samples=n_ortho_samples,
+        lines=n_ortho_lines,
+        bands=n_bands,
+        interleave="bip",
+        map_info=map_info,
+        coord_sys=glt["spatial_ref"],
+        description="EMIT GLT orthorectified output",
+    )
+    return ortho_path, hdr_path
diff --git a/neon_h5.py b/neon_h5.py
new file mode 100644
index 0000000..344edfc
--- /dev/null
+++ b/neon_h5.py
@@ -0,0 +1,523 @@
+"""
+Convert NEON AOP HDF5 radiance files to ISOFIT-ready ENVI format.
+
+NEON stores at-sensor radiance in a split-integer encoding to save disk space:
+
+    radiance = RadianceIntegerPart + RadianceDecimalPart / Scale_Factor
+
+This module reconstructs float32 radiance and writes three ENVI binary files
+that ISOFIT's apply_oe expects:
+
+- ``.rad``  — at-sensor radiance, Band-Interleaved-by-Line (BIL)
+- ``.loc``  — per-pixel location (easting, northing, elevation), BIL
+- ``.obs``  — per-pixel observation geometry (angles, path length, time), BIP
+
+Call graph:
+----------
+::
+
+    convert_neon_h5                         # main entry point
+    ├── convert_rad                         # reconstruct & write radiance
+    │       ├── _get_scale                  # find Scale_Factor in H5 attrs
+    │       │       └── _to_float           # safe attr → float cast
+    │       ├── _to_float                   # safe attr → float cast
+    │       ├── _decode                     # bytes/numpy scalar → str
+    │       └── _write_envi_header          # write .hdr text file
+    ├── convert_loc                         # write location cube
+    │       └── _write_envi_header
+    └── convert_obs                         # write geometry cube
+            ├── _to_decimal_hours           # normalize time band
+            └── _write_envi_header
+"""
+
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Union
+
+import h5py
+import numpy as np
+
+# Band name constants exposed so callers can reference them without magic strings
+LOC_BAND_NAMES = ["easting_utm", "northing_utm", "elevation_m"]
+
+OBS_BAND_NAMES = [
+    "path_length",
+    "sensor_azimuth",
+    "sensor_zenith",
+    "solar_azimuth",
+    "solar_zenith",
+    "toa_azimuth",
+    "toa_zenith",
+    "view_azimuth",
+    "cos_incidence",
+    "time_decimal_hours",
+]
+
+
+# ---------------------------------------------------------------------------
+# Private helpers
+# ---------------------------------------------------------------------------
+
+
+# --------------------------------
+# Called by: convert_rad
+# Calls:     _decode (recursive)
+# --------------------------------
+def _decode(x) -> str:
+    """Return a plain string from an H5 scalar that may be bytes or a numpy void."""
+    if isinstance(x, bytes):
+        return x.decode("utf-8", errors="ignore")
+    if hasattr(x, "shape") and x.shape == () and x.dtype.kind in ("S", "O"):
+        return _decode(x[()])
+    return str(x)
+
+
+# --------------------------------
+# Called by: convert_rad, _get_scale
+# Calls:     (none)
+# --------------------------------
+def _to_float(x, default=None) -> float | None:
+    """Safely cast an H5 attribute value to float, returning *default* on failure."""
+    try:
+        if x is None:
+            return default
+        if np.isscalar(x):
+            return float(x)
+        if getattr(x, "shape", ()) == ():
+            return float(x[()])
+        if isinstance(x, (bytes, str)):
+            return float(x)
+    except Exception:
+        pass
+    return default
+
+
+# --------------------------------
+# Called by: convert_rad
+# Calls:     _to_float
+# --------------------------------
+def _get_scale(dec_ds, group) -> Union[float, np.ndarray]:
+    """
+    Return the radiance scale factor from H5 attributes.
+
+    Searches both the decimal-part dataset and the parent group because the
+    attribute name is inconsistent across NEON file versions.  Falls back to
+    1.0 (no scaling) if nothing is found.
+
+    Parameters
+    ----------
+    dec_ds :
+        HDF5 dataset for ``RadianceDecimalPart``.
+    group :
+        HDF5 group containing the radiance datasets.
+
+    Returns
+    -------
+    float or np.ndarray
+        Scalar scale factor, or a 1-D per-band array.
+    """
+    for obj in (dec_ds, group):
+        if obj is None:
+            continue
+        for key in ("Scale_Factor", "Scale", "scale_factor"):
+            if key not in getattr(obj, "attrs", {}):
+                continue
+            val = obj.attrs[key]
+            if np.isscalar(val) or getattr(val, "shape", ()) == ():
+                s = _to_float(val)
+                if s:
+                    return s
+            else:
+                arr = np.asarray(val)
+                if arr.ndim == 1:
+                    return arr.astype(np.float32)
+    return 1.0
+
+
+# --------------------------------
+# Called by: convert_obs
+# Calls:     (none)
+# --------------------------------
+def _to_decimal_hours(arr: np.ndarray) -> np.ndarray:
+    """
+    Normalise a time array to decimal hours.
+
+    NEON OBS_Data band 10 (time) may arrive as decimal hours, seconds since
+    midnight, or HHMMSS numeric format depending on the file version.
+
+    Parameters
+    ----------
+    arr : np.ndarray
+        Raw time values from OBS_Data.
+
+    Returns
+    -------
+    np.ndarray
+        Time in decimal hours, float32.
+    """
+    x = arr.astype(np.float64)
+    finite = np.isfinite(x)
+    if not np.any(finite):
+        return arr.astype(np.float32)
+    m = np.nanmax(x[finite])
+    if m <= 24.5:
+        return x.astype(np.float32)  # already decimal hours
+    if m < 86_400:
+        return (x / 3600.0).astype(np.float32)  # seconds → hours
+    # HHMMSS numeric → hours
+    h = np.floor(x / 10_000.0)
+    mm = np.floor((x - h * 10_000.0) / 100.0)
+    ss = x - h * 10_000.0 - mm * 100.0
+    return (h + mm / 60.0 + ss / 3600.0).astype(np.float32)
+
+
+# --------------------------------
+# Called by: convert_rad, convert_loc, convert_obs
+# Calls:     (none)
+# --------------------------------
+def _write_envi_header(
+    path: str,
+    *,
+    samples: int,
+    lines: int,
+    bands: int,
+    interleave: str,
+    no_data: float = -9999.0,
+    wavelengths: np.ndarray | None = None,
+    fwhm: np.ndarray | None = None,
+    map_info: str | None = None,
+    coord_sys: str | None = None,
+    band_names: list[str] | None = None,
+    description: str = "NEON HDF5 → ENVI",
+) -> None:
+    """
+    Write a minimal ENVI ``.hdr`` file.
+
+    Parameters
+    ----------
+    path : str
+        Output path for the header (typically ``<data_file>.hdr``).
+    samples, lines, bands : int
+        Spatial and spectral dimensions.
+    interleave : str
+        ``"bil"`` or ``"bip"``.
+    no_data : float
+        Fill value written to the ``data ignore value`` field.
+    wavelengths : np.ndarray, optional
+        Centre wavelengths in nm (written for ``.rad`` files).
+    fwhm : np.ndarray, optional
+        Full-width half-maximum values in nm.
+    map_info : str, optional
+        ENVI map info string.
+    coord_sys : str, optional
+        Coordinate system string.
+    band_names : list[str], optional
+        Per-band labels.
+    description : str
+        Short description written into the header.
+    """
+    rows = [
+        "ENVI",
+        f"description = {{{description}}}",
+        f"samples = {samples}",
+        f"lines   = {lines}",
+        f"bands   = {bands}",
+        "header offset = 0",
+        "file type = ENVI Standard",
+        "data type = 4",  # float32
+        "byte order = 0",  # little-endian
+        f"interleave = {interleave.lower()}",
+        f"data ignore value = {no_data}",
+    ]
+    if wavelengths is not None:
+        rows.append("wavelength = {" + ", ".join(f"{w:.6f}" for w in wavelengths) + "}")
+    if fwhm is not None:
+        rows.append("fwhm = {" + ", ".join(f"{v:.6f}" for v in fwhm) + "}")
+    if map_info is not None:
+        rows.append(f"map info = {{{map_info}}}")
+    if coord_sys is not None:
+        rows.append(f"coordinate system string = {{{coord_sys}}}")
+    if band_names is not None:
+        rows.append("band names = {" + ", ".join(band_names) + "}")
+
+    Path(path).write_text("\n".join(rows) + "\n")
+
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+
+# --------------------------------
+# Called by: convert_neon_h5
+# Calls:     _get_scale, _to_float, _decode, _write_envi_header
+# --------------------------------
+def convert_rad(
+    h5_path: str, out_basename: str, chunk_rows: int = 64
+) -> tuple[str, str]:
+    """
+    Convert a NEON HDF5 radiance file to an ENVI BIL ``.rad`` file.
+
+    NEON stores radiance as two integer datasets to reduce file size.  The
+    physical radiance (W m⁻² sr⁻¹ nm⁻¹) is reconstructed as::
+
+        R = RadianceIntegerPart + RadianceDecimalPart / Scale_Factor
+
+    The output is written as float32, Band-Interleaved-by-Line (BIL), which
+    is the format expected by ISOFIT's ``apply_oe``.
+
+    Parameters
+    ----------
+    h5_path : str
+        Path to the NEON L1 radiance HDF5 file.
+    out_basename : str
+        Output path without extension.  A ``.rad`` data file and a ``.hdr``
+        header are written next to each other.
+    chunk_rows : int, default 64
+        Number of image rows processed at once.  Reduce if memory is limited.
+
+    Returns
+    -------
+    rad_path : str
+    hdr_path : str
+    """
+    rad_path = out_basename + ".rad"
+    hdr_path = out_basename + ".rad.hdr"
+
+    with h5py.File(h5_path, "r") as f:
+        site = next(iter(f.keys()))
+        rad_grp = f[site]["Radiance"]
+
+        ds_int = rad_grp["RadianceIntegerPart"]  # (lines, samples, bands)
+        ds_dec = rad_grp["RadianceDecimalPart"]
+        n_lines, n_samples, n_bands = ds_int.shape
+
+        scale = _get_scale(ds_dec, rad_grp)
+        if not (
+            np.isscalar(scale)
+            or (isinstance(scale, np.ndarray) and scale.shape == (n_bands,))
+        ):
+            raise RuntimeError(
+                f"Unexpected scale shape {getattr(scale, 'shape', type(scale))}; "
+                "expected scalar or 1-D per-band array."
+            )
+
+        int_nd = _to_float(ds_int.attrs.get("Data_Ignore_Value"))
+        dec_nd = _to_float(ds_dec.attrs.get("Data_Ignore_Value"))
+
+        meta = rad_grp["Metadata"]
+        wavelengths = meta["Spectral_Data"]["Wavelength"][:]
+        fwhm = meta["Spectral_Data"]["FWHM"][:]
+        map_info = _decode(meta["Coordinate_System"]["Map_Info"][()])
+        try:
+            coord_sys = _decode(
+                meta["Coordinate_System"]["Coordinate_System_String"][()]
+            )
+        except Exception:
+            coord_sys = None
+
+        with open(rad_path, "wb") as out:
+            row = 0
+            while row < n_lines:
+                r1 = min(row + chunk_rows, n_lines)
+                I = ds_int[row:r1].astype(np.float32)
+                D = ds_dec[row:r1].astype(np.float32)
+
+                nodata_mask = np.zeros(I.shape, dtype=bool)
+                if int_nd is not None:
+                    nodata_mask |= I == int_nd
+                if dec_nd is not None:
+                    nodata_mask |= D == dec_nd
+
+                if np.isscalar(scale):
+                    R = I + D / float(scale)
+                else:
+                    R = I + D / scale.reshape(1, 1, -1)
+
+                R = R.astype(np.float32)
+                R[nodata_mask] = -9999.0
+
+                # BIL: each row on disk is (bands, samples) — transpose (samples, bands)
+                for i in range(R.shape[0]):
+                    out.write(R[i].T.tobytes(order="C"))
+                row = r1
+
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bil",
+        wavelengths=wavelengths,
+        fwhm=fwhm,
+        map_info=map_info,
+        coord_sys=coord_sys,
+        description="NEON HDF5 → ENVI radiance",
+    )
+    return rad_path, hdr_path
+
+
+# --------------------------------
+# Called by: convert_neon_h5
+# Calls:     _write_envi_header
+# --------------------------------
+def convert_loc(h5_path: str, out_basename: str) -> tuple[str, str]:
+    """
+    Export NEON IGM_Data (easting, northing, elevation) to an ENVI BIL ``.loc`` file.
+
+    Parameters
+    ----------
+    h5_path : str
+        Path to the NEON L1 radiance HDF5 file.
+    out_basename : str
+        Output path without extension.
+
+    Returns
+    -------
+    loc_path : str
+    hdr_path : str
+    """
+    with h5py.File(h5_path, "r") as f:
+        site = next(iter(f.keys()))
+        loc = f[site]["Radiance"]["Metadata"]["Ancillary_Rasters"]["IGM_Data"][:]
+
+    loc = loc.astype(np.float32)
+    loc[loc == -9999] = -9999.0
+    n_lines, n_samples, n_bands = loc.shape
+
+    loc_path = out_basename + ".loc"
+    with open(loc_path, "wb") as out:
+        for i in range(n_lines):
+            out.write(loc[i].T.tobytes(order="C"))  # BIL: (bands, samples) per row
+
+    hdr_path = loc_path + ".hdr"
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bil",
+        band_names=LOC_BAND_NAMES,
+        description="NEON HDF5 → ENVI location (IGM)",
+    )
+    return loc_path, hdr_path
+
+
+# --------------------------------
+# Called by: convert_neon_h5
+# Calls:     _to_decimal_hours, _write_envi_header
+# --------------------------------
+def convert_obs(h5_path: str, out_basename: str) -> tuple[str, str]:
+    """
+    Export NEON OBS_Data (observation geometry) to an ENVI BIP ``.obs`` file.
+
+    The ten output bands are:
+
+    1. path_length — slant range in meters
+    2. sensor_azimuth — degrees clockwise from north
+    3. sensor_zenith — degrees from nadir
+    4. solar_azimuth — degrees clockwise from north
+    5. solar_zenith — degrees from zenith
+    6. toa_azimuth
+    7. toa_zenith
+    8. view_azimuth
+    9. cos_incidence — cosine of solar incidence angle on the surface
+    10. time_decimal_hours — UTC acquisition time in decimal hours
+
+    ISOFIT expects observation geometry in BIP (Band-Interleaved-by-Pixel)
+    format, which is why this file uses a different interleave from ``.rad``
+    and ``.loc``.
+
+    Parameters
+    ----------
+    h5_path : str
+        Path to the NEON L1 radiance HDF5 file.
+    out_basename : str
+        Output path without extension.
+
+    Returns
+    -------
+    obs_path : str
+    hdr_path : str
+    """
+    with h5py.File(h5_path, "r") as f:
+        site = next(iter(f.keys()))
+        obs = f[site]["Radiance"]["Metadata"]["Ancillary_Rasters"]["OBS_Data"][
+            :, :, :10
+        ]
+
+    obs = obs.astype(np.float32)
+    obs[:, :, 9] = _to_decimal_hours(obs[:, :, 9])
+    obs[obs == -9999] = -9999.0
+    n_lines, n_samples, n_bands = obs.shape
+
+    obs_path = out_basename + ".obs"
+    with open(obs_path, "wb") as out:
+        for i in range(n_lines):
+            # BIP: each row is (samples, bands) in C order — no transpose needed
+            out.write(obs[i].tobytes(order="C"))
+
+    hdr_path = obs_path + ".hdr"
+    _write_envi_header(
+        hdr_path,
+        samples=n_samples,
+        lines=n_lines,
+        bands=n_bands,
+        interleave="bip",
+        band_names=OBS_BAND_NAMES,
+        description="NEON HDF5 → ENVI observation geometry",
+    )
+    return obs_path, hdr_path
+
+
+# --------------------------------
+# Called by: (user / external callers)
+# Calls:     convert_rad, convert_loc, convert_obs
+# --------------------------------
+def convert_neon_h5(
+    h5_path: str, out_root: str, *, subdir: bool = True
+) -> dict[str, str]:
+    """
+    Convert a NEON HDF5 file to the full ISOFIT input triplet (.rad, .loc, .obs).
+
+    Parameters
+    ----------
+    h5_path : str
+        Path to the NEON L1 radiance HDF5 file.
+    out_root : str
+        Root output directory.  With ``subdir=True`` (default), files are
+        placed in ``<out_root>/<stem>/``; with ``subdir=False`` directly in
+        ``<out_root>/``.
+    subdir : bool, default True
+        Whether to create a per-file subdirectory named after the H5 stem.
+
+    Returns
+    -------
+    dict[str, str]
+        Mapping of ``{"rad", "rad_hdr", "loc", "loc_hdr", "obs", "obs_hdr"}``
+        to the written file paths.
+
+    Examples
+    --------
+    >>> paths = convert_neon_h5("flight.h5", "./envi_output")
+    >>> paths["rad"]
+    './envi_output/flight/flight.rad'
+    """
+    h5_path = Path(h5_path)
+    out_dir = Path(out_root) / h5_path.stem if subdir else Path(out_root)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    base = str(out_dir / h5_path.stem)
+
+    rad, rad_hdr = convert_rad(str(h5_path), base)
+    loc, loc_hdr = convert_loc(str(h5_path), base)
+    obs, obs_hdr = convert_obs(str(h5_path), base)
+
+    return {
+        "rad": rad,
+        "rad_hdr": rad_hdr,
+        "loc": loc,
+        "loc_hdr": loc_hdr,
+        "obs": obs,
+        "obs_hdr": obs_hdr,
+    }
diff --git a/tests/test_emit_nc.py b/tests/test_emit_nc.py
new file mode 100644
index 0000000..2efca77
--- /dev/null
+++ b/tests/test_emit_nc.py
@@ -0,0 +1,436 @@
+"""
+Tests for emit_nc:
+
+Synthetic EMIT L1B RAD and OBS netCDF (HDF5-backed) files are constructed in
+a temporary directory for each test so no real EMIT granule is required. The
+structure mirrors the actual EMIT L1B product layout (verified against a
+real EMIT_L1B_RAD/OBS granule pair).
+"""
+
+import importlib.util
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+
+# ---------------------------------------------------------------------------
+# Import the module directly by path (repo root is flat, no package install
+# required to run these tests).
+# ---------------------------------------------------------------------------
+_spec = importlib.util.spec_from_file_location(
+    "emit_nc",
+    Path(__file__).resolve().parents[1] / "emit_nc.py",
+)
+emit_nc = importlib.util.module_from_spec(_spec)
+_spec.loader.exec_module(emit_nc)
+
+convert_rad = emit_nc.convert_rad
+convert_loc = emit_nc.convert_loc
+convert_obs = emit_nc.convert_obs
+convert_emit_nc = emit_nc.convert_emit_nc
+LOC_BAND_NAMES = emit_nc.LOC_BAND_NAMES
+apply_glt = emit_nc.apply_glt
+read_glt = emit_nc.read_glt
+orthorectify = emit_nc.orthorectify
+_write_envi_header = emit_nc._write_envi_header
+_read_envi = emit_nc._read_envi
+
+# ---------------------------------------------------------------------------
+# Synthetic netCDF (HDF5) fixtures
+# ---------------------------------------------------------------------------
+
+N_LINES, N_SAMPLES, N_BANDS = 4, 5, 6
+NO_DATA = -9999.0
+
+WAVELENGTHS = np.array([400.0, 420.0, 440.0, 460.0, 480.0, 500.0], dtype=np.float32)
+FWHM = np.full(N_BANDS, 8.5, dtype=np.float32)
+
+OBS_BAND_NAMES = [
+    "Path length (sensor-to-ground in meters)",
+    "To-sensor azimuth (0 to 360 degrees CW from N)",
+    "To-sensor zenith (0 to 90 degrees from zenith)",
+    "To-sun azimuth (0 to 360 degrees CW from N)",
+    "To-sun zenith (0 to 90 degrees from zenith)",
+    "Solar phase",
+    "Slope",
+    "Aspect",
+    "Cosine(i)",
+    "UTC Time (decimal hours for mid-line pixels)",
+    "Earth-sun distance (AU)",
+]
+N_OBS_BANDS = len(OBS_BAND_NAMES)
+
+
+def _make_location_group(f: h5py.File, lon, lat, elev) -> None:
+    loc = f.create_group("location")
+    loc.create_dataset("lon", data=lon)
+    loc.create_dataset("lat", data=lat)
+    loc.create_dataset("elev", data=elev)
+    # glt_x/glt_y intentionally omitted: the converter does not use them.
+
+
+def _make_rad_nc(path: Path, lon, lat, elev) -> Path:
+    """Write a minimal synthetic EMIT L1B RAD netCDF-like H5 file."""
+    rng = np.random.default_rng(0)
+    radiance = rng.uniform(0, 50, size=(N_LINES, N_SAMPLES, N_BANDS)).astype(np.float32)
+
+    with h5py.File(path, "w") as f:
+        ds = f.create_dataset("radiance", data=radiance)
+        ds.attrs["_FillValue"] = np.float32(NO_DATA)
+        ds.attrs["units"] = b"uW/cm^2/SR/nm"
+
+        sbp = f.create_group("sensor_band_parameters")
+        sbp.create_dataset("wavelengths", data=WAVELENGTHS)
+        sbp.create_dataset("fwhm", data=FWHM)
+
+        _make_location_group(f, lon, lat, elev)
+
+    return path
+
+
+def _make_obs_nc(path: Path, lon, lat, elev) -> Path:
+    """Write a minimal synthetic EMIT L1B OBS netCDF-like H5 file."""
+    rng = np.random.default_rng(1)
+    obs = rng.uniform(0, 100, size=(N_LINES, N_SAMPLES, N_OBS_BANDS)).astype(np.float32)
+    obs[:, :, 9] = 16.25  # UTC time already in decimal hours
+
+    with h5py.File(path, "w") as f:
+        ds = f.create_dataset("obs", data=obs)
+        ds.attrs["_FillValue"] = np.float32(NO_DATA)
+
+        sbp = f.create_group("sensor_band_parameters")
+        sbp.create_dataset(
+            "observation_bands",
+            data=np.array([n.encode() for n in OBS_BAND_NAMES], dtype=object),
+        )
+
+        _make_location_group(f, lon, lat, elev)
+
+    return path
+
+
+@pytest.fixture
+def nc_pair(tmp_path):
+    rng = np.random.default_rng(2)
+    lon = rng.uniform(-120.5, -120.0, size=(N_LINES, N_SAMPLES)).astype(np.float64)
+    lat = rng.uniform(34.5, 35.0, size=(N_LINES, N_SAMPLES)).astype(np.float64)
+    elev = rng.uniform(100, 500, size=(N_LINES, N_SAMPLES)).astype(np.float64)
+
+    rad_path = _make_rad_nc(tmp_path / "EMIT_L1B_RAD_001_test.nc", lon, lat, elev)
+    obs_path = _make_obs_nc(tmp_path / "EMIT_L1B_OBS_001_test.nc", lon, lat, elev)
+    return rad_path, obs_path, lon, lat, elev
+
+
+# ---------------------------------------------------------------------------
+# convert_rad
+# ---------------------------------------------------------------------------
+
+
+def test_convert_rad_files_created(tmp_path, nc_pair):
+    rad_nc, _, *_ = nc_pair
+    base = str(tmp_path / "out")
+    rad_path, hdr_path = convert_rad(str(rad_nc), base)
+    assert Path(rad_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_rad_file_size(tmp_path, nc_pair):
+    """BIL float32: file size = lines x bands x samples x 4 bytes."""
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    rad_path, _ = convert_rad(str(rad_nc), base)
+    expected = N_LINES * N_BANDS * N_SAMPLES * 4
+    assert Path(rad_path).stat().st_size == expected
+
+
+def test_convert_rad_values_passthrough(tmp_path, nc_pair):
+    """EMIT radiance is already physical — output should equal the source array exactly."""
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    convert_rad(str(rad_nc), base)
+
+    raw = np.fromfile(base + ".rad", dtype=np.float32)
+    bil = raw.reshape(N_LINES, N_BANDS, N_SAMPLES).transpose(0, 2, 1)
+
+    with h5py.File(rad_nc, "r") as f:
+        expected = f["radiance"][:]
+
+    np.testing.assert_array_equal(bil, expected)
+
+
+def test_convert_rad_header_content(tmp_path, nc_pair):
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_rad(str(rad_nc), base)
+    hdr = Path(hdr_path).read_text()
+    assert "interleave = bil" in hdr
+    assert "data type = 4" in hdr
+    assert "400.000000" in hdr  # first wavelength
+
+
+def test_convert_rad_chunking_matches_single_shot(tmp_path, nc_pair):
+    """A small chunk_rows should produce identical output to the default."""
+    rad_nc, *_ = nc_pair
+    base_a = str(tmp_path / "a")
+    base_b = str(tmp_path / "b")
+    convert_rad(str(rad_nc), base_a, chunk_rows=64)
+    convert_rad(str(rad_nc), base_b, chunk_rows=1)
+
+    a = np.fromfile(base_a + ".rad", dtype=np.float32)
+    b = np.fromfile(base_b + ".rad", dtype=np.float32)
+    np.testing.assert_array_equal(a, b)
+
+
+# ---------------------------------------------------------------------------
+# convert_loc
+# ---------------------------------------------------------------------------
+
+
+def test_convert_loc_files_created(tmp_path, nc_pair):
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    loc_path, hdr_path = convert_loc(str(rad_nc), base)
+    assert Path(loc_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_loc_file_size(tmp_path, nc_pair):
+    """BIL float32: lines x 3 bands x samples x 4 bytes."""
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    loc_path, _ = convert_loc(str(rad_nc), base)
+    expected = N_LINES * 3 * N_SAMPLES * 4
+    assert Path(loc_path).stat().st_size == expected
+
+
+def test_convert_loc_band_names_in_header(tmp_path, nc_pair):
+    rad_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_loc(str(rad_nc), base)
+    hdr = Path(hdr_path).read_text()
+    for name in LOC_BAND_NAMES:
+        assert name in hdr
+
+
+def test_convert_loc_values(tmp_path, nc_pair):
+    """LOC band order should be [longitude, latitude, elevation], matching the source."""
+    rad_nc, _, lon, lat, elev = nc_pair
+    base = str(tmp_path / "out")
+    loc_path, _ = convert_loc(str(rad_nc), base)
+
+    raw = np.fromfile(loc_path, dtype=np.float32)
+    bil = raw.reshape(N_LINES, 3, N_SAMPLES).transpose(0, 2, 1)
+
+    expected = np.stack([lon, lat, elev], axis=-1).astype(np.float32)
+    np.testing.assert_allclose(bil, expected, rtol=1e-5)
+
+
+def test_convert_loc_reads_from_obs_file_too(tmp_path, nc_pair):
+    """The location group is duplicated in the OBS file; convert_loc should work from either."""
+    _, obs_nc, lon, lat, elev = nc_pair
+    base = str(tmp_path / "out")
+    loc_path, _ = convert_loc(str(obs_nc), base)
+
+    raw = np.fromfile(loc_path, dtype=np.float32)
+    bil = raw.reshape(N_LINES, 3, N_SAMPLES).transpose(0, 2, 1)
+    expected = np.stack([lon, lat, elev], axis=-1).astype(np.float32)
+    np.testing.assert_allclose(bil, expected, rtol=1e-5)
+
+
+# ---------------------------------------------------------------------------
+# convert_obs
+# ---------------------------------------------------------------------------
+
+
+def test_convert_obs_files_created(tmp_path, nc_pair):
+    _, obs_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    obs_path, hdr_path = convert_obs(str(obs_nc), base)
+    assert Path(obs_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_obs_file_size(tmp_path, nc_pair):
+    """BIP float32: lines x samples x 11 bands x 4 bytes."""
+    _, obs_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    obs_path, _ = convert_obs(str(obs_nc), base)
+    expected = N_LINES * N_SAMPLES * N_OBS_BANDS * 4
+    assert Path(obs_path).stat().st_size == expected
+
+
+def test_convert_obs_interleave_is_bip(tmp_path, nc_pair):
+    _, obs_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_obs(str(obs_nc), base)
+    assert "interleave = bip" in Path(hdr_path).read_text()
+
+
+def test_convert_obs_band_names_from_file(tmp_path, nc_pair):
+    """Band names should be read from the file's own observation_bands dataset."""
+    _, obs_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_obs(str(obs_nc), base)
+    hdr = Path(hdr_path).read_text()
+    for name in OBS_BAND_NAMES:
+        assert name in hdr
+
+
+def test_convert_obs_values_passthrough(tmp_path, nc_pair):
+    """EMIT OBS time band is already decimal hours — values pass through unchanged."""
+    _, obs_nc, *_ = nc_pair
+    base = str(tmp_path / "out")
+    obs_path, _ = convert_obs(str(obs_nc), base)
+
+    raw = np.fromfile(obs_path, dtype=np.float32).reshape(N_LINES, N_SAMPLES, N_OBS_BANDS)
+    time_band = raw[:, :, 9]
+    np.testing.assert_allclose(time_band, 16.25, rtol=1e-5)
+
+    with h5py.File(obs_nc, "r") as f:
+        expected = f["obs"][:]
+    np.testing.assert_array_equal(raw, expected)
+
+
+# ---------------------------------------------------------------------------
+# convert_emit_nc (end-to-end triplet)
+# ---------------------------------------------------------------------------
+
+
+def test_convert_emit_nc_all_files_created(tmp_path, nc_pair):
+    rad_nc, obs_nc, *_ = nc_pair
+    paths = convert_emit_nc(str(rad_nc), str(obs_nc), str(tmp_path))
+    for key in ("rad", "rad_hdr", "loc", "loc_hdr", "obs", "obs_hdr"):
+        assert key in paths
+        assert Path(paths[key]).exists(), f"Missing: {key}"
+
+
+def test_convert_emit_nc_subdir(tmp_path, nc_pair):
+    """With subdir=True (default) files should land in <out>/<rad_stem>/."""
+    rad_nc, obs_nc, *_ = nc_pair
+    paths = convert_emit_nc(str(rad_nc), str(obs_nc), str(tmp_path))
+    stem = Path(rad_nc).stem
+    assert Path(paths["rad"]).parent.name == stem
+
+
+def test_convert_emit_nc_no_subdir(tmp_path, nc_pair):
+    """With subdir=False files should land directly in out_root."""
+    rad_nc, obs_nc, *_ = nc_pair
+    paths = convert_emit_nc(str(rad_nc), str(obs_nc), str(tmp_path), subdir=False)
+    assert Path(paths["rad"]).parent == tmp_path
+
+
+# ---------------------------------------------------------------------------
+# apply_glt / read_glt / orthorectify
+# ---------------------------------------------------------------------------
+
+GEOTRANSFORM = np.array([-120.5, 0.0005, 0.0, 35.0, 0.0, -0.0005])
+SPATIAL_REF = 'GEOGCS["WGS 84",DATUM["WGS_1984"]]'
+
+
+def test_apply_glt_places_pixels_and_fills_nodata():
+    # 2x3 native cube, 1 band, distinctive values so placement is obvious.
+    cube = np.arange(6, dtype=np.float32).reshape(2, 3, 1)
+
+    # 2x2 ortho grid: three cells map to native pixels (1-based), one is nodata.
+    glt_x = np.array([[1, 3], [0, 2]], dtype=np.int32)
+    glt_y = np.array([[1, 2], [0, 1]], dtype=np.int32)
+
+    ortho = apply_glt(cube, glt_x, glt_y, fill_value=-9999.0)
+
+    assert ortho.shape == (2, 2, 1)
+    np.testing.assert_array_equal(ortho[0, 0], cube[0, 0])  # (y=1,x=1) -> native (0,0)
+    np.testing.assert_array_equal(ortho[0, 1], cube[1, 2])  # (y=2,x=3) -> native (1,2)
+    np.testing.assert_array_equal(ortho[1, 0], [-9999.0])   # glt 0 -> nodata
+    np.testing.assert_array_equal(ortho[1, 1], cube[0, 1])  # (y=1,x=2) -> native (0,1)
+
+
+def test_apply_glt_multiband_consistent():
+    cube = np.random.default_rng(3).uniform(0, 1, size=(4, 5, 3)).astype(np.float32)
+    glt_x = np.array([[2, 5], [0, 3]], dtype=np.int32)
+    glt_y = np.array([[1, 4], [0, 2]], dtype=np.int32)
+
+    ortho = apply_glt(cube, glt_x, glt_y)
+
+    np.testing.assert_array_equal(ortho[0, 0], cube[0, 1])
+    np.testing.assert_array_equal(ortho[0, 1], cube[3, 4])
+    np.testing.assert_array_equal(ortho[1, 1], cube[1, 2])
+
+
+def test_read_glt(tmp_path):
+    glt_x = np.array([[0, 1], [2, 3]], dtype=np.int32)
+    glt_y = np.array([[0, 1], [1, 2]], dtype=np.int32)
+
+    nc_path = tmp_path / "glt_test.nc"
+    with h5py.File(nc_path, "w") as f:
+        loc = f.create_group("location")
+        loc.create_dataset("glt_x", data=glt_x)
+        loc.create_dataset("glt_y", data=glt_y)
+        f.attrs["geotransform"] = GEOTRANSFORM
+        f.attrs["spatial_ref"] = SPATIAL_REF.encode()
+
+    glt = read_glt(str(nc_path))
+    np.testing.assert_array_equal(glt["glt_x"], glt_x)
+    np.testing.assert_array_equal(glt["glt_y"], glt_y)
+    np.testing.assert_allclose(glt["geotransform"], GEOTRANSFORM)
+    assert glt["spatial_ref"] == SPATIAL_REF
+
+
+@pytest.mark.parametrize("interleave", ["bil", "bip", "bsq"])
+def test_read_envi_roundtrip(tmp_path, interleave):
+    lines, samples, bands = 3, 4, 2
+    cube = np.random.default_rng(4).uniform(0, 100, size=(lines, samples, bands)).astype(
+        np.float32
+    )
+
+    data_path = tmp_path / f"test.{interleave}"
+    hdr_path = str(data_path) + ".hdr"
+    with open(data_path, "wb") as out:
+        if interleave == "bsq":
+            out.write(np.ascontiguousarray(cube.transpose(2, 0, 1)).tobytes())
+        elif interleave == "bil":
+            for i in range(lines):
+                out.write(cube[i].T.tobytes(order="C"))
+        else:  # bip
+            out.write(np.ascontiguousarray(cube).tobytes())
+
+    _write_envi_header(
+        hdr_path, samples=samples, lines=lines, bands=bands, interleave=interleave
+    )
+
+    result = _read_envi(hdr_path)
+    np.testing.assert_array_equal(result, cube)
+
+
+def test_orthorectify_end_to_end(tmp_path):
+    # A tiny "isofit output"-style native-grid ENVI file, 2 bands, BIP.
+    native = np.arange(2 * 3 * 2, dtype=np.float32).reshape(2, 3, 2)
+    native_path = tmp_path / "reflectance.bip"
+    native_hdr = str(native_path) + ".hdr"
+    with open(native_path, "wb") as out:
+        out.write(np.ascontiguousarray(native).tobytes())
+    _write_envi_header(native_hdr, samples=3, lines=2, bands=2, interleave="bip")
+
+    glt_x = np.array([[1, 2], [0, 3]], dtype=np.int32)
+    glt_y = np.array([[1, 1], [0, 2]], dtype=np.int32)
+    nc_path = tmp_path / "glt.nc"
+    with h5py.File(nc_path, "w") as f:
+        loc = f.create_group("location")
+        loc.create_dataset("glt_x", data=glt_x)
+        loc.create_dataset("glt_y", data=glt_y)
+        f.attrs["geotransform"] = GEOTRANSFORM
+        f.attrs["spatial_ref"] = SPATIAL_REF.encode()
+
+    out_base = str(tmp_path / "ortho_out")
+    ortho_path, hdr_path = orthorectify(native_hdr, str(nc_path), out_base)
+
+    raw = np.fromfile(ortho_path, dtype=np.float32).reshape(2, 2, 2)
+    np.testing.assert_array_equal(raw[0, 0], native[0, 0])
+    np.testing.assert_array_equal(raw[0, 1], native[0, 1])
+    np.testing.assert_array_equal(raw[1, 0], [-9999.0, -9999.0])
+    np.testing.assert_array_equal(raw[1, 1], native[1, 2])
+
+    hdr = Path(hdr_path).read_text()
+    assert "map info" in hdr
+    assert "-120.500000" in hdr
+    assert SPATIAL_REF in hdr
diff --git a/tests/test_neon_h5.py b/tests/test_neon_h5.py
new file mode 100644
index 0000000..47e0d08
--- /dev/null
+++ b/tests/test_neon_h5.py
@@ -0,0 +1,290 @@
+"""
+Tests for isofit.utils.neon_h5:
+
+A synthetic NEON HDF5 file is constructed in a temporary directory for each
+test so no real data file is required.  The structure mirrors the actual NEON
+L1 radiance product layout.
+"""
+
+import importlib.util
+import struct
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+
+# ---------------------------------------------------------------------------
+# Import the module directly to avoid pulling in the full isofit package,
+# which requires optional heavy dependencies (torch, etc.).
+# ---------------------------------------------------------------------------
+_spec = importlib.util.spec_from_file_location(
+    "neon_h5",
+    Path(__file__).resolve().parents[1] / "utils" / "neon_h5.py",
+)
+neon_h5 = importlib.util.module_from_spec(_spec)
+_spec.loader.exec_module(neon_h5)
+
+convert_rad = neon_h5.convert_rad
+convert_loc = neon_h5.convert_loc
+convert_obs = neon_h5.convert_obs
+convert_neon_h5 = neon_h5.convert_neon_h5
+LOC_BAND_NAMES = neon_h5.LOC_BAND_NAMES
+OBS_BAND_NAMES = neon_h5.OBS_BAND_NAMES
+_to_decimal_hours = neon_h5._to_decimal_hours
+_get_scale = neon_h5._get_scale
+
+# ---------------------------------------------------------------------------
+# Synthetic H5 fixture
+# ---------------------------------------------------------------------------
+
+SITE = "NEON_TEST"
+N_LINES, N_SAMPLES, N_BANDS = 4, 5, 6
+SCALE = 10_000.0
+
+WAVELENGTHS = np.array([400.0, 420.0, 440.0, 460.0, 480.0, 500.0], dtype=np.float32)
+FWHM = np.full(N_BANDS, 10.0, dtype=np.float32)
+MAP_INFO = "UTM, 1, 1, 480000.0, 4442000.0, 1.0, 1.0, 13, North, WGS-84, units=Meters"
+COORD_SYS = "PROJCS[WGS 84 / UTM zone 13N]"
+
+
+def _make_h5(path: Path) -> Path:
+    """Write a minimal synthetic NEON H5 file to *path* and return it."""
+    rng = np.random.default_rng(0)
+
+    int_part = rng.integers(50, 200, size=(N_LINES, N_SAMPLES, N_BANDS), dtype=np.int16)
+    dec_part = rng.integers(0, 9999, size=(N_LINES, N_SAMPLES, N_BANDS), dtype=np.int16)
+
+    igm = rng.uniform(
+        [480000, 4_442_000, 1600],
+        [480100, 4_442_100, 1700],
+        size=(N_LINES, N_SAMPLES, 3),
+    ).astype(np.float32)
+
+    obs = rng.uniform(0, 1, size=(N_LINES, N_SAMPLES, 10)).astype(np.float32)
+    obs[:, :, 9] = 58_500.0  # seconds since midnight → 16.25 decimal hours
+
+    with h5py.File(path, "w") as f:
+        rad = f.create_group(f"{SITE}/Radiance")
+
+        ds_int = rad.create_dataset("RadianceIntegerPart", data=int_part)
+        ds_int.attrs["Data_Ignore_Value"] = -9999
+
+        ds_dec = rad.create_dataset("RadianceDecimalPart", data=dec_part)
+        ds_dec.attrs["Data_Ignore_Value"] = -9999
+        ds_dec.attrs["Scale_Factor"] = SCALE
+
+        spec = rad.create_group("Metadata/Spectral_Data")
+        spec.create_dataset("Wavelength", data=WAVELENGTHS)
+        spec.create_dataset("FWHM", data=FWHM)
+
+        crs = rad.create_group("Metadata/Coordinate_System")
+        crs.create_dataset("Map_Info", data=MAP_INFO.encode())
+        crs.create_dataset("Coordinate_System_String", data=COORD_SYS.encode())
+
+        anc = rad.create_group("Metadata/Ancillary_Rasters")
+        anc.create_dataset("IGM_Data", data=igm)
+        anc.create_dataset("OBS_Data", data=obs)
+
+    return path
+
+
+@pytest.fixture
+def h5_file(tmp_path):
+    return _make_h5(tmp_path / "test_flight.h5")
+
+
+# ---------------------------------------------------------------------------
+# Helper tests
+# ---------------------------------------------------------------------------
+
+
+def test_to_decimal_hours_passthrough():
+    """Values already in decimal hours should pass through unchanged."""
+    arr = np.array([16.25, 0.0, 23.99], dtype=np.float32)
+    result = _to_decimal_hours(arr)
+    np.testing.assert_allclose(result, arr, rtol=1e-5)
+
+
+def test_to_decimal_hours_from_seconds():
+    """Seconds since midnight should convert to decimal hours."""
+    arr = np.array([58_500.0], dtype=np.float32)  # 16.25 h × 3600
+    result = _to_decimal_hours(arr)
+    np.testing.assert_allclose(result, [16.25], rtol=1e-4)
+
+
+def test_to_decimal_hours_from_hhmmss():
+    """HHMMSS numeric format (e.g. 162500 = 16:25:00) should convert correctly."""
+    arr = np.array([162_500.0], dtype=np.float64)
+    result = _to_decimal_hours(arr)
+    np.testing.assert_allclose(result, [16.0 + 25.0 / 60.0], rtol=1e-4)
+
+
+def test_get_scale_reads_attribute(h5_file):
+    """_get_scale should find Scale_Factor on the decimal-part dataset."""
+    with h5py.File(h5_file, "r") as f:
+        rad = f[SITE]["Radiance"]
+        scale = _get_scale(rad["RadianceDecimalPart"], rad)
+    assert scale == SCALE
+
+
+# ---------------------------------------------------------------------------
+# convert_rad
+# ---------------------------------------------------------------------------
+
+
+def test_convert_rad_files_created(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    rad_path, hdr_path = convert_rad(str(h5_file), base)
+    assert Path(rad_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_rad_file_size(tmp_path, h5_file):
+    """BIL float32: file size = lines × bands × samples × 4 bytes."""
+    base = str(tmp_path / "out")
+    rad_path, _ = convert_rad(str(h5_file), base)
+    expected = N_LINES * N_BANDS * N_SAMPLES * 4
+    assert Path(rad_path).stat().st_size == expected
+
+
+def test_convert_rad_values(tmp_path, h5_file):
+    """Reconstructed radiance should equal integer + decimal / scale."""
+    base = str(tmp_path / "out")
+    convert_rad(str(h5_file), base)
+
+    # Read back the BIL binary manually
+    raw = np.fromfile(base + ".rad", dtype=np.float32)
+    # BIL layout: (lines, bands, samples) → reshape and transpose to (lines, samples, bands)
+    bil = raw.reshape(N_LINES, N_BANDS, N_SAMPLES).transpose(0, 2, 1)
+
+    with h5py.File(h5_file, "r") as f:
+        rad = f[SITE]["Radiance"]
+        int_part = rad["RadianceIntegerPart"][:].astype(np.float32)
+        dec_part = rad["RadianceDecimalPart"][:].astype(np.float32)
+
+    expected = int_part + dec_part / SCALE
+    np.testing.assert_allclose(bil, expected, rtol=1e-5)
+
+
+def test_convert_rad_header_content(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_rad(str(h5_file), base)
+    hdr = Path(hdr_path).read_text()
+    assert "interleave = bil" in hdr
+    assert "data type = 4" in hdr
+    assert "400.000000" in hdr  # first wavelength
+
+
+# ---------------------------------------------------------------------------
+# convert_loc
+# ---------------------------------------------------------------------------
+
+
+def test_convert_loc_files_created(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    loc_path, hdr_path = convert_loc(str(h5_file), base)
+    assert Path(loc_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_loc_file_size(tmp_path, h5_file):
+    """BIL float32: lines × 3 bands × samples × 4 bytes."""
+    base = str(tmp_path / "out")
+    loc_path, _ = convert_loc(str(h5_file), base)
+    expected = N_LINES * 3 * N_SAMPLES * 4
+    assert Path(loc_path).stat().st_size == expected
+
+
+def test_convert_loc_band_names_in_header(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_loc(str(h5_file), base)
+    hdr = Path(hdr_path).read_text()
+    for name in LOC_BAND_NAMES:
+        assert name in hdr
+
+
+def test_convert_loc_values(tmp_path, h5_file):
+    """LOC values should match the IGM_Data stored in the H5 file."""
+    base = str(tmp_path / "out")
+    loc_path, _ = convert_loc(str(h5_file), base)
+
+    raw = np.fromfile(loc_path, dtype=np.float32)
+    bil = raw.reshape(N_LINES, 3, N_SAMPLES).transpose(0, 2, 1)
+
+    with h5py.File(h5_file, "r") as f:
+        igm = f[SITE]["Radiance"]["Metadata"]["Ancillary_Rasters"]["IGM_Data"][:]
+
+    np.testing.assert_allclose(bil, igm, rtol=1e-5)
+
+
+# ---------------------------------------------------------------------------
+# convert_obs
+# ---------------------------------------------------------------------------
+
+
+def test_convert_obs_files_created(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    obs_path, hdr_path = convert_obs(str(h5_file), base)
+    assert Path(obs_path).exists()
+    assert Path(hdr_path).exists()
+
+
+def test_convert_obs_file_size(tmp_path, h5_file):
+    """BIP float32: lines × samples × 10 bands × 4 bytes."""
+    base = str(tmp_path / "out")
+    obs_path, _ = convert_obs(str(h5_file), base)
+    expected = N_LINES * N_SAMPLES * 10 * 4
+    assert Path(obs_path).stat().st_size == expected
+
+
+def test_convert_obs_interleave_is_bip(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_obs(str(h5_file), base)
+    assert "interleave = bip" in Path(hdr_path).read_text()
+
+
+def test_convert_obs_band_names_in_header(tmp_path, h5_file):
+    base = str(tmp_path / "out")
+    _, hdr_path = convert_obs(str(h5_file), base)
+    hdr = Path(hdr_path).read_text()
+    for name in OBS_BAND_NAMES:
+        assert name in hdr
+
+
+def test_convert_obs_time_converted(tmp_path, h5_file):
+    """Band 10 (index 9) should be converted from seconds to decimal hours."""
+    base = str(tmp_path / "out")
+    obs_path, _ = convert_obs(str(h5_file), base)
+
+    # BIP layout: (lines, samples, bands) in C order
+    raw = np.fromfile(obs_path, dtype=np.float32).reshape(N_LINES, N_SAMPLES, 10)
+    time_band = raw[:, :, 9]
+
+    # 58500 seconds → 16.25 hours
+    np.testing.assert_allclose(time_band, 16.25, rtol=1e-3)
+
+
+# ---------------------------------------------------------------------------
+# convert_neon_h5 (end-to-end triplet)
+# ---------------------------------------------------------------------------
+
+
+def test_convert_neon_h5_all_files_created(tmp_path, h5_file):
+    paths = convert_neon_h5(str(h5_file), str(tmp_path))
+    for key in ("rad", "rad_hdr", "loc", "loc_hdr", "obs", "obs_hdr"):
+        assert key in paths
+        assert Path(paths[key]).exists(), f"Missing: {key}"
+
+
+def test_convert_neon_h5_subdir(tmp_path, h5_file):
+    """With subdir=True (default) files should land in <out>/<stem>/."""
+    paths = convert_neon_h5(str(h5_file), str(tmp_path))
+    stem = Path(h5_file).stem
+    assert Path(paths["rad"]).parent.name == stem
+
+
+def test_convert_neon_h5_no_subdir(tmp_path, h5_file):
+    """With subdir=False files should land directly in out_root."""
+    paths = convert_neon_h5(str(h5_file), str(tmp_path), subdir=False)
+    assert Path(paths["rad"]).parent == tmp_path