Update with new GeoUtils (casting single-band to 2D, implicit loading…

…, new module names) (#341)

* Fix satimg import with new geoutils

* Make NMAD applicable to Raster objects

* Fix BlockWise bug when cadrant full of nodatas

* Incremental commit on fixing behaviour with new Mask class

* Fix behaviour of several functions with GeoUtils update

* Refactor with module renaming in geoutils

* Update with latest changes in GeoUtils, and shorten modules path for function imports

* Fix behaviour with 2D .data for DEMs

* Linting

README.md

@@ -77,7 +77,7 @@ Please see our [contribution page](CONTRIBUTING.md) for more detailed instructio
 See the documentation at https://xdem.readthedocs.io
 
 ### Testing - again please read!
-These tools are only valuable if we can rely on them to perform exactly as we expect. So, we need testing. Please create tests for every function that you make, as much as you are able. Guidance/examples here for the moment: https://github.com/GeoUtils/georaster/blob/master/test/test_georaster.py
+These tools are only valuable if we can rely on them to perform exactly as we expect. So, we need testing. Please create tests for every function that you make, as much as you are able. Guidance/examples here for the moment: https://github.com/GeoUtils/raster/blob/master/test/test_raster.py
 https://github.com/corteva/rioxarray/blob/master/test/integration/test_integration__io.py
 
 

dev-environment.yml

@@ -37,5 +37,5 @@ dependencies:
   - richdem
 
   - pip:
-    # - git+https://github.com/GlacioHack/GeoUtils.git
-    - geoutils==0.0.10
+     - git+https://github.com/GlacioHack/GeoUtils.git
+#    - geoutils==0.0.10

docs/source/code/spatialstats.py

@@ -5,20 +5,20 @@
 import xdem
 
 # Load data
-dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+dh = gu.Raster(xdem.examples.get_path("longyearbyen_ddem"))
 ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
-glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+glacier_mask = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 mask = glacier_mask.create_mask(dh)
 
 slope = xdem.terrain.get_terrain_attribute(ref_dem, attribute=["slope"])
 
 # Keep only stable terrain data
 dh.set_mask(mask)
-dh_arr = gu.spatial_tools.get_array_and_mask(dh)[0]
-slope_arr = gu.spatial_tools.get_array_and_mask(slope)[0]
+dh_arr = gu.raster.get_array_and_mask(dh)[0]
+slope_arr = gu.raster.get_array_and_mask(slope)[0]
 
 # Subsample to run the snipped code faster
-indices = gu.spatial_tools.subsample_raster(dh_arr, subsample=10000, return_indices=True, random_state=42)
+indices = gu.raster.subsample_array(dh_arr, subsample=10000, return_indices=True, random_state=42)
 dh_arr = dh_arr[indices]
 slope_arr = slope_arr[indices]
 

docs/source/code/spatialstats_heterosc_slope.py

@@ -5,9 +5,9 @@
 import xdem
 
 # Load data
-dh = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ddem"))
+dh = gu.Raster(xdem.examples.get_path("longyearbyen_ddem"))
 ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
-glacier_mask = gu.geovector.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
+glacier_mask = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
 mask = glacier_mask.create_mask(dh)
 
 # Get slope for non-stationarity

environment.yml

@@ -23,5 +23,5 @@ dependencies:
   - pip
 
   - pip:
-    # - git+https://github.com/GlacioHack/GeoUtils.git
-    - geoutils==0.0.10
+     - git+https://github.com/GlacioHack/GeoUtils.git
+#    - geoutils==0.0.10

examples/advanced/plot_heterosc_estimation_modelling.py

@@ -44,11 +44,11 @@
 
 # %%
 # We convert to arrays and keep only stable terrain for the analysis of variability
-dh_arr = gu.spatial_tools.get_array_and_mask(dh)[0][~mask_glacier]
-slope_arr = gu.spatial_tools.get_array_and_mask(slope)[0][~mask_glacier]
-aspect_arr = gu.spatial_tools.get_array_and_mask(aspect)[0][~mask_glacier]
-planc_arr = gu.spatial_tools.get_array_and_mask(planc)[0][~mask_glacier]
-profc_arr = gu.spatial_tools.get_array_and_mask(profc)[0][~mask_glacier]
+dh_arr = dh.get_nanarray()[~mask_glacier]
+slope_arr = slope.get_nanarray()[~mask_glacier]
+aspect_arr = aspect.get_nanarray()[~mask_glacier]
+planc_arr = planc.get_nanarray()[~mask_glacier]
+profc_arr = profc.get_nanarray()[~mask_glacier]
 
 # %%
 # We use :func:`xdem.spatialstats.nd_binning` to perform N-dimensional binning on all those terrain variables, with uniform

examples/basic/plot_dem_subtraction.py

@@ -38,7 +38,7 @@
 # To hide this prompt, add ``.reproject(..., silent=True)``.
 # By default, :func:`xdem.DEM.reproject` uses "bilinear" resampling (assuming resampling is needed).
 # Other options are "nearest" (fast but inaccurate), "cubic_spline", "lanczos" and others.
-# See `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.georaster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.
+# See `geoutils.Raster.reproject() <https://geoutils.readthedocs.io/en/latest/api.html#geoutils.raster.Raster.reproject>`_ and `rasterio.enums.Resampling <https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Resampling>`_ for more information about reprojection.
 #
 # Now, we are ready to generate the dDEM:
 

tests/test_coreg.py

@@ -14,7 +14,7 @@
 import pytest
 import rasterio as rio
 from geoutils import Raster, Vector
-from geoutils.georaster.raster import RasterType
+from geoutils.raster import RasterType
 
 with warnings.catch_warnings():
     warnings.simplefilter("ignore")
@@ -264,12 +264,12 @@ def test_deramping(self) -> None:
         deramp.fit(**self.fit_params)
 
         # Apply the deramping to a DEM
-        deramped_dem, _ = deramp.apply(self.tba.data, self.ref.transform, self.ref.crs)
+        deramped_dem = deramp.apply(self.tba)
 
         # Get the periglacial offset after deramping
-        periglacial_offset = (self.ref.data.squeeze() - deramped_dem)[self.inlier_mask.squeeze()]
+        periglacial_offset = (self.ref - deramped_dem)[self.inlier_mask]
         # Get the periglacial offset before deramping
-        pre_offset = ((self.ref.data - self.tba.data)[self.inlier_mask]).squeeze()
+        pre_offset = (self.ref - self.tba)[self.inlier_mask]
 
         # Check that the error improved
         assert np.abs(np.mean(periglacial_offset)) < np.abs(np.mean(pre_offset))
@@ -434,15 +434,15 @@ def test_z_scale_corr(self) -> None:
         assert np.abs(np.nanmedian(diff)) < 0.01
 
         # Create a spatially correlated error field to mess with the algorithm a bit.
-        corr_size = int(self.ref.data.shape[2] / 100)
+        corr_size = int(self.ref.data.shape[1] / 100)
         error_field = cv2.resize(
             cv2.GaussianBlur(
                 np.repeat(
                     np.repeat(
                         np.random.randint(
                             0,
                             255,
-                            (self.ref.data.shape[1] // corr_size, self.ref.data.shape[2] // corr_size),
+                            (self.ref.data.shape[0] // corr_size, self.ref.data.shape[1] // corr_size),
                             dtype="uint8",
                         ),
                         corr_size,
@@ -455,7 +455,7 @@ def test_z_scale_corr(self) -> None:
                 sigmaX=corr_size,
             )
             / 255,
-            dsize=(self.ref.data.shape[2], self.ref.data.shape[1]),
+            dsize=(self.ref.data.shape[1], self.ref.data.shape[0]),
         )
 
         # Create 50000 random nans
@@ -517,13 +517,13 @@ def test_blockwise_coreg(self, pipeline: coreg.Coreg, subdivision: int) -> None:
         chunk_numbers = [m["i"] for m in blockwise._meta["coreg_meta"]]
         assert np.unique(chunk_numbers).shape[0] == len(chunk_numbers)
 
-        transformed_dem, _ = blockwise.apply(self.tba.data, self.tba.transform, self.tba.crs)
+        transformed_dem = blockwise.apply(self.tba)
 
-        ddem_pre = (self.ref.data - self.tba.data)[~self.inlier_mask].squeeze().filled(np.nan)
-        ddem_post = (self.ref.data.squeeze() - transformed_dem)[~self.inlier_mask.squeeze()].filled(np.nan)
+        ddem_pre = (self.ref - self.tba)[~self.inlier_mask]
+        ddem_post = (self.ref - transformed_dem)[~self.inlier_mask]
 
         # Check that the periglacial difference is lower after coregistration.
-        assert abs(np.nanmedian(ddem_post)) < abs(np.nanmedian(ddem_pre))
+        assert abs(np.ma.median(ddem_post)) < abs(np.ma.median(ddem_pre))
 
         stats = blockwise.stats()
 
@@ -554,7 +554,7 @@ def test_blockwise_coreg_large_gaps(self) -> None:
         # Copy the TBA DEM and set a square portion to nodata
         tba = self.tba.copy()
         mask = np.zeros(np.shape(tba.data), dtype=bool)
-        mask[0, 450:500, 450:500] = True
+        mask[450:500, 450:500] = True
         tba.set_mask(mask=mask)
 
         blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), 8, warn_failures=False)
@@ -583,7 +583,7 @@ def test_coreg_raster_and_ndarray_args(self) -> None:
             nodata=-9999,
         )
         # Assign a funny value to one particular pixel. This is to validate that reprojection works perfectly.
-        dem1.data[0, 1, 1] = 100
+        dem1.data[1, 1] = 100
 
         # Translate the DEM 1 "meter" right and add a bias
         dem2 = dem1.reproject(dst_bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
@@ -822,7 +822,7 @@ def test_coreg_oneliner(self) -> None:
 def test_apply_matrix() -> None:
     warnings.simplefilter("error")
     ref, tba, outlines = load_examples()  # Load example reference, to-be-aligned and mask.
-    ref_arr = gu.spatial_tools.get_array_and_mask(ref)[0]
+    ref_arr = gu.raster.get_array_and_mask(ref)[0]
 
     # Test only bias (it should just apply the bias and not make anything else)
     bias = 5
@@ -1043,7 +1043,7 @@ def test_create_inlier_mask() -> None:
 
     # - Assert that without filtering create_inlier_mask behaves as if calling Vector.create_mask - #
     # Masking inside - using Vector
-    inlier_mask_comp = ~outlines.create_mask(ref)
+    inlier_mask_comp = ~outlines.create_mask(ref, as_array=True)
     inlier_mask = xdem.coreg.create_inlier_mask(
         tba,
         ref,
@@ -1253,7 +1253,7 @@ def test_dem_coregistration() -> None:
         ],
         resample=True,
     )
-    gl_mask = outlines.create_mask(dem_coreg)
+    gl_mask = outlines.create_mask(dem_coreg, as_array=True)
     assert np.all(~inlier_mask[gl_mask])
 
     # Testing with plot

tests/test_ddem.py

@@ -29,7 +29,7 @@ def test_copy(self) -> None:
 
         ddem2.data += 1
 
-        assert self.ddem != ddem2
+        assert not self.ddem.raster_equal(ddem2)
 
     def test_filled_data(self) -> None:
         """Test that the filled_data property points to the right data."""

tests/test_dem.py

@@ -2,13 +2,13 @@
 import os
 import warnings
 
-import geoutils.georaster as gr
-import geoutils.satimg as si
+import geoutils.raster as gr
+import geoutils.raster.satimg as si
 import numpy as np
 import pyproj
 import pytest
 import rasterio as rio
-from geoutils.georaster.raster import _default_rio_attrs
+from geoutils.raster.raster import _default_rio_attrs
 
 import xdem
 from xdem.dem import DEM
@@ -94,7 +94,7 @@ def test_copy(self) -> None:
         assert isinstance(r2, xdem.dem.DEM)
 
         # Check all immutable attributes are equal
-        # georaster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
+        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
         #                    'res', 'shape', 'transform', 'width']
         # satimg_attrs = ['satellite', 'sensor', 'product', 'version', 'tile_name', 'datetime']
         # dem_attrs = ['vref', 'vref_grid', 'ccrs']

tests/test_demcollection.py

@@ -24,10 +24,10 @@ def test_init(self) -> None:
         # Make sure the glacier was bigger in 1990, since this is assumed later.
         assert scott_1990.ds.area.sum() > scott_2010.ds.area.sum()
 
-        mask_2010 = scott_2010.create_mask(self.dem_2009).reshape(self.dem_2009.data.shape)
+        mask_2010 = scott_2010.create_mask(self.dem_2009)
 
         dem_2060 = self.dem_2009.copy()
-        dem_2060.data[mask_2010] -= 30
+        dem_2060[mask_2010] -= 30
 
         dems = xdem.DEMCollection(
             [self.dem_1990, self.dem_2009, dem_2060],
@@ -36,7 +36,7 @@ def test_init(self) -> None:
             reference_dem=1,
         )
 
-        # Check that the first raster is the oldest one and
+        # Check that the first raster is the oldest one
         assert dems.dems[0].data.max() == self.dem_1990.data.max()
         assert dems.reference_dem.data.max() == self.dem_2009.data.max()
 

tests/test_examples.py

@@ -62,6 +62,6 @@ def test_array_nodata(self, rst_and_truenodata: tuple[Raster, int]) -> None:
 
         rst = rst_and_truenodata[0]
         truenodata = rst_and_truenodata[1]
-        mask = gu.spatial_tools.get_array_and_mask(rst)[1]
+        mask = gu.raster.get_array_and_mask(rst)[1]
 
         assert np.sum(mask) == truenodata

tests/test_filters.py

@@ -12,8 +12,8 @@ class TestFilters:
     """Test cases for the filter functions."""
 
     # Load example data.
-    dem_2009 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
-    dem_1990 = gu.georaster.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
+    dem_2009 = gu.Raster(xdem.examples.get_path("longyearbyen_ref_dem"))
+    dem_1990 = gu.Raster(xdem.examples.get_path("longyearbyen_tba_dem")).reproject(dem_2009, silent=True)
 
     def test_gauss(self) -> None:
         """Test applying the various Gaussian filters on DEMs with/without NaNs"""
@@ -51,15 +51,15 @@ def test_gauss(self) -> None:
         assert np.nanmax(dem_with_nans) > np.max(dem_sm)
 
         # Test that it works with 3D arrays
-        array_3d = np.vstack((dem_array, dem_array))
+        array_3d = np.vstack((dem_array[np.newaxis, :], dem_array[np.newaxis, :]))
         dem_sm = xdem.filters.gaussian_filter_scipy(array_3d, sigma=5)
         assert array_3d.shape == dem_sm.shape
 
         # 3D case not implemented with OpenCV
         pytest.raises(NotImplementedError, xdem.filters.gaussian_filter_cv, array_3d, sigma=5)
 
         # Tests that it fails with 1D arrays with appropriate error
-        data = dem_array[0, :, 0]
+        data = dem_array[:, 0]
         pytest.raises(ValueError, xdem.filters.gaussian_filter_scipy, data, sigma=5)
         pytest.raises(ValueError, xdem.filters.gaussian_filter_cv, data, sigma=5)
 
@@ -73,19 +73,19 @@ def test_dist_filter(self) -> None:
         count = 1000
         cols = np.random.randint(0, high=self.dem_1990.width - 1, size=count, dtype=int)
         rows = np.random.randint(0, high=self.dem_1990.height - 1, size=count, dtype=int)
-        ddem.data[0, rows, cols] = 5000
+        ddem.data[rows, cols] = 5000
 
         # Filter gross outliers
         filtered_ddem = xdem.filters.distance_filter(ddem.data, radius=20, outlier_threshold=50)
 
         # Check that all outliers were properly filtered
-        assert np.all(np.isnan(filtered_ddem[0, rows, cols]))
+        assert np.all(np.isnan(filtered_ddem[rows, cols]))
 
         # Assert that non filtered pixels remain the same
         assert ddem.data.shape == filtered_ddem.shape
         assert np.all(ddem.data[np.isfinite(filtered_ddem)] == filtered_ddem[np.isfinite(filtered_ddem)])
 
         # Check that it works with NaNs too
-        ddem.data[0, rows[:500], cols[:500]] = np.nan
+        ddem.data[rows[:500], cols[:500]] = np.nan
         filtered_ddem = xdem.filters.distance_filter(ddem.data, radius=20, outlier_threshold=50)
-        assert np.all(np.isnan(filtered_ddem[0, rows, cols]))
+        assert np.all(np.isnan(filtered_ddem[rows, cols]))

tests/test_misc.py

@@ -105,3 +105,54 @@ def useless_func() -> int:
         else:
             with pytest.raises(ValueError, match="^" + text + "$"):
                 useless_func()
+
+    def test_diff_environment_yml(self, capsys) -> None:  # type: ignore
+
+        # Test with synthetic environment
+        env = {"dependencies": ["python==3.9", "numpy", "fiona"]}
+        devenv = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv"]}
+
+        # This should print the difference between the two
+        xdem.misc.diff_environment_yml(env, devenv, input_dict=True, print_dep="conda")
+
+        # Capture the stdout and check it is indeed the right diff
+        captured = capsys.readouterr().out
+        assert captured == "opencv\n"
+
+        # This should print the difference including pip
+        xdem.misc.diff_environment_yml(env, devenv, input_dict=True, print_dep="both")
+
+        captured = capsys.readouterr().out
+        assert captured == "opencv\nNone\n"
+
+        env2 = {"dependencies": ["python==3.9", "numpy", "fiona"]}
+        devenv2 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils", "-e ./"]}]}
+
+        # The diff function should not account for -e ./ that is the local install for developers
+        xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="both")
+        captured = capsys.readouterr().out
+
+        assert captured == "opencv\ngeoutils\n"
+
+        # This should print only pip
+        xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="pip")
+        captured = capsys.readouterr().out
+
+        assert captured == "geoutils\n"
+
+        # This should raise an error because print_dep is not well defined
+        with pytest.raises(ValueError, match='The argument "print_dep" can only be "conda", "pip" or "both".'):
+            xdem.misc.diff_environment_yml(env2, devenv2, input_dict=True, print_dep="lol")
+
+        # When the dependencies are not defined in dev-env but in env, it should raise an error
+        # For normal dependencies
+        env3 = {"dependencies": ["python==3.9", "numpy", "fiona", "lol"]}
+        devenv3 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils"]}]}
+        with pytest.raises(ValueError, match="The following dependencies are listed in env but not dev-env: lol"):
+            xdem.misc.diff_environment_yml(env3, devenv3, input_dict=True, print_dep="pip")
+
+        # For pip dependencies
+        env4 = {"dependencies": ["python==3.9", "numpy", "fiona", {"pip": ["lol"]}]}
+        devenv4 = {"dependencies": ["python==3.9", "numpy", "fiona", "opencv", {"pip": ["geoutils"]}]}
+        with pytest.raises(ValueError, match="The following pip dependencies are listed in env but not dev-env: lol"):
+            xdem.misc.diff_environment_yml(env4, devenv4, input_dict=True, print_dep="pip")