Bug in xisoslice

[bjornaa]

Hi,

I am generally very impressed by the xisoslice function. It is a good example of what can be done
within the xarray framework. The new examples show that it is widely useful.

However, if the iso-value is found exactly in the iso-array it returns not-a-number.
I consider this a bug.

My attempt for an easy fix was to change the test for negativity to <= 0.
Then it returns a value, unfortunately in error.

There is also a more arguable bug if the iso-value is an upper or lower boundary of the
iso-array. Also here NaN is returned. This may be the expected behaviour.

Below is the test function a wrote to understand the behaviour of xisoslice.
The bug is caught by the 2nd test, test_no_interpolation.
The questionable boundary handling is exposed by the 3rd test.

------------ 8< ------------------------

import numpy as np
import xarray as xr
import xroms
import pytest

Z = np.linspace(-50, 0, 51)
A = np.sin(Z)
Zdims = ["Z"]
Zcoords = dict(Z=Z)

1D tests

def test_linear_OK():
"""Interpolation is ordinay linear interpolation"""
iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = -22.3
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == projected_array.interp(Z=iso_value)

def test_no_interpolation():
"""Works correctly if the iso_value is already an element in the iso_array"""
iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = Z[11]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == projected_array[11]

def test_boundary():
"""Bounday values are handled correctly"""
iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
# Upper boundary
iso_value = Z[-1]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == A[-1]
# Lower boundary
iso_value = Z[0]
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == A[0]

def test_value_out_of_range():
"""Should return nan if value is out of range"""
iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
iso_value = 3.14 # Too high
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)
iso_value = -70 # Too low
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)

def test_non_monotone():
"""Interpolate where unique, not-a-number if not unique"""

iso_array = xr.DataArray(-np.abs(Z + 10), dims=Zdims, coords=Zcoords)
projected_array = xr.DataArray(np.arange(51), dims="Z", coords=Zcoords)
# Unique iso-"surface", index = 40 + iso_value
# other "solution" 40 - iso_value is out of range
iso_value = -10.5
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert proj_val == 40 + iso_value
# Non-unique iso-"surface", index = 40+iso_value and 40-iso_value
# Refuse to guess
iso_value = -5.5
proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
assert np.isnan(proj_val)

[bjornaa]

The code in test_xisoslice.py came out badly.
Here is a new try,

Bjørn

import numpy as np
import xarray as xr
import xroms
import pytest


Z = np.linspace(-50, 0, 51)
A = np.sin(Z)
Zdims = ["Z"]
Zcoords = dict(Z=Z)

# 1D tests

def test_linear_OK():
    """Interpolation is ordinay linear interpolation"""
    iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
    projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
    iso_value = -22.3
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert proj_val == projected_array.interp(Z=iso_value)


def test_no_interpolation():
    """Works correctly if the iso_value is already an element in the iso_array"""
    iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
    projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
    iso_value = Z[11]
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert proj_val == projected_array[11]


def test_boundary():
    """Bounday values are handled correctly"""
    iso_array = xr.DataArray(Z, dims=Zdims, coords=Zcoords)
    projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
    # Upper boundary
    iso_value = Z[-1]
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert proj_val == A[-1]
    # Lower boundary
    iso_value = Z[0]
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert proj_val == A[0]


def test_value_out_of_range():
    """Should return nan if value is out of range"""
    iso_array = xr.DataArray(Z, dims="Z", coords=Zcoords)
    projected_array = xr.DataArray(A, dims="Z", coords=Zcoords)
    iso_value = 3.14  # Too high
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert np.isnan(proj_val)
    iso_value = -70  # Too low
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert np.isnan(proj_val)


def test_non_monotone():
    """Interpolate where unique, not-a-number if not unique"""

    iso_array = xr.DataArray(-np.abs(Z + 10), dims=Zdims, coords=Zcoords)
    projected_array = xr.DataArray(np.arange(51), dims="Z", coords=Zcoords)
    # Unique iso-"surface", index = 40 + iso_value
    # other "solution" 40 - iso_value is out of range
    iso_value = -10.5
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert proj_val == 40 + iso_value
    # Non-unique iso-"surface", index = 40+iso_value and 40-iso_value
    # Refuse to guess
    iso_value = -5.5
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")
    assert np.isnan(proj_val)


# 2D


def test_2D():
    """Works correctly in two dimensions"""

    iso_array = xr.DataArray(
        np.array([Z, 2 * Z, 3.14 * Z]).transpose(),
        dims=["Z", "X"],
        coords=dict(X=[1, 2, 4], Z=Z),
    )
    projected_array = np.sin(iso_array)
    iso_value = -60.5
    proj_val = xroms.xisoslice(iso_array, iso_value, projected_array, "Z")

    # Result is an xarray of correct shape
    assert type(proj_val) == xr.DataArray
    assert proj_val.dims == ("X",)

    # Out of range for first column
    assert np.isnan(proj_val[0])
    # Linear interpolation for the others
    assert np.isclose(
        proj_val.isel(X=1),
        np.interp(iso_value, iso_array.isel(X=1), projected_array.isel(X=1)),
    )
    assert np.isclose(
        proj_val.isel(X=2),
        np.interp(iso_value, iso_array.isel(X=2), projected_array.isel(X=2)),
    )

[kthyng]

Hi @bjornaa! I am sorry that I missed this issue, and the other couple too. I didn't have notifications set up and I was unaware that they occurred. The past 1.5 years or so I have been at a new job and it's been hard to make time for xroms, but I'm starting to carve out some time! So, please know that I intend to start chipping away at issues here and improve the package. Thanks for your patience, interest, and details issue reports.

[kthyng]

Hi again! I've recently put out a new version to PyPI (v0.2.4, still coming through on conda-forge). Please check to see if this issue is still present in the new version so that over time we can work to address these. Thanks!