makepath / xarray-spatial

Raster-based Spatial Analytics for Python
https://xarray-spatial.readthedocs.io/
MIT License
805 stars 81 forks source link

Overflow error with NDVI calculation #726

Closed jessjaco closed 1 year ago

jessjaco commented 2 years ago

Describe the bug ndvi is giving incorrect results when input arrays are unsigned integer due to overflow in subtraction.

Expected behavior Consider the following

import numpy as np
import xarray as xr
import xrspatial

a = xr.DataArray(np.array([[1,1,1],[1,1,1]], dtype='uint16'))
b = xr.DataArray(np.array([[0,1,2],[0,1,2]], dtype='uint16'))
xrspatial.ndvi(a,b)

Results:

<xarray.DataArray 'ndvi' (dim_0: 2, dim_1: 3)>
array([[1.000000e+00, 0.000000e+00, 6.148915e+18],
       [1.000000e+00, 0.000000e+00, 6.148915e+18]], dtype=float32)
Dimensions without coordinates: dim_0, dim_1

Values in third column are due to overflow of unsigned integer during subtraction, see https://github.com/numpy/numpy/issues/21237. Using the above data, this can be illustrated by

>>> a[0,2] - b[0,2]
<xarray.DataArray ()>
array(65535, dtype=uint16)

A solution could be to call np.subtract directly in e.g. _normalized_ratio_cpu, though that may not be robust to all cases. I think there are some alternatives in the numpy issues linked above

>>> np.subtract(a[0,2], b[0,2], dtype='int32')
<xarray.DataArray ()>
array(-1, dtype=int32)
brendancol commented 1 year ago

@thuydotm can you take a look at this?

brendancol commented 1 year ago

@thuydotm let's put this as high-priority