Test against QGIS, GDAL (#744)

Add QA tests against QGIS / GDAL

Author: thuydotm

xrspatial/tests/conftest.py

@@ -8,3 +8,48 @@ def random_data(size, dtype):
     data = rng.integers(-100, 100, size=size)
     data = data.astype(dtype)
     return data
+
+
+@pytest.fixture
+def elevation_raster():
+    elevation = np.array([
+        [   np.nan,    np.nan,    np.nan,    np.nan,    np.nan,    np.nan],
+        [704.237  , 242.24084, 429.3324 , 779.8816 , 193.29506, 984.6926 ],
+        [226.56795, 815.7483 , 290.6041 ,  76.49687, 820.89716,  32.27882],
+        [344.8238 , 256.34998, 806.8326 , 602.0442 , 721.1633 , 496.95636],
+        [185.43515, 834.10425, 387.0871 , 716.0262 ,  49.61273, 752.95483],
+        [302.4271 , 151.49211, 442.32797, 358.4702 , 659.8187 , 447.1241 ],
+        [148.04834, 819.2133 , 468.97913, 977.11694, 597.69666, 999.14185],
+        [268.1575 , 625.96466, 840.26483, 448.28333, 859.2699 , 528.04095]
+    ], dtype=np.float32)
+    return elevation
+
+
+@pytest.fixture
+def elevation_raster_no_nans():
+    elevation = np.array([
+        [870.5345 , 283.04907, 845.2779 ,  51.21859, 990.8278 , 600.64545],
+        [704.237  , 242.24084, 429.3324 , 779.8816 , 193.29506, 984.6926 ],
+        [226.56795, 815.7483 , 290.6041 ,  76.49687, 820.89716,  32.27882],
+        [344.8238 , 256.34998, 806.8326 , 602.0442 , 721.1633 , 496.95636],
+        [185.43515, 834.10425, 387.0871 , 716.0262 ,  49.61273, 752.95483],
+        [302.4271 , 151.49211, 442.32797, 358.4702 , 659.8187 , 447.1241 ],
+        [148.04834, 819.2133 , 468.97913, 977.11694, 597.69666, 999.14185],
+        [268.1575 , 625.96466, 840.26483, 448.28333, 859.2699 , 528.04095]
+    ], dtype=np.float32)
+    return elevation
+
+
+@pytest.fixture
+def raster():
+    data = np.array([
+        [6., 7., 3., 4., 8., 1.],
+        [4., 9., 7., 5., 6., 9.],
+        [4., 3., 3., 1., 3., 7.],
+        [3., 4., 9., 3., 7., 0.],
+        [2., 1., 6., 5., 6., 2.],
+        [4., 2., 4., 3., 8., 5.],
+        [4., 1., 8., 5., 7., 0.],
+        [7., 4., 6., 4., 1., 1.]
+    ], dtype=np.float32)
+    return data

xrspatial/tests/general_checks.py

@@ -11,12 +11,26 @@
 
 
 def create_test_raster(
-    data, backend='numpy', name='myraster', dims=['y', 'x'], attrs=None, chunks=(3, 3)
+        data,
+        backend='numpy',
+        name='myraster',
+        dims=['y', 'x'],
+        attrs={'res': (0.5, 0.5), 'crs': 'EPSG: 4326'},
+        chunks=(3, 3)
 ):
     raster = xr.DataArray(data, name=name, dims=dims, attrs=attrs)
-    # set coords for test raster
-    for i, dim in enumerate(dims):
-        raster[dim] = np.linspace(0, data.shape[i] - 1, data.shape[i])
+
+    # default res if none provided
+    res = (0.5, 0.5)
+    if attrs is not None:
+        if 'res' in attrs:
+            res = attrs['res']
+    # set coords for test raster, 2D coords only
+    raster[dims[0]] = np.linspace((data.shape[0] - 1) * res[0], 0, data.shape[0])
+    raster[dims[1]] = np.linspace(0, (data.shape[1] - 1) * res[1], data.shape[1])
+
+    raster[dims[0]] = np.linspace((data.shape[0] - 1)/2, 0, data.shape[0])
+    raster[dims[1]] = np.linspace(0, (data.shape[1] - 1)/2, data.shape[1])
 
     if has_cuda_and_cupy() and 'cupy' in backend:
         import cupy

xrspatial/tests/test_aspect.py

@@ -7,56 +7,46 @@
                                             cuda_and_cupy_available, general_output_checks)
 
 
-def input_data(backend='numpy'):
-    data = np.asarray([
-        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-        [1584.8767, 1584.8767, 1585.0546, 1585.2324, 1585.2324, 1585.2324],
-        [1585.0546, 1585.0546, 1585.2324, 1585.588, 1585.588, 1585.588],
-        [1585.2324, 1585.4102, 1585.588, 1585.588, 1585.588, 1585.588],
-        [1585.588, 1585.588, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
-        [1585.7659, 1585.9437, 1585.7659, 1585.7659, 1585.7659, 1585.7659],
-        [1585.9437, 1585.9437, 1585.9437, 1585.7659, 1585.7659, 1585.7659]],
-        dtype=np.float32
-    )
-    raster = create_test_raster(data, backend, attrs={'res': (10.0, 10.0)})
+def input_data(data, backend='numpy'):
+    raster = create_test_raster(data, backend)
     return raster
 
 
 @pytest.fixture
-def qgis_output():
+def qgis_aspect():
     result = np.array([
-        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-        [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-        [330.94687, 335.55496, 320.70786, 330.94464, 0., 0.],
-        [333.43494, 333.43494, 329.03394, 341.56897, 0., 18.434948],
-        [338.9621, 338.20062, 341.56506, 0., 0., 45.],
-        [341.56506, 351.8699, 26.56505, 45., -1., 90.],
-        [351.86676, 11.306906, 45., 45., 45., 108.431015]], dtype=np.float32
-    )
+        [    np.nan,     np.nan,     np.nan,     np.nan,     np.nan,    np.nan],
+        [    np.nan,     np.nan,     np.nan,     np.nan,     np.nan,    np.nan],
+        [233.19478 , 278.358   ,  45.18813 , 306.6476  , 358.34296 , 106.45898 ],
+        [267.7002  , 274.42487 ,  11.035832, 357.9641  , 129.98279 , 50.069843],
+        [263.18484 , 238.47426 , 196.37103 , 149.25227 , 187.85748 , 263.684   ],
+        [266.63937 , 271.05124 , 312.09726 , 348.89136 , 351.618   , 315.59424 ],
+        [279.90872 , 314.11356 , 345.76315 , 327.5568  , 339.5455  , 312.9249  ],
+        [271.93985 , 268.81046 ,  24.793104, 185.978   , 299.82904 ,159.0188  ]], dtype=np.float32)
     return result
 
 
-def test_numpy_equals_qgis(qgis_output):
-    numpy_agg = input_data()
+def test_numpy_equals_qgis(elevation_raster, qgis_aspect):
+    numpy_agg = input_data(elevation_raster, backend='numpy')
     xrspatial_aspect = aspect(numpy_agg, name='numpy_aspect')
 
     general_output_checks(numpy_agg, xrspatial_aspect, verify_dtype=True)
     assert xrspatial_aspect.name == 'numpy_aspect'
 
     xrspatial_vals = xrspatial_aspect.data[1:-1, 1:-1]
-    qgis_vals = qgis_output[1:-1, 1:-1]
+    qgis_vals = qgis_aspect[1:-1, 1:-1]
     # aspect is nan if nan input
     # aspect is invalid (-1) if slope equals 0
     # otherwise aspect are from 0 to 360
-    np.testing.assert_allclose(xrspatial_vals, qgis_vals, equal_nan=True)
+    np.testing.assert_allclose(xrspatial_vals, qgis_vals, rtol=1e-05, equal_nan=True)
     # nan edge effect
     assert_nan_edges_effect(xrspatial_aspect)
 
 
-def test_numpy_equals_dask_qgis_data():
+def test_numpy_equals_dask_qgis_data(elevation_raster):
     # compare using the data run through QGIS
-    numpy_agg = input_data('numpy')
-    dask_agg = input_data('dask+numpy')
+    numpy_agg = input_data(elevation_raster, 'numpy')
+    dask_agg = input_data(elevation_raster, 'dask+numpy')
     assert_numpy_equals_dask_numpy(numpy_agg, dask_agg, aspect)
 
 
@@ -72,8 +62,8 @@ def test_numpy_equals_dask_random_data(random_data):
 @cuda_and_cupy_available
 def test_numpy_equals_cupy_qgis_data():
     # compare using the data run through QGIS
-    numpy_agg = input_data()
-    cupy_agg = input_data('cupy')
+    numpy_agg = input_data(elevation_raster)
+    cupy_agg = input_data(elevation_raster, 'cupy')
     assert_numpy_equals_cupy(numpy_agg, cupy_agg, aspect)
 
 

xrspatial/tests/test_classify.py

@@ -14,7 +14,7 @@ def input_data(backend='numpy'):
         [10., 11., 12., 13., 14.],
         [15., 16., 17., 18., np.inf],
     ])
-    raster = create_test_raster(elevation, backend, attrs={'res': (10.0, 10.0)})
+    raster = create_test_raster(elevation, backend)
     return raster
 
 
@@ -218,7 +218,7 @@ def test_natural_breaks_numpy_num_sample(result_natural_breaks_num_sample):
 def test_natural_breaks_cpu_deterministic():
     results = []
     elevation = np.arange(100).reshape(10, 10)
-    agg = xr.DataArray(elevation, attrs={'res': (10.0, 10.0)})
+    agg = xr.DataArray(elevation)
 
     k = 5
     numIters = 3

xrspatial/tests/test_focal.py

@@ -270,15 +270,15 @@ def test_2d_convolution_gpu(
 
 
 def test_calc_cellsize_unit_input_attrs(convolve_2d_data):
-    agg = create_test_raster(convolve_2d_data, attrs={'res': 1, 'unit': 'km'})
+    agg = create_test_raster(convolve_2d_data, attrs={'res': (1, 1), 'unit': 'km'})
     cellsize = calc_cellsize(agg)
     assert cellsize == (1000, 1000)
 
 
 def test_calc_cellsize_no_attrs(convolve_2d_data):
     agg = create_test_raster(convolve_2d_data)
     cellsize = calc_cellsize(agg)
-    assert cellsize == (1, 1)
+    assert cellsize == (0.5, 0.5)
 
 
 @pytest.fixture