Multispectral tools: convert data to float32 dtype before doing calculations (#755)

* convert to f4 dtype before doing calculations

* add tests

Author: thuydotm

xrspatial/multispectral.py

@@ -146,7 +146,9 @@ def arvi(nir_agg: xr.DataArray,
                                       cupy_func=_arvi_cupy,
                                       dask_cupy_func=_arvi_dask_cupy)
 
-    out = mapper(red_agg)(nir_agg.data, red_agg.data, blue_agg.data)
+    out = mapper(red_agg)(
+        nir_agg.data.astype('f4'), red_agg.data.astype('f4'), blue_agg.data.astype('f4')
+    )
 
     return DataArray(out,
                      name=name,
@@ -311,8 +313,10 @@ def evi(nir_agg: xr.DataArray,
                                       cupy_func=_evi_cupy,
                                       dask_cupy_func=_evi_dask_cupy)
 
-    out = mapper(red_agg)(nir_agg.data, red_agg.data, blue_agg.data, c1, c2,
-                          soil_factor, gain)
+    out = mapper(red_agg)(
+        nir_agg.data.astype('f4'), red_agg.data.astype('f4'), blue_agg.data.astype('f4'),
+        c1, c2, soil_factor, gain
+    )
 
     return DataArray(out,
                      name=name,
@@ -431,7 +435,7 @@ def gci(nir_agg: xr.DataArray,
                                       cupy_func=_gci_cupy,
                                       dask_cupy_func=_gci_dask_cupy)
 
-    out = mapper(nir_agg)(nir_agg.data, green_agg.data)
+    out = mapper(nir_agg)(nir_agg.data.astype('f4'), green_agg.data.astype('f4'))
 
     return DataArray(out,
                      name=name,
@@ -510,7 +514,7 @@ def nbr(nir_agg: xr.DataArray,
         dask_cupy_func=_run_normalized_ratio_dask_cupy,
     )
 
-    out = mapper(nir_agg)(nir_agg.data, swir2_agg.data)
+    out = mapper(nir_agg)(nir_agg.data.astype('f4'), swir2_agg.data.astype('f4'))
 
     return DataArray(out,
                      name=name,
@@ -594,7 +598,7 @@ def nbr2(swir1_agg: xr.DataArray,
         dask_cupy_func=_run_normalized_ratio_dask_cupy,
     )
 
-    out = mapper(swir1_agg)(swir1_agg.data, swir2_agg.data)
+    out = mapper(swir1_agg)(swir1_agg.data.astype('f4'), swir2_agg.data.astype('f4'))
 
     return DataArray(out,
                      name=name,
@@ -671,7 +675,7 @@ def ndvi(nir_agg: xr.DataArray,
         dask_cupy_func=_run_normalized_ratio_dask_cupy,
     )
 
-    out = mapper(nir_agg)(nir_agg.data, red_agg.data)
+    out = mapper(nir_agg)(nir_agg.data.astype('f4'), red_agg.data.astype('f4'))
 
     return DataArray(out,
                      name=name,
@@ -753,7 +757,7 @@ def ndmi(nir_agg: xr.DataArray,
         dask_cupy_func=_run_normalized_ratio_dask_cupy,
     )
 
-    out = mapper(nir_agg)(nir_agg.data, swir1_agg.data)
+    out = mapper(nir_agg)(nir_agg.data.astype('f4'), swir1_agg.data.astype('f4'))
 
     return DataArray(out,
                      name=name,
@@ -937,7 +941,7 @@ def savi(nir_agg: xr.DataArray,
                                       cupy_func=_savi_cupy,
                                       dask_cupy_func=_savi_dask_cupy)
 
-    out = mapper(red_agg)(nir_agg.data, red_agg.data, soil_factor)
+    out = mapper(red_agg)(nir_agg.data.astype('f4'), red_agg.data.astype('f4'), soil_factor)
 
     return DataArray(out,
                      name=name,
@@ -1071,7 +1075,9 @@ def sipi(nir_agg: xr.DataArray,
                                       cupy_func=_sipi_cupy,
                                       dask_cupy_func=_sipi_dask_cupy)
 
-    out = mapper(red_agg)(nir_agg.data, red_agg.data, blue_agg.data)
+    out = mapper(red_agg)(
+        nir_agg.data.astype('f4'), red_agg.data.astype('f4'), blue_agg.data.astype('f4')
+    )
 
     return DataArray(out,
                      name=name,
@@ -1238,7 +1244,9 @@ def ebbi(red_agg: xr.DataArray,
                                       cupy_func=_ebbi_cupy,
                                       dask_cupy_func=_ebbi_dask_cupy)
 
-    out = mapper(red_agg)(red_agg.data, swir_agg.data, tir_agg.data)
+    out = mapper(red_agg)(
+        red_agg.data.astype('f4'), swir_agg.data.astype('f4'), tir_agg.data.astype('f4')
+    )
 
     return DataArray(out,
                      name=name,
@@ -1298,7 +1306,7 @@ def _normalize_data(agg, pixel_max, c, th):
                                       dask_func=_normalize_data_dask,
                                       cupy_func=_normalize_data_cupy,
                                       dask_cupy_func=_normalize_data_dask_cupy)
-    out = mapper(agg)(agg.data, pixel_max, c, th)
+    out = mapper(agg)(agg.data.astype('f4'), pixel_max, c, th)
     return out
 
 

xrspatial/tests/test_multispectral.py

@@ -281,6 +281,60 @@ def qgis_ebbi():
     return result
 
 
+@pytest.fixture
+def data_uint_dtype_normalized_ratio(dtype):
+    # test data for input data array of uint dtype
+    # normalized ratio is applied with different bands for NBR, NBR2, NDVI, NDMI.
+    band1 = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    band2 = xr.DataArray(np.array([[0, 2], [1, 2]], dtype=dtype))
+    result = np.array([[1, -0.33333334], [0, -0.33333334]], dtype=np.float32)
+    return band1, band2, result
+
+
+@pytest.fixture
+def data_uint_dtype_arvi(dtype):
+    nir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    red = xr.DataArray(np.array([[0, 1], [0, 2]], dtype=dtype))
+    blue = xr.DataArray(np.array([[0, 2], [1, 2]], dtype=dtype))
+    result = np.array([[1, 0.2], [1, -0.14285715]], dtype=np.float32)
+    return nir, red, blue, result
+
+
+@pytest.fixture
+def data_uint_dtype_evi(dtype):
+    nir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    red = xr.DataArray(np.array([[0, 1], [0, 2]], dtype=dtype))
+    blue = xr.DataArray(np.array([[0, 2], [1, 2]], dtype=dtype))
+    result = np.array([[1.25, 0.], [-0.45454547, 2.5]], dtype=np.float32)
+    return nir, red, blue, result
+
+
+@pytest.fixture
+def data_uint_dtype_savi(dtype):
+    nir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    red = xr.DataArray(np.array([[0, 1], [0, 2]], dtype=dtype))
+    result = np.array([[0.25, 0.], [0.25, -0.125]], dtype=np.float32)
+    return nir, red, result
+
+
+@pytest.fixture
+def data_uint_dtype_sipi(dtype):
+    nir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    red = xr.DataArray(np.array([[0, 0], [0, 2]], dtype=dtype))
+    blue = xr.DataArray(np.array([[0, 2], [1, 2]], dtype=dtype))
+    result = np.array([[1, -1], [0, 1]], dtype=np.float32)
+    return nir, red, blue, result
+
+
+@pytest.fixture
+def data_uint_dtype_ebbi(dtype):
+    red = xr.DataArray(np.array([[0, 0], [0, 2]], dtype=dtype))
+    swir = xr.DataArray(np.array([[1, 1], [1, 1]], dtype=dtype))
+    tir = xr.DataArray(np.array([[0, 2], [1, 2]], dtype=dtype))
+    result = np.array([[0.1, 0.05773503], [0.07071068, -0.05773503]], dtype=np.float32)
+    return red, swir, tir, result
+
+
 # NDVI -------------
 def test_ndvi_data_contains_valid_values():
     _x = np.mgrid[1:0:21j]
@@ -310,6 +364,13 @@ def test_ndvi_cpu_against_qgis(nir_data, red_data, qgis_ndvi):
     general_output_checks(nir_data, result, qgis_ndvi, verify_dtype=True)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_ndvi_uint_dtype(data_uint_dtype_normalized_ratio):
+    nir_data, red_data, result_ndvi = data_uint_dtype_normalized_ratio
+    result = ndvi(nir_data, red_data)
+    general_output_checks(nir_data, result, result_ndvi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_ndvi_gpu(nir_data, red_data, qgis_ndvi):
@@ -325,6 +386,7 @@ def test_savi_zero_soil_factor_cpu_against_qgis(nir_data, red_data, qgis_ndvi):
     general_output_checks(nir_data, qgis_savi, qgis_ndvi, verify_dtype=True)
 
 
+@cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_savi_zero_soil_factor_gpu(nir_data, red_data, qgis_ndvi):
     # savi should be same as ndvi at soil_factor=0
@@ -339,6 +401,13 @@ def test_savi_cpu_against_qgis(nir_data, red_data, qgis_savi):
     general_output_checks(nir_data, result, qgis_savi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_savi_uint_dtype(data_uint_dtype_savi):
+    nir_data, red_data, result_savi = data_uint_dtype_savi
+    result = savi(nir_data, red_data)
+    general_output_checks(nir_data, result, result_savi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_savi_gpu(nir_data, red_data, qgis_savi):
@@ -354,6 +423,13 @@ def test_arvi_cpu_against_qgis(nir_data, red_data, blue_data, qgis_arvi):
     general_output_checks(nir_data, result, qgis_arvi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_arvi_uint_dtype(data_uint_dtype_arvi):
+    nir_data, red_data, blue_data, result_arvi = data_uint_dtype_arvi
+    result = arvi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_arvi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_arvi_gpu(nir_data, red_data, blue_data, qgis_arvi):
@@ -368,6 +444,13 @@ def test_evi_cpu_against_qgis(nir_data, red_data, blue_data, qgis_evi):
     general_output_checks(nir_data, result, qgis_evi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_evi_uint_dtype(data_uint_dtype_evi):
+    nir_data, red_data, blue_data, result_evi = data_uint_dtype_evi
+    result = evi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_evi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_evi_gpu(nir_data, red_data, blue_data, qgis_evi):
@@ -396,6 +479,13 @@ def test_sipi_cpu_against_qgis(nir_data, red_data, blue_data, qgis_sipi):
     general_output_checks(nir_data, result, qgis_sipi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_sipi_uint_dtype(data_uint_dtype_sipi):
+    nir_data, red_data, blue_data, result_sipi = data_uint_dtype_sipi
+    result = sipi(nir_data, red_data, blue_data)
+    general_output_checks(nir_data, result, result_sipi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_sipi_gpu(nir_data, red_data, blue_data, qgis_sipi):
@@ -410,6 +500,13 @@ def test_nbr_cpu_against_qgis(nir_data, swir2_data, qgis_nbr):
     general_output_checks(nir_data, result, qgis_nbr)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_nbr_uint_dtype(data_uint_dtype_normalized_ratio):
+    nir_data, red_data, result_nbr = data_uint_dtype_normalized_ratio
+    result = nbr(nir_data, red_data)
+    general_output_checks(nir_data, result, result_nbr, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_nbr_gpu(nir_data, swir2_data, qgis_nbr):
@@ -424,6 +521,13 @@ def test_nbr2_cpu_against_qgis(swir1_data, swir2_data, qgis_nbr2):
     general_output_checks(swir1_data, result, qgis_nbr2)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_nbr2_uint_dtype(data_uint_dtype_normalized_ratio):
+    nir_data, red_data, result_nbr2 = data_uint_dtype_normalized_ratio
+    result = nbr2(nir_data, red_data)
+    general_output_checks(nir_data, result, result_nbr2, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_nbr2_gpu(swir1_data, swir2_data, qgis_nbr2):
@@ -438,6 +542,13 @@ def test_ndmi_cpu_against_qgis(nir_data, swir1_data, qgis_ndmi):
     general_output_checks(nir_data, result, qgis_ndmi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_ndmi_uint_dtype(data_uint_dtype_normalized_ratio):
+    nir_data, red_data, result_ndmi = data_uint_dtype_normalized_ratio
+    result = ndmi(nir_data, red_data)
+    general_output_checks(nir_data, result, result_ndmi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_ndmi_gpu(nir_data, swir1_data, qgis_ndmi):
@@ -452,6 +563,13 @@ def test_ebbi_cpu_against_qgis(red_data, swir1_data, tir_data, qgis_ebbi):
     general_output_checks(red_data, result, qgis_ebbi)
 
 
+@pytest.mark.parametrize("dtype", ["uint8", "uint16"])
+def test_ebbi_uint_dtype(data_uint_dtype_ebbi):
+    red_data, swir_data, tir_data, result_ebbi = data_uint_dtype_ebbi
+    result = ebbi(red_data, swir_data, tir_data)
+    general_output_checks(red_data, result, result_ebbi, verify_dtype=True)
+
+
 @cuda_and_cupy_available
 @pytest.mark.parametrize("backend", ["cupy", "dask+cupy"])
 def test_ebbi_gpu(red_data, swir1_data, tir_data, qgis_ebbi):