Griddata Module¶

The georeader.griddata module provides functions for interpolating scattered geographic data (irregularly-sampled points with per-pixel coordinates) onto regular grids. This is essential for orthorectifying swath-based satellite data like hyperspectral sensors.

Overview¶

Many satellite sensors, particularly pushbroom and whiskbroom scanners, produce data where each pixel has its own geographic coordinates rather than following a regular grid. To analyze this data in GIS software or combine with other datasets, you need to resample it to a regular grid.

┌─────────────────────────────────────────────────────────────────────────────┐
│              IRREGULAR vs REGULAR GRID REPRESENTATION                        │
│                                                                              │
│   Sensor Swath (Irregular)              Orthorectified (Regular)            │
│   ────────────────────────              ────────────────────────            │
│                                                                              │
│       ●  ●   ●  ●                        ┌──┬──┬──┬──┬──┐                   │
│     ●    ●  ●    ●                       ├──┼──┼──┼──┼──┤                   │
│      ●   ● ●  ●                          ├──┼──┼──┼──┼──┤                   │
│    ●   ●    ●   ●                        ├──┼──┼──┼──┼──┤                   │
│                                          └──┴──┴──┴──┴──┘                   │
│                                                                              │
│   Each pixel has (lon, lat)              Fixed affine transform             │
│   from attitude/ephemeris data           pixel (i,j) → (x,y) = T × (i,j)   │
│                                                                              │
│   Causes:                                Benefits:                           │
│   - Sensor scan geometry                 - GIS compatible                   │
│   - Platform motion                      - Easy reprojection                │
│   - Terrain relief                       - Stack multiple images            │
│   - Earth curvature                      - Standard analysis tools          │
│                                                                              │
└─────────────────────────────────────────────────────────────────────────────┘

When to Use This Module¶

Scenario	Use griddata module?	Alternative
Hyperspectral with per-pixel coords	✅ Yes	-
Swath data with lat/lon arrays	✅ Yes	-
Point observations (weather stations)	✅ Yes	-
Regular grid → different CRS	❌ No	`georeader.read.read_reproject`
EMIT with GLT file	⚠️ Use GLT	`georreference()` is faster

Interpolation Methods¶

The module uses scipy.interpolate.griddata internally, which supports three interpolation methods:

┌─────────────────────────────────────────────────────────────────────────────┐
│                        INTERPOLATION METHODS                                 │
│                                                                              │
│   "nearest"                  "linear"                   "cubic" (default)   │
│   ──────────                 ────────                   ─────────────────   │
│                                                                              │
│   ■ ■ ■ ■ ■ ■               ╱╲  ╱╲  ╱╲                   ∿∿∿∿∿∿∿∿∿∿       │
│   ■ ■ ■ ■ ■ ■              ╱  ╲╱  ╲╱  ╲                                     │
│   ■ ■ ■ ■ ■ ■                                                              │
│                                                                              │
│   Voronoi cells            Barycentric on              Clough-Tocher       │
│   Nearest neighbor         Delaunay triangles          C² smooth surface   │
│                                                                              │
│   Continuity: C⁰           Continuity: C⁰              Continuity: C²      │
│   Speed: Fast              Speed: Medium               Speed: Slow         │
│                                                                              │
│   Use for:                 Use for:                    Use for:             │
│   - Classification maps    - Quick previews            - Radiance/Refl     │
│   - Masks                  - When speed matters        - Smooth data       │
│   - Categorical data       - Large datasets            - Final products    │
│                                                                              │
└─────────────────────────────────────────────────────────────────────────────┘

Key Functions¶

`read_to_crs` - Simple Orthorectification¶

The easiest way to orthorectify data when you just want a specific resolution:

from georeader.griddata import read_to_crs
import numpy as np

# Hyperspectral radiance with per-pixel coordinates
radiance = np.random.rand(1000, 1000, 285)  # (H, W, bands)
lons = np.load("pixel_longitudes.npy")       # (1000, 1000)
lats = np.load("pixel_latitudes.npy")        # (1000, 1000)

# Orthorectify to 30m UTM grid (auto-detects UTM zone)
ortho = read_to_crs(
    radiance, lons, lats,
    resolution_dst=30.0,      # 30 meters
    method="cubic"            # Smooth interpolation
)

print(f"Input shape: {radiance.shape}")   # (1000, 1000, 285) - HWC
print(f"Output shape: {ortho.shape}")     # (285, H_out, W_out) - CHW
print(f"Output CRS: {ortho.crs}")         # e.g., EPSG:32610 (UTM Zone 10N)

`read_reproject_like` - Match Existing Grid¶

Orthorectify to match an existing dataset's grid exactly:

from georeader.griddata import read_reproject_like

# Load reference dataset (e.g., Sentinel-2)
reference = GeoTensor.load_file("sentinel2_tile.tif")

# Orthorectify hyperspectral to match Sentinel-2 grid
ortho = read_reproject_like(
    radiance, lons, lats,
    data_like=reference,     # Match this grid
    method="cubic"
)

# ortho now has same CRS, resolution, and extent as reference

`reproject` - Full Control¶

When you need complete control over output parameters:

import rasterio
from georeader.griddata import reproject

# Define exact output grid
transform = rasterio.transform.from_origin(
    west=550000,    # UTM easting
    north=4200000,  # UTM northing  
    xsize=30,       # 30m pixel width
    ysize=30        # 30m pixel height
)

ortho = reproject(
    radiance, lons, lats,
    width=1000,
    height=1000,
    transform=transform,
    dst_crs="EPSG:32610",
    method="cubic",
    fill_value_default=-9999
)

`georreference` - GLT-Based (Fast)¶

For sensors that provide Geolocation Lookup Tables (like EMIT), use this for exact pixel mapping without interpolation:

from georeader.griddata import georreference

# GLT maps output pixel → sensor pixel
# glt[0, i, j] = source column
# glt[1, i, j] = source row  
glt = GeoTensor(glt_array, transform=output_transform, crs="EPSG:32610")

# Fast exact orthorectification (no interpolation)
ortho = georreference(glt, radiance)

GLT vs Interpolation¶

┌─────────────────────────────────────────────────────────────────────────────┐
│              GLT ORTHORECTIFICATION vs INTERPOLATION                         │
│                                                                              │
│   Geolocation Lookup Table (GLT)         Interpolation (griddata)           │
│   ──────────────────────────────         ────────────────────────           │
│                                                                              │
│   ┌───────────────┐                      ┌───────────────┐                  │
│   │ Sensor Array  │                      │ Sensor Array  │                  │
│   │  ┌─┬─┬─┐      │                      │  ●  ●  ●      │                  │
│   │  │A│B│C│      │                      │   ●  ●  ●     │                  │
│   │  └─┴─┴─┘      │                      │  ●  ●  ●      │                  │
│   └───────┬───────┘                      └───────┬───────┘                  │
│           │ GLT lookup                           │ Interpolate              │
│           ▼                                      ▼                          │
│   ┌───────────────┐                      ┌───────────────┐                  │
│   │ Output Grid   │                      │ Output Grid   │                  │
│   │  ┌─┬─┬─┐      │                      │  ┌─┬─┬─┐      │                  │
│   │  │ │A│B│      │                      │  │≈│≈│≈│      │                  │
│   │  └─┴─┴─┘      │                      │  └─┴─┴─┘      │                  │
│   └───────────────┘                      └───────────────┘                  │
│                                                                              │
│   Pros:                                  Pros:                              │
│   ✓ Exact pixel values preserved         ✓ Works without GLT               │
│   ✓ Very fast (array indexing)           ✓ Smooth output                   │
│   ✓ No resampling artifacts              ✓ Any output resolution           │
│                                                                              │
│   Cons:                                  Cons:                              │
│   ✗ Requires GLT from data provider      ✗ Changes pixel values            │
│   ✗ Fixed output resolution              ✗ Slower (O(n log n) Delaunay)    │
│   ✗ May have gaps                        ✗ Edge artifacts possible         │
│                                                                              │
│   Use when:                              Use when:                          │
│   - GLT available (EMIT, etc.)           - No GLT available                │
│   - Exact values needed                  - Custom output resolution        │
│   - Processing derived products          - Point data sources              │
│                                                                              │
└─────────────────────────────────────────────────────────────────────────────┘

Complete Example: EMIT-style Workflow¶

import numpy as np
from georeader.griddata import read_to_crs, georreference, get_shape_transform_crs
from georeader.geotensor import GeoTensor

def process_hyperspectral_swath(radiance, lons, lats, glt=None, 
                                 resolution=30.0, method="cubic"):
    """
    Orthorectify hyperspectral swath data.

    Args:
        radiance: Raw radiance (H, W, C) or (H, W)
        lons: Per-pixel longitudes (H, W)
        lats: Per-pixel latitudes (H, W)
        glt: Optional GLT array (2, H_out, W_out) 
        resolution: Output resolution in meters
        method: Interpolation method if no GLT

    Returns:
        Orthorectified GeoTensor
    """
    if glt is not None:
        # Fast path: use GLT
        print("Using GLT-based orthorectification (exact)")

        # Transpose to (C, H, W) if needed
        if len(radiance.shape) == 3 and radiance.shape[-1] < radiance.shape[0]:
            radiance = np.transpose(radiance, (2, 0, 1))

        return georreference(glt, radiance)

    else:
        # Slow path: interpolation
        print(f"Using {method} interpolation")
        return read_to_crs(
            radiance, lons, lats,
            resolution_dst=resolution,
            method=method,
            fill_value_default=np.nan
        )


# Example usage
radiance = np.random.rand(1000, 1000, 100)
lons = np.linspace(-122.5, -122.0, 1000)[None, :].repeat(1000, axis=0)
lats = np.linspace(37.5, 38.0, 1000)[:, None].repeat(1000, axis=1)

# Add some irregularity (simulating real sensor geometry)
lons += np.random.normal(0, 0.001, lons.shape)
lats += np.random.normal(0, 0.001, lats.shape)

ortho = process_hyperspectral_swath(radiance, lons, lats, resolution=30.0)
print(f"Output: {ortho.shape}, CRS: {ortho.crs}")

Performance Tips¶

Use cubic interpolation sparingly: It's O(n log n) for Delaunay + O(n) per query. For large arrays (>10M points), consider downsampling first.
GLT is always faster: If available, georreference() is a simple array lookup.
Match grids efficiently: Use read_reproject_like instead of computing output parameters manually.
Handle fill values: Areas outside the convex hull of input points will be filled with fill_value_default. Consider using NaN for easy masking.

API Reference¶

`georeader.griddata` ¶

Irregular Grid Interpolation and Georeferencing Module.

This module provides functions for interpolating scattered (non-gridded) geographic data onto regular grids, and for applying geolocation lookup tables (GLT).

Coordinate Systems & Grid Types¶

::

┌─────────────────────────────────────────────────────────────────────────┐
│   IRREGULAR vs REGULAR GRIDS                                            │
│                                                                         │
│   Irregular (Swath/Sensor)           Regular (Orthorectified)           │
│   ─────────────────────────           ──────────────────────            │
│                                                                         │
│       ●  ●   ●  ●                     ┌──┬──┬──┬──┐                     │
│     ●    ●  ●    ●                    ├──┼──┼──┼──┤                     │
│      ●   ● ●  ●                       ├──┼──┼──┼──┤                     │
│    ●   ●    ●   ●                     ├──┼──┼──┼──┤                     │
│                                       └──┴──┴──┴──┘                     │
│                                                                         │
│   Each pixel has unique (lon, lat)    Fixed transform: pixel → geo      │
│   Spacing varies with scan angle      Uniform spacing, axis-aligned     │
│   Common in: pushbroom sensors        Required for: GIS, web maps       │
│                                                                         │
└─────────────────────────────────────────────────────────────────────────┘

Interpolation Methods¶

The module uses :func:scipy.interpolate.griddata for interpolation:

::

┌────────────────────────────────────────────────────────────────────────┐
│  INTERPOLATION METHOD COMPARISON                                       │
│                                                                        │
│  Method      │ Continuity │ Speed  │ Best For                          │
│  ────────────┼────────────┼────────┼─────────────────────────────────  │
│  "nearest"   │ C⁰         │ Fast   │ Categorical data, masks           │
│  "linear"    │ C⁰         │ Medium │ Simple surfaces, quick preview    │
│  "cubic"     │ C²         │ Slow   │ Smooth continuous data (default)  │
│                                                                        │
│  C⁰ = continuous but not differentiable (may have sharp edges)         │
│  C² = smooth, twice differentiable (recommended for radiance/refl)     │
│                                                                        │
└────────────────────────────────────────────────────────────────────────┘

Geolocation Lookup Tables (GLT)¶

Some sensors (like NASA EMIT) provide a GLT array that maps output grid pixels to input sensor pixels. This is faster than interpolation for orthorectification:

::

┌────────────────────────────────────────────────────────────────────────┐
│  GLT-BASED ORTHORECTIFICATION                                          │
│                                                                        │
│  Sensor Array (irregular)              Output Grid (regular)           │
│  ┌───────────────────────┐             ┌──┬──┬──┬──┬──┐                │
│  │  0   1   2   3   ...  │             │  │  │  │  │  │                │
│  │                       │             ├──┼──┼──┼──┼──┤                │
│  │ [r,c] = radiance      │     GLT     │  │██│██│  │  │                │
│  │                       │  ────────►  ├──┼──┼──┼──┼──┤                │
│  │                       │             │  │██│██│██│  │                │
│  └───────────────────────┘             └──┴──┴──┴──┴──┘                │
│                                                                        │
│  GLT[0, i, j] = column in sensor array                                 │
│  GLT[1, i, j] = row in sensor array                                    │
│  output[i, j] = sensor[GLT[1,i,j], GLT[0,i,j]]                         │
│                                                                        │
└────────────────────────────────────────────────────────────────────────┘

Module Functions¶

Grid Interpolation

:func:reproject: Core interpolation from lon/lat arrays to regular grid
:func:read_reproject_like: Match grid to existing GeoData
:func:read_to_crs: Auto-compute grid for given resolution

Grid Utilities

:func:meshgrid: Generate coordinate arrays from transform
:func:get_shape_transform_crs: Compute output grid parameters
:func:footprint: Bounding polygon from lon/lat arrays

GLT Operations

:func:georreference: Apply GLT for fast orthorectification

Example Workflow¶

Orthorectify PRISMA-style data with per-pixel coordinates::

import numpy as np
from georeader.griddata import read_to_crs

# Hyperspectral radiance with irregular coordinates
radiance = np.random.rand(1000, 1000, 285)  # (H, W, bands)
lons = np.random.uniform(-122.5, -122.3, (1000, 1000))  # irregular
lats = np.random.uniform(37.7, 37.9, (1000, 1000))      # irregular

# Interpolate to regular 30m UTM grid
ortho = read_to_crs(
    radiance, lons, lats,
    resolution_dst=30.0,  # 30 meters
    method="cubic"        # smooth interpolation
)
# ortho.shape: (285, H_out, W_out) - regular grid
# ortho.crs: auto-detected UTM zone
# ortho.transform: proper affine transform

References¶

SciPy griddata: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
NASA EMIT L2A Products: https://lpdaac.usgs.gov/products/emitl2arflv001/

`reproject(data, lons, lats, width, height, transform, dst_crs, crs='EPSG:4326', fill_value_default=-1, method=METHOD_DEFAULT)` ¶

Interpolate scattered data to a regular georeferenced grid.

This is the core function for converting irregularly-sampled geographic data (e.g., from pushbroom sensors or point observations) to a regular grid suitable for analysis and visualization.

Algorithm Overview¶

::

┌────────────────────────────────────────────────────────────────────┐
│  INTERPOLATION WORKFLOW                                            │
│                                                                    │
│  1. Flatten inputs                                                 │
│     data: (H, W, C) → (H×W, C)  [or (H, W) → (H×W,)]               │
│     lons/lats: (H, W) → (H×W,)                                     │
│                                                                    │
│  2. Generate output coordinate grid                                │
│     meshgrid(transform, width, height) → (xs, ys)                  │
│     Transform xs, ys from dst_crs to input crs if different        │
│                                                                    │
│  3. Call scipy.interpolate.griddata                                │
│     points = (lons_flat, lats_flat)                                │
│     values = data_flat                                             │
│     xi = (xs_grid, ys_grid)                                        │
│     result = griddata(points, values, xi, method=method)           │
│                                                                    │
│  4. Reshape and handle nodata                                      │
│     Fill NaN regions with fill_value_default                       │
│     Transpose to (C, H, W) if multi-band                           │
│                                                                    │
└────────────────────────────────────────────────────────────────────┘

Interpolation Methods¶

"nearest": Voronoi cell assignment. Fast but produces blocky output. Use for categorical data or masks.
"linear": Barycentric interpolation on Delaunay triangulation. Continuous but not smooth (C⁰ continuity).
"cubic" (default): Clough-Tocher scheme on Delaunay triangulation. Smooth and twice-differentiable (C² continuity). Best for continuous data like radiance or reflectance. Slower than linear.

Parameters:

Name	Type	Description	Default
`data`	`NDArray`	Input array, either: - 2D: (H, W) single-band image - 3D: (H, W, C) multi-band image with C channels Note: This is height × width × channels order, not CHW!	required
`lons`	`NDArray`	Longitude coordinates for each pixel, shape (H, W). Must be same spatial shape as data.	required
`lats`	`NDArray`	Latitude coordinates for each pixel, shape (H, W). Must be same spatial shape as data.	required
`width`	`int`	Output grid width in pixels.	required
`height`	`int`	Output grid height in pixels.	required
`transform`	`Affine`	Output affine transform mapping pixel to CRS coordinates.	required
`dst_crs`	`Any`	Output coordinate reference system (e.g., "EPSG:32610").	required
`crs`	`Optional[Any]`	CRS of input lon/lat arrays. Default "EPSG:4326" (WGS84).	`'EPSG:4326'`
`fill_value_default`	`float`	Value for pixels outside convex hull of input points. Default -1.	`-1`
`method`	`str`	Interpolation method: "nearest", "linear", or "cubic". Default "cubic".	`METHOD_DEFAULT`

Returns:

Type	Description
`GeoTensor`	GeoTensor with shape (H, W) or (C, H, W), georeferenced with
`GeoTensor`	`transform` and `dst_crs`.

Raises:

Type	Description
`ValueError`	If data is not 2D or 3D.

Examples¶

Orthorectify single-band thermal data::

>>> import numpy as np
>>> from georeader.griddata import reproject
>>> import rasterio
>>> 
>>> # Simulated thermal data with irregular coords
>>> temperature = np.random.uniform(280, 320, (100, 100))  # Kelvin
>>> lons = np.linspace(-122.5, -122.3, 100)[None, :] + np.random.normal(0, 0.001, (100, 100))
>>> lats = np.linspace(37.7, 37.9, 100)[:, None] + np.random.normal(0, 0.001, (100, 100))
>>> 
>>> # Define output grid (UTM Zone 10N, 100m resolution)
>>> transform = rasterio.transform.from_origin(550000, 4200000, 100, 100)
>>> ortho = reproject(temperature, lons, lats, 
...                   width=200, height=200, 
...                   transform=transform, 
...                   dst_crs="EPSG:32610")
>>> print(ortho.shape, ortho.crs)
(200, 200) EPSG:32610

Multi-band hyperspectral orthorectification::

>>> # Shape: (H, W, bands) - note the axis order!
>>> radiance = np.random.rand(1000, 1000, 285)
>>> lons = np.load("pixel_lons.npy")  # (1000, 1000)
>>> lats = np.load("pixel_lats.npy")  # (1000, 1000)
>>> 
>>> ortho = reproject(radiance, lons, lats,
...                   width=500, height=500,
...                   transform=my_transform,
...                   dst_crs="EPSG:32610",
...                   method="cubic")
>>> # Output shape: (285, 500, 500) - transposed to (C, H, W)

Warning¶

Input data must be (H, W) or (H, W, C), NOT (C, H, W)
Cubic interpolation can be slow for large arrays (>1M points)
Output pixels outside the convex hull of input points get fill_value

`read_to_crs(data, lons, lats, resolution_dst, dst_crs=None, fill_value_default=-1, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

Reprojects data to the given dst_crs figuring out the transform and shape.

Parameters:

Name	Type	Description	Default
`data`	`Array`	2D or 3D in the form (H, W, bands)	required
`lons`	`Array`	2D array of longitudes (H, W).	required
`lats`	`Array`	2D array of latitudes (H, W).	required
`resolution_dst`	`Union[float, Tuple[float, float]]`	Output resolution	required
`dst_crs`	`Optional[Any]`	Output crs. If None, the dst_crs will be the UTM crs of the center of the data. Defaults to None.	`None`
`fill_value_default`	`float`	fill value. Defaults to -1.	`-1`
`crs`	`_type_`	Input crs. Defaults to "EPSG:4326".	`'EPSG:4326'`
`method`	`str`	Interpolation method. Defaults to "cubic". One of "nearest", "linear", "cubic". (See https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata)	`METHOD_DEFAULT`

Returns:

Name	Type	Description
`GeoTensor`	`GeoTensor`	with reprojected data

Source code in georeader/griddata.py

def read_to_crs(data:NDArray, lons: NDArray, lats:NDArray, 
                resolution_dst:Union[float, Tuple[float,float]], 
                dst_crs:Optional[Any]=None,fill_value_default:float=-1,
                crs:Optional[Any]="EPSG:4326",
                method:str=METHOD_DEFAULT) -> GeoTensor:
    """
    Reprojects data to the given dst_crs figuring out the transform and shape.

    Args:
        data (Array): 2D or 3D in the form (H, W, bands)
        lons (Array): 2D array of longitudes (H, W).
        lats (Array): 2D array of latitudes (H, W).
        resolution_dst (Union[float, Tuple[float,float]]): Output resolution
        dst_crs (Optional[Any], optional): Output crs. If None, 
            the dst_crs will be the UTM crs of the center of the data. Defaults to None.
        fill_value_default (float, optional): fill value. Defaults to -1.
        crs (_type_, optional): Input crs. Defaults to "EPSG:4326".
        method (str, optional): Interpolation method. Defaults to "cubic". One of
            "nearest", "linear", "cubic". (See https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata)

    Returns:
        GeoTensor: with reprojected data
    """
    width, height, transform, dst_crs = get_shape_transform_crs(lons, lats, 
                                                                resolution_dst=resolution_dst, 
                                                                dst_crs=dst_crs, crs=crs)

    return reproject(data, lons=lons, lats=lats, width=width, 
                     height=height, transform=transform, dst_crs=dst_crs,
                     fill_value_default=fill_value_default,
                     crs=crs, method=method)

`read_reproject_like(data, lons, lats, data_like, resolution_dst=None, fill_value_default=None, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

Reprojects data to the same crs, transform and shape as data_like

Parameters:

Name	Type	Description	Default
`data`	`Array`	input data 2D or 3D in the form (height, width, bands)	required
`lons`	`Array`	2D array of longitudes	required
`lats`	`Array`	2D array of latitudes	required
`data_like`	`GeoData`	GeoData to reproject to	required
`resolution_dst`	`Optional[Union[float, Tuple[float, float]]]`	If provided, the output resolution will be set to this value. Otherwise, the output resolution will be the same as data_like. Defaults to None.	`None`
`fill_value_default`	`Optional[float]`	fill value. Defaults to None.	`None`
`crs`	`Optional[Any]`	Input crs. Defaults to "EPSG:4326".	`'EPSG:4326'`
`method`	`str`	Interpolation method. Defaults to "cubic". One of "nearest", "linear", "cubic".	`METHOD_DEFAULT`

Returns:

Name	Type	Description
`GeoTensor`	`GeoTensor`	with reprojected data

Source code in georeader/griddata.py

def read_reproject_like(data:NDArray, lons: NDArray, lats:NDArray, 
                        data_like:GeoData, resolution_dst:Optional[Union[float, Tuple[float,float]]]=None,
                        fill_value_default:Optional[float]=None,
                        crs:Optional[Any]="EPSG:4326",
                        method:str=METHOD_DEFAULT) -> GeoTensor:
    """
    Reprojects data to the same crs, transform and shape as data_like

    Args:
        data (Array): input data 2D or 3D in the form (height, width, bands)
        lons (Array): 2D array of longitudes
        lats (Array): 2D array of latitudes
        data_like (GeoData): GeoData to reproject to
        resolution_dst (Optional[Union[float, Tuple[float,float]]], optional): If provided, the output
            resolution will be set to this value. Otherwise, the output resolution will be the same
            as data_like. Defaults to None.
        fill_value_default (Optional[float], optional): fill value. Defaults to None.
        crs (Optional[Any], optional): Input crs. Defaults to "EPSG:4326".
        method (str, optional): Interpolation method. Defaults to "cubic". One of
            "nearest", "linear", "cubic".

    Returns:
        GeoTensor: with reprojected data
    """
    width = data_like.shape[-1]
    height = data_like.shape[-2]
    transform = data_like.transform
    dst_crs = data_like.crs
    if resolution_dst is not None:
        transform = transform_to_resolution_dst(transform, resolution_dst)

    fill_value_default = fill_value_default or data_like.fill_value_default
    return reproject(data, lons, lats, width, height, transform, dst_crs, 
                     fill_value_default=fill_value_default, crs=crs,
                     method=method)

`georreference(glt, data, valid_glt=None, fill_value_default=None)` ¶

Apply a Geolocation Lookup Table (GLT) to orthorectify sensor data.

This function performs fast, exact orthorectification by using a pre-computed lookup table that maps output grid pixels to input sensor pixels. Unlike interpolation-based methods, GLT orthorectification preserves original pixel values without resampling artifacts.

GLT Structure¶

::

┌────────────────────────────────────────────────────────────────────┐
│  GLT ARRAY STRUCTURE                                               │
│                                                                    │
│  glt.shape = (2, H_out, W_out)                                     │
│                                                                    │
│  glt[0, i, j] = source column (x-index in sensor array)            │
│  glt[1, i, j] = source row    (y-index in sensor array)            │
│                                                                    │
│  For each output pixel (i, j):                                     │
│    output[..., i, j] = data[..., glt[1,i,j], glt[0,i,j]]           │
│                                                                    │
│  ┌──────────────────────┐      ┌──────────────────────┐            │
│  │    Sensor Array     │       │   Output Grid        │            │
│  │    (raw data)       │       │   (orthorectified)   │            │
│  │   ┌───┬───┬───┐     │       │  ┌──┬──┬──┬──┐       │            │
│  │   │ A │ B │ C │     │       │  │  │ A│ B│  │       │            │
│  │   ├───┼───┼───┤     │ GLT   │  ├──┼──┼──┼──┤       │            │
│  │   │ D │ E │ F │  ───────►   │  │  │ D│ E│ F│       │            │
│  │   ├───┼───┼───┤     │       │  ├──┼──┼──┼──┤       │            │
│  │   │ G │ H │ I │     │       │  │ G│ H│ I│  │       │            │
│  │   └───┴───┴───┘     │       │  └──┴──┴──┴──┘       │            │
│  └──────────────────────┘      └──────────────────────┘            │
│                                                                    │
│  GLT handles: terrain distortion, sensor geometry, Earth curvature │
│  Invalid pixels: glt values = fill_value_default (typically -1)    │
│                                                                    │
└────────────────────────────────────────────────────────────────────┘

Common Use Cases¶

Post-processing orthorectified products: If you compute spectral indices or run ML inference on sensor-geometry data, use this to orthorectify the results without re-processing the full cube.
Custom band math: Calculate a derived product from raw bands, then apply GLT for geographic alignment.
Mask application: Create masks in sensor space, then georeference to match the ortho grid for overlay.

Parameters:

Name	Type	Description	Default
`glt`	`GeoTensor`	GLT GeoTensor with shape (2, H_out, W_out). Contains the (column, row) indices mapping each output pixel to the sensor array. Must have valid transform and CRS for the output grid.	required
`data`	`NDArray`	Sensor-space data array with shape (H_sensor, W_sensor) or (C, H_sensor, W_sensor). Will be indexed by GLT values.	required
`valid_glt`	`Optional[NDArray]`	Optional boolean mask of shape (H_out, W_out) indicating valid GLT pixels. If None, auto-computed as pixels where both GLT channels differ from fill_value_default.	`None`
`fill_value_default`	`Optional[Union[int, float]]`	Fill value for output pixels with invalid GLT. If None, defaults to 0.	`None`

Returns:

Type	Description
`GeoTensor`	GeoTensor with shape (H_out, W_out) or (C, H_out, W_out) matching
`GeoTensor`	the GLT's spatial dimensions, with transform and CRS from the GLT.

Raises:

Type	Description
`ValueError`	If data shape is not 2D or 3D.

Examples¶

Orthorectify a spectral index computed in sensor space::

>>> import numpy as np
>>> from georeader.griddata import georreference
>>> 
>>> # Load EMIT data (sensor geometry)
>>> emit = EMITImage("EMIT_L2A_file.nc")
>>> radiance = emit.load_raw()  # (285, 1242, 1280) - sensor space
>>> 
>>> # Compute NDVI in sensor space (faster than on ortho grid)
>>> nir = radiance[120]  # ~850nm
>>> red = radiance[60]   # ~665nm  
>>> ndvi = (nir - red) / (nir + red + 1e-6)  # (1242, 1280)
>>> 
>>> # Get GLT and orthorectify
>>> glt = emit.load_glt()  # (2, H_ortho, W_ortho)
>>> ndvi_ortho = georreference(glt, ndvi)
>>> 
>>> print(f"Sensor: {ndvi.shape}, Ortho: {ndvi_ortho.shape}")
Sensor: (1242, 1280), Ortho: (1500, 1600)

Orthorectify multi-band processed data::

>>> # ML model output in sensor space
>>> class_probs = model.predict(radiance)  # (10, 1242, 1280)
>>> 
>>> # Orthorectify all probability bands at once
>>> probs_ortho = georreference(glt, class_probs, fill_value_default=0)
>>> print(probs_ortho.shape)  # (10, H_ortho, W_ortho)

Notes¶

This function does NOT interpolate - it performs exact pixel lookup
Much faster than :func:reproject for sensor→ortho conversion
The GLT must be pre-computed (provided by data producer or computed from RPC/sensor model)
Invalid GLT values (outside valid sensor range) must be marked with fill_value_default in the GLT

`meshgrid(transform, width, height, source_crs=None, dst_crs=None)` ¶

Generate the meshgrid of geographic coordinates from the transform. If source_crs and dst_crs are provided, the meshgrid will be transformed to the dst_crs.

Parameters:

Name	Type	Description	Default
`transform`	`Affine`	transform	required
`width`	`int`	width	required
`height`	`int`	height	required
`source_crs`	`Optional[Any]`	source crs. Defaults to None.	`None`
`dst_crs`	`Optional[Any]`	destination crs. Defaults to None.	`None`

Returns:

Type	Description
`Tuple[NDArray, NDArray]`	Tuple[NDArray, NDArray]: 2D arrays of xs and ys coordinates

Source code in georeader/griddata.py

def meshgrid(transform:rasterio.transform.Affine, width:int, height:int, 
             source_crs:Optional[Any]=None, dst_crs:Optional[Any]=None) -> Tuple[NDArray, NDArray]:
    """
    Generate the meshgrid of geographic coordinates from the transform.
    If source_crs and dst_crs are provided, the meshgrid will be transformed to the dst_crs.

    Args:
        transform (rasterio.transform.Affine): transform
        width (int): width
        height (int): height
        source_crs (Optional[Any], optional): source crs. Defaults to None.
        dst_crs (Optional[Any], optional): destination crs. Defaults to None.

    Returns:
        Tuple[NDArray, NDArray]: 2D arrays of xs and ys coordinates
    """
    cols, rows = np.meshgrid(np.arange(width), np.arange(height))
    xs, ys = rasterio.transform.xy(transform, rows, cols)
    xs= np.array(xs)
    ys = np.array(ys)

    if dst_crs is not None:
        assert source_crs is not None, "source_crs must be provided if dst_crs is provided"
        xs, ys = rasterio.warp.transform(source_crs, dst_crs, xs.ravel(),ys.ravel())
        xs = np.array(xs).reshape(height, width)
        ys = np.array(ys).reshape(height, width)

    return xs, ys

`footprint(lons, lats)` ¶

Returns the Polygon surrounding the given longitudes and latitudes

Parameters:

Name	Type	Description	Default
`lons`	`array`	2D array of longitudes	required
`lats`	`array`	2D array of latitudes	required

Returns:

Name	Type	Description
`Polygon`	`Polygon`	Polygon surrounding the given longitudes and latitudes

Source code in georeader/griddata.py

def footprint(lons:NDArray, lats:NDArray) -> Polygon:
    """
    Returns the Polygon surrounding the given longitudes and latitudes

    Args:
        lons (np.array): 2D array of longitudes
        lats (np.array): 2D array of latitudes

    Returns:
        Polygon: Polygon surrounding the given longitudes and latitudes
    """
    lonsrav = lons.ravel()
    latsrav = lats.ravel()
    idxminlon = np.argmin(lonsrav)
    idxminlat = np.argmin(latsrav)
    idxmaxlon = np.argmax(lonsrav)
    idxmaxlat = np.argmax(latsrav)

    return Polygon([(lonsrav[idx],latsrav[idx]) for idx in [idxminlon, idxminlat, idxmaxlon, idxmaxlat]])

`get_shape_transform_crs(lons, lats, resolution_dst, dst_crs=None, crs='EPSG:4326')` ¶

Get the shape, transform and crs for the given lons and lats and resolution_dst.

Parameters:

Name	Type	Description	Default
`lons`	`NDArray`	2D array of longitudes (H, W).	required
`lats`	`NDArray`	2D array of latitudes (H, W).	required
`resolution_dst`	`Union[float, Tuple[float, float]]`	Output resolution. If a single float is provided, the resolution will be (resolution_dst, resolution_dst).	required
`dst_crs`	`Optional[Any]`	Output crs. If None, the dst_crs will be the UTM crs of the center of the data. Defaults to None.	`None`
`crs`	`Any`	Input crs. Defaults to "EPSG:4326".	`'EPSG:4326'`

Returns:

Type	Description
`Tuple[int, int, Affine, Any]`	Tuple[int, int, rasterio.transform.Affine, Any]: width, height, transform and dst_crs.

Source code in georeader/griddata.py

def get_shape_transform_crs(lons: NDArray, lats:NDArray, 
                            resolution_dst:Union[float, Tuple[float,float]], 
                            dst_crs:Optional[Any]=None,
                            crs:Optional[Any]="EPSG:4326") -> Tuple[int, int, rasterio.transform.Affine, Any]:
    """
    Get the shape, transform and crs for the given lons and lats and resolution_dst.    

    Args:
        lons (NDArray): 2D array of longitudes (H, W).
        lats (NDArray): 2D array of latitudes (H, W).
        resolution_dst (Union[float, Tuple[float,float]]): Output resolution.
            If a single float is provided, the resolution will be (resolution_dst, resolution_dst).
        dst_crs (Optional[Any], optional): Output crs. If None,
            the dst_crs will be the UTM crs of the center of the data. Defaults to None.
        crs (Any, optional): Input crs. Defaults to "EPSG:4326".

    Returns:
        Tuple[int, int, rasterio.transform.Affine, Any]: width, height, transform and dst_crs.
    """
    if isinstance(resolution_dst, numbers.Number):
        resolution_dst = (abs(resolution_dst), abs(resolution_dst))

    # Figure out UTM crs
    if dst_crs is None:
        mean_lat = np.nanmean(lats)
        mean_lon = np.nanmean(lons)
        dst_crs = georeader.get_utm_epsg((mean_lon, mean_lat), 
                                         crs_point_or_geom=crs)

    # Figure out transform
    pol = footprint(lons, lats)
    pol_dst_crs = polygon_to_crs(pol, crs_polygon=crs, dst_crs=dst_crs)
    minx, miny, maxx, maxy = pol_dst_crs.bounds

    # Add the resolution to the max values to get the correct shape.
    maxx = maxx + resolution_dst[0]
    miny = miny - resolution_dst[1]
    transform = rasterio.transform.from_origin(minx, maxy, resolution_dst[0], resolution_dst[1])

    # resolution_dst= res(transform)
    width = math.ceil(abs((maxx -minx) / resolution_dst[0]))
    height = math.ceil(abs((maxy - miny) / resolution_dst[1]))
    return width, height, transform, dst_crs

Griddata Module¶

Overview¶

When to Use This Module¶

Interpolation Methods¶

Key Functions¶

`read_to_crs` - Simple Orthorectification¶

`read_reproject_like` - Match Existing Grid¶

`reproject` - Full Control¶

`georreference` - GLT-Based (Fast)¶

GLT vs Interpolation¶

Complete Example: EMIT-style Workflow¶

Performance Tips¶

API Reference¶

`georeader.griddata` ¶

Coordinate Systems & Grid Types¶

Interpolation Methods¶

Geolocation Lookup Tables (GLT)¶

Module Functions¶

Example Workflow¶

See Also¶

References¶

`reproject(data, lons, lats, width, height, transform, dst_crs, crs='EPSG:4326', fill_value_default=-1, method=METHOD_DEFAULT)` ¶

Algorithm Overview¶

Interpolation Methods¶

Examples¶

Warning¶

See Also¶

`read_to_crs(data, lons, lats, resolution_dst, dst_crs=None, fill_value_default=-1, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

`read_reproject_like(data, lons, lats, data_like, resolution_dst=None, fill_value_default=None, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

`georreference(glt, data, valid_glt=None, fill_value_default=None)` ¶

GLT Structure¶

Common Use Cases¶

Examples¶

Notes¶

See Also¶

`meshgrid(transform, width, height, source_crs=None, dst_crs=None)` ¶

`footprint(lons, lats)` ¶

`get_shape_transform_crs(lons, lats, resolution_dst, dst_crs=None, crs='EPSG:4326')` ¶

See Also¶

Griddata Module¶

Overview¶

When to Use This Module¶

Interpolation Methods¶

Key Functions¶

read_to_crs - Simple Orthorectification¶

read_reproject_like - Match Existing Grid¶

reproject - Full Control¶

georreference - GLT-Based (Fast)¶

GLT vs Interpolation¶

Complete Example: EMIT-style Workflow¶

Performance Tips¶

API Reference¶

georeader.griddata ¶

Coordinate Systems & Grid Types¶

Interpolation Methods¶

Geolocation Lookup Tables (GLT)¶

Module Functions¶

Example Workflow¶

See Also¶

References¶

reproject(data, lons, lats, width, height, transform, dst_crs, crs='EPSG:4326', fill_value_default=-1, method=METHOD_DEFAULT) ¶

Algorithm Overview¶

Interpolation Methods¶

Examples¶

Warning¶

See Also¶

read_to_crs(data, lons, lats, resolution_dst, dst_crs=None, fill_value_default=-1, crs='EPSG:4326', method=METHOD_DEFAULT) ¶

read_reproject_like(data, lons, lats, data_like, resolution_dst=None, fill_value_default=None, crs='EPSG:4326', method=METHOD_DEFAULT) ¶

georreference(glt, data, valid_glt=None, fill_value_default=None) ¶

GLT Structure¶

Common Use Cases¶

Examples¶

Notes¶

See Also¶

meshgrid(transform, width, height, source_crs=None, dst_crs=None) ¶

footprint(lons, lats) ¶

get_shape_transform_crs(lons, lats, resolution_dst, dst_crs=None, crs='EPSG:4326') ¶

See Also¶

`read_to_crs` - Simple Orthorectification¶

`read_reproject_like` - Match Existing Grid¶

`reproject` - Full Control¶

`georreference` - GLT-Based (Fast)¶

`georeader.griddata` ¶

`reproject(data, lons, lats, width, height, transform, dst_crs, crs='EPSG:4326', fill_value_default=-1, method=METHOD_DEFAULT)` ¶

`read_to_crs(data, lons, lats, resolution_dst, dst_crs=None, fill_value_default=-1, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

`read_reproject_like(data, lons, lats, data_like, resolution_dst=None, fill_value_default=None, crs='EPSG:4326', method=METHOD_DEFAULT)` ¶

`georreference(glt, data, valid_glt=None, fill_value_default=None)` ¶

`meshgrid(transform, width, height, source_crs=None, dst_crs=None)` ¶

`footprint(lons, lats)` ¶

`get_shape_transform_crs(lons, lats, resolution_dst, dst_crs=None, crs='EPSG:4326')` ¶