Feature/circular difference obsfunction

[ReubenHill]

Description

Add test files for the circular difference obsfunction added in https://github.com/JCSDA-internal/ufo/pull/4058. The file was created with the below python code:

"""
Test file generator for CircularDifference obsfunction.

This script generates NetCDF test files for testing a CircularDifference
obsfunction that calculates elementwise circular differences between two vectors.

Test cases cover wind direction observations (circular period 360 degrees):
- Two vectors: windDirectionA and windDirectionB
- Testing unsigned circular distance: abs((A - B + P/2) % P - P/2)
- Testing signed circular distance: (B - A + P/2) % P - P/2 (from A to B)
- Cases include wrap-around at 0/360 boundary

The circular difference between two values x and y with period P is:
- Unsigned: abs((x - y + P/2) % P - P/2)
- Signed (from x to y): ((y - x + P/2) % P) - P/2
  This gives the shortest signed rotation from x to y

Test data includes various cases:
- Small differences within same quadrant
- Large differences requiring wrap-around
- Boundary cases near 0 and period
- Missing values
"""

import netCDF4 as nc
import numpy as np
from datetime import datetime as dt, timezone as tz

missing = -3.3687952621450501176e38  # jedi's missing value


def circular_difference_unsigned(a, b, period):
    """Compute the minimal circular distance between two values."""
    half = period / 2.0
    return abs((b - a + half) % period - half)


def circular_difference_signed(a, b, period):
    """Calculate signed circular difference from a to b (b - a)."""
    half = period / 2.0
    return ((b - a + half) % period) - half


def write_circular_difference_obs_file(
    name,
    station_id,
    datetime,
    lat,
    lon,
    wind_direction_a,
    wind_direction_b,
    wind_direction_error,
    unsigned_wind_direction_diff_ref,
    signed_wind_direction_diff_ref,
):
    """
    Generate netCDF test file for CircularDifference obsfunction.

    Args:
        name (str): Name of the output netCDF file.
        station_id (np.ndarray): Array of station IDs (strings).
        datetime (np.ndarray): Array of observation datetimes
            (seconds since epoch).
        lat (np.ndarray): Array of latitudes.
        lon (np.ndarray): Array of longitudes.
        wind_direction_a (np.ndarray): Wind direction A values in
            degrees (0-365).
        wind_direction_b (np.ndarray): Wind direction B values in
            degrees (0-365).
        wind_direction_error (np.ndarray): Wind direction error values.
        unsigned_wind_direction_diff_ref (np.ndarray): Reference
            unsigned circular difference for wind direction.
        signed_wind_direction_diff_ref (np.ndarray): Reference signed
            circular difference for wind direction.
    """
    file = nc.Dataset(name, "w", format="NETCDF4")
    nlocs = len(station_id)
    file.createDimension("Location", nlocs)
    file._ioda_layout = "ObsGroup"
    file._ioda_layout_version = 0
    file.date_time = "20200101T0000Z"

    # Metadata variables
    lat_var = file.createVariable("MetaData/latitude", "f4", ("Location"))
    lat_var[:] = lat

    lon_var = file.createVariable("MetaData/longitude", "f4", ("Location"))
    lon_var[:] = lon

    datetime_var = file.createVariable(
        "MetaData/dateTime", "i8", ("Location")
    )
    datetime_var.units = "seconds since 1970-01-01T00:00:00Z"
    datetime_var[:] = datetime

    station_id_var = file.createVariable(
        "MetaData/stationIdentification", str, ("Location")
    )
    station_id_var[:] = station_id

    # ObsValue variables - wind direction A and B
    wind_direction_a_var = file.createVariable(
        "ObsValue/windDirectionA", "f4", ("Location")
    )
    wind_direction_a_var[:] = wind_direction_a

    wind_direction_b_var = file.createVariable(
        "ObsValue/windDirectionB", "f4", ("Location")
    )
    wind_direction_b_var[:] = wind_direction_b

    wind_direction_error_a_var = file.createVariable(
        "ObsError/windDirectionA", "f4", ("Location")
    )
    wind_direction_error_a_var[:] = wind_direction_error

    wind_direction_error_b_var = file.createVariable(
        "ObsError/windDirectionB", "f4", ("Location")
    )
    wind_direction_error_b_var[:] = wind_direction_error

    # TestReference variables - wind direction circular differences
    unsigned_wind_diff_var = file.createVariable(
        "TestReference/unsignedWindDirectionCircularDiff",
        "f4",
        ("Location"),
    )
    unsigned_wind_diff_var[:] = unsigned_wind_direction_diff_ref

    signed_wind_diff_var = file.createVariable(
        "TestReference/signedWindDirectionCircularDiff", "f4", ("Location")
    )
    signed_wind_diff_var[:] = signed_wind_direction_diff_ref

    file.close()


if __name__ == "__main__":
    # Create observations for elementwise circular difference testing
    # Each observation compares windDirectionA[i] with windDirectionB[i]

    nlocs = 15  # 15 test cases covering various scenarios

    # Station identifiers (all same station for simplicity)
    station_id = np.array(["TestStation"] * nlocs)

    # Create observation datetime (same for all observations)
    start_time = int(dt(2020, 1, 1, 12, 0, 0, tzinfo=tz.utc).timestamp())
    datetime = np.full(nlocs, start_time, dtype=np.int64)

    # Create lat/lon arrays (arbitrary, same for all)
    lat = np.zeros(nlocs)
    lon = np.zeros(nlocs)

    # Wind direction test cases (period = 360 degrees)
    # Test various scenarios for circular difference
    wind_direction_a = np.array([
        10.0,      # Case 1: Small difference, no wrap
        350.0,     # Case 2: Wrap-around near boundary
        0.0,       # Case 3: At lower boundary
        360.0,     # Case 4: At upper boundary (same as 0)
        100.0,     # Case 5: Mid-range
        200.0,     # Case 6: Large forward difference
        50.0,      # Case 7: Large backward difference
        180.0,     # Case 8: Half period (180->0)
        missing,   # Case 9: Missing A value
        150.0,     # Case 10: Missing B value
        270.0,     # Case 11: Quarter period difference
        90.0,      # Case 12: Quarter period difference (opposite)
        1.0,       # Case 13: Near-boundary small difference
        359.0,     # Case 14: Near-boundary small difference
        0.0,       # Case 15: Half period (0->180)
    ])

    wind_direction_b = np.array([
        20.0,      # Case 1: B-A = +10 degrees (unsigned: 10, signed: +10)
        10.0,      # Case 2: B-A = -340 = +20 (unsigned: 20, signed: +20)
        360.0,     # Case 3: B-A = 360 (unsigned: 0, signed: 0)
        5.0,       # Case 4: B-A = -355 = +5 (unsigned: 5, signed: +5)
        250.0,     # Case 5: B-A = +150 (unsigned: 150, signed: +150)
        50.0,      # Case 6: B-A = -150 (unsigned: 150, signed: -150)
        300.0,     # Case 7: B-A = +250 = -110 (unsigned: 110, signed: -110)
        0.0,       # Case 8: B-A = -180 (unsigned: 180, signed: +180 or -180)
        100.0,     # Case 9: Should be missing in result
        missing,   # Case 10: Should be missing in result
        0.0,       # Case 11: B-A = -270 = +90 (unsigned: 90, signed: +90)
        360.0,     # Case 12: B-A = +270 = -90 (unsigned: 90, signed: -90)
        359.0,     # Case 13: B-A = +358 = -2 (unsigned: 2, signed: -2)
        1.0,       # Case 14: B-A = -358 = +2 (unsigned: 2, signed: +2)
        180.0,     # Case 15: B-A = +180 (unsigned: 180, signed: +180 or -180)
    ])

    # Calculate reference circular differences
    wind_period = 360.0

    # Unsigned circular differences
    unsigned_wind_diff = circular_difference_unsigned(
        wind_direction_a, wind_direction_b, wind_period
    )

    # Signed circular differences
    signed_wind_diff = circular_difference_signed(
        wind_direction_a, wind_direction_b, wind_period
    )

    # Handle missing values
    wind_missing_mask = (wind_direction_a == missing) | (
        wind_direction_b == missing
    )
    unsigned_wind_diff[wind_missing_mask] = missing
    signed_wind_diff[wind_missing_mask] = missing

    # Error arrays (arbitrary values)
    wind_direction_error = np.ones(nlocs)

    # Write to netCDF file (no permutation needed for this simple test)
    write_circular_difference_obs_file(
        name="circular_difference_test_obs.nc4",
        station_id=station_id,
        datetime=datetime,
        lat=lat,
        lon=lon,
        wind_direction_a=wind_direction_a,
        wind_direction_b=wind_direction_b,
        wind_direction_error=wind_direction_error,
        unsigned_wind_direction_diff_ref=unsigned_wind_diff,
        signed_wind_direction_diff_ref=signed_wind_diff,
    )

    print("Generated circular_difference_test_obs.nc4")
    print(f"Total observations: {nlocs}")
    print(f"Wind direction period: {wind_period} degrees")
    print("\nTest cases cover:")
    print("- Small and large differences")
    print("- Wrap-around at boundaries")
    print("- Exactly half-period differences")
    print("- Missing value handling")

Issue(s) addressed

N/A

Dependencies

List the other PRs that this PR is dependent on:

• https://github.com/JCSDA-internal/ufo/pull/4058

build-group=https://github.com/JCSDA-internal/ufo/pull/4058

Impact

Expected impact on downstream repositories: None

Checklist

• I have performed a self-review of my own code
• I have made corresponding changes to the documentation
• I have run the unit tests before creating the PR

[ctgh]

Thanks for providing the python code!