Select to view content in your preferred language

Really close points

303
0
2 weeks ago
Labels (1)
DanPatterson
MVP Esteemed Contributor
1 0 303

Back to the simple point example of comparing two shapes for common points where there might be floating point issues.  Note that in the b1 array, I introduces a really tiny difference in the second point so that it differs from its previous incarnation.... b

a = np.array([[0., 0.], [1., 4.], [4., 3.], [5., 0.],[0., 0.]])

b = np.array([[1., 0.], [1., 4.], [5., 0.], [1., 0.]])

b1 = np.array([[1., 0.], [1., 4.000001], [5., 0.], [1., 0.]])

So if we compare a to b for equality, then we would have an issue because of the small floating point difference.  To account for this, we need to step up the check to account for this. 

Here is the code.  Read the header.

def a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False):
    r"""Return indices of, or points where, two point arrays are equal.

    Parameters
    ----------
    a, b : 2d arrays
        No error checking so ensure that a and b are 2d, but their shape need
        not be the same.

    Notes
    -----
    Modified from `np.isclose`, stripping out the nan and ma checks.
    Adjust atol and rtol appropriately.

    >>> np.nonzero(a_eq_b(a, b))[0]

    One could use::

    >>> np.equal(A, B).all(1).nonzero()[0]
    >>> np.where((a[:, None] == b).all(-1).any(1))[0]

    If you aren't worried about floating-point issues in equality checks.
    """
    a = np.atleast_2d(a)
    b = np.atleast_2d(b)
    if b.size > 2:
        b = b[:, None]
    w = np.less_equal(np.abs(a - b), atol + rtol * abs(b)).all(-1)
    if w.ndim > 1:
        w = w.any(-1).squeeze()
    if return_pnts:
        return b[w].squeeze()
    return w.squeeze()

The magic is to account for tolerances (covered in the NumPy documentation).  This allows one to identify points that are ridiculously close that they should be considered the same.  If the difference is due to errors, then they can be fixed.

So compare a to b and a to b1 with and without tolerance checks.

# -- a to b comparison, both are correct
a_eq_b(a, b, atol=1.0e-8, rtol=1.0e-5, return_pnts=False)
array([0, 1, 1, 0])

a_eq_b(a, b, atol=0., rtol=0., return_pnts=False)
array([0, 1, 1, 0])

# -- a to b1 comparison, tolerance check is correct
a_eq_b(a, b1, atol=1.0e-8, rtol=1.0e-5, return_pnts=False)
Out[160]: array([0, 1, 1, 0])

# -- exactitude check only identifies one 'equality'
a_eq_b(a, b1, atol=0., rtol=0., return_pnts=False)
Out[161]: array([0, 0, 1, 0])

So when you things are all good, make sure that you account for the small differences.

If coordinates/values are the result of a calculation, it is really important to round/truncate your results to the expected precision that is associated with the rest of your data.

Enough for now.

 

Addendum

see Common points for a fuller discussion

Common points

Tags (3)
Contributors
About the Author
Retired Geomatics Instructor (also DanPatterson_Retired). Currently working on geometry projects (various) as they relate to GIS and spatial analysis. I use NumPy, python and kin and interface with ArcGIS Pro.
Labels