healpix_analyse.directional_filter#

Geographical directional filtering on NESTED HEALPix fields.

This module provides a deterministic directional filtering operator whose spatial semantics are defined by:

physical distance in metres,
geographical forward azimuth,
a requested geographical direction.

It is intended for Earth-observation processing where an operation has a physical direction such as solar azimuth, illumination direction, wind direction, or another geographical bearing.

This is deliberately different from HealPixConv.

HealPixConv represents a gauge-oriented spherical convolution on a rotated stencil. The operator implemented here instead works directly with the real WGS84 geometry between HEALPix cell centres.

For an output/target cell i and a contributing neighbour/source cell j:

                 geographic North
                       ^
                       |
                       |
                 j  *  |
                   /   |
                  /    |
                 /     |
                *------+
                i

distance_ij
    = WGS84 geodesic distance from i to j

bearing_ij
    = WGS84 forward azimuth from i to j

relative_bearing_ij
    = wrap(bearing_ij - requested_azimuth)

The user-supplied kernel is evaluated as:

weight_ij = kernel(
    distance_ij,
    relative_bearing_ij,
)

Azimuth follows the geographical convention used by the gradient API:

0          = North
pi / 2     = East
pi         = South
3 * pi / 2 = West

Angles increase clockwise from geographic North.

Neighbour selection uses a physical radius in metres. HEALPix ring number, Cartesian pixel windows such as 3x3 or 5x5, and array row/column directions are not part of the public spatial semantics.

The geometry is independent of signal values. Neighbour selection and WGS84 relative geometry are therefore delegated to the shared neighbourhood infrastructure in healpix_analyse._neighbourhood.

Functions#

directional_filter(→ numpy.ndarray | torch.Tensor)

Apply a geographical directional filter to a HEALPix scalar field.

Module Contents#

healpix_analyse.directional_filter.directional_filter(values: numpy.ndarray | torch.Tensor, cell_ids: numpy.ndarray, refinement_level: int, *, max_distance_m: float, azimuth_rad: float, kernel: Kernel, normalize: bool = False, domain: numpy.ndarray | None = None, ellipsoid: str = 'WGS84') → numpy.ndarray | torch.Tensor[source]#

Apply a geographical directional filter to a HEALPix scalar field.

Parameters:

values – NumPy array or PyTorch tensor containing the input signal.

The final dimension corresponds one-to-one with cell_ids. Leading dimensions are preserved, so input may represent one field, several bands, or arbitrary batches.

Examples:
```
(N,)
(bands, N)
(time, bands, N)
```
cell_ids – One-dimensional array of unique NESTED HEALPix cell IDs for which input values are available.
refinement_level – HEALPix refinement level.
max_distance_m – Maximum centre-to-centre geographical distance, in metres, over which cells may contribute.

This is a physical-distance definition. It is intentionally independent of HEALPix ring number, kernel_sz, or Cartesian pixel-window dimensions.
azimuth_rad – Requested geographical azimuth in radians.

The convention is:
```
0          North
pi / 2     East
pi         South
3*pi / 2   West
```
Angles increase clockwise from geographic North.

Values outside [0, 2*pi) are accepted and wrapped.

kernel – Callable defining the directional spatial weight.

It is evaluated as:

kernel(
    distance_m,
    relative_bearing_rad,
)

where relative_bearing_rad is:

wrap(
    bearing_from_target_to_neighbour
    - azimuth_rad
)

in [-pi, +pi).

Therefore:

relative_bearing = 0
    exactly follows the requested direction.

The callable receives one-dimensional NumPy arrays for all valid target-neighbour pairs and must return either one scalar weight or an array broadcastable to the same shape.

Example of a simple forward angular sector:

def forward_sector(distance_m, relative_bearing_rad):
    return (
        np.abs(relative_bearing_rad)
        <= np.deg2rad(20.0)
    ).astype(float)

A distance-dependent directional kernel can combine both arguments:

def directional_gaussian(
    distance_m,
    relative_bearing_rad,
):
    radial = np.exp(
        -0.5 * (distance_m / 200.0) ** 2
    )

    angular = np.exp(
        -0.5
        * (
            relative_bearing_rad
            / np.deg2rad(15.0)
        ) ** 2
    )

    return radial * angular

normalize –

If False (default), return the unnormalised weighted sum:
```
sum(weight * value)
```
If True, divide by the sum of effective weights:
```
sum(weight * value)
-------------------
    sum(weight)
```
Cells outside domain and NaN-valued samples do not contribute to either numerator or denominator.

If the effective weight sum is zero, the normalized result is NaN.
domain – Optional processing and output domain.

domain must be a subset of cell_ids.

If omitted:
```
domain == cell_ids
```
Output ordering follows the exact ordering supplied by domain.

Cells outside domain are absent from the directional operation. They are not interpreted as zero, NaN padding, periodic continuation, or wrapped array indices.
ellipsoid – Geographical ellipsoid used by the shared neighbourhood geometry.

The current relative-geometry implementation supports "WGS84".

Returns:

numpy.ndarray or torch.Tensor – Directionally filtered values.

The output has shape:

values.shape[:-1] + (len(domain),)

NumPy input returns NumPy output.

Torch input returns Torch output on the original device. The weighted signal operation remains differentiable with respect to values.

Notes

The processing pipeline is:

physical neighbourhood
       |
       v
HEALPix neighbour cells
       |
       v
WGS84 relative geometry
distance + forward azimuth
       |
       v
relative bearing
bearing - requested azimuth
       |
       v
user directional kernel
       |
       v
weighted aggregation

The forward bearing is always defined from the output/target cell to the contributing neighbour/source cell:

target i  ----->  source j

Reversing this convention would change an asymmetric directional kernel by approximately 180 degrees and is therefore not equivalent.

The geographical bearing of a cell to itself is undefined. For zero-distance self contributions, relative_bearing_rad is defined as zero by convention.

NaN input samples are treated as unavailable observations rather than propagating through the entire neighbourhood. Their corresponding weights are removed. With normalize=True the remaining valid weights are renormalized.

This operator does not attempt to reproduce Cartesian roll, shifted arrays, image rotation, or a gauge-oriented convolution stencil. Its contract is the underlying geographical operation:

physical distance + geographical direction + HEALPix geometry.

Examples

Apply a forward-looking directional mean over a 500 metre physical neighbourhood:

import numpy as np

from healpix_analyse import directional_filter

def forward_kernel(
    distance_m,
    relative_bearing_rad,
):
    forward = (
        np.abs(relative_bearing_rad)
        <= np.deg2rad(20.0)
    )

    return forward.astype(float)

filtered = directional_filter(
    values,
    cell_ids,
    refinement_level=17,
    max_distance_m=500.0,
    azimuth_rad=np.deg2rad(135.0),
    kernel=forward_kernel,
    normalize=True,
)

In the example above, 135 degrees has its normal geographical meaning: southeast, measured clockwise from North.