Mirror geoutils API changes (#462)

Co-authored-by: Romain Hugonnet <romain.hugonnet@gmail.com>

.gitignore

@@ -153,6 +153,7 @@ examples/data/*.csv
 doc/source/basic_examples/
 doc/source/advanced_examples/
 doc/source/gen_modules/
+doc/source/sg_execution_times.rst
 examples/basic/temp.tif
 
 # Directory where myst_nb executes jupyter code

dev-environment.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python>=3.9,<3.12
-  - geopandas>=0.12.0,<0.14
+  - geopandas>=0.12.0
   - numba=0.*
   - numpy=1.*
   - matplotlib=3.*
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0
-  - geoutils=0.0.15
+  - geoutils=0.1.*
 
   # Development-specific, to mirror manually in setup.cfg [options.extras_require].
   - pip
@@ -51,4 +51,4 @@ dependencies:
     - noisyopt
 
     # To run CI against latest GeoUtils
-#    - git+https://github.com/GlacioHack/geoutils.git
+    # - git+https://github.com/GlacioHack/geoutils.git

environment.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python>=3.9,<3.12
-  - geopandas>=0.12.0,<0.14
+  - geopandas>=0.12.0
   - numba=0.*
   - numpy=1.*
   - matplotlib=3.*
@@ -13,7 +13,7 @@ dependencies:
   - tqdm
   - scikit-image=0.*
   - scikit-gstat>=1.0
-  - geoutils=0.0.15
+  - geoutils=0.1.*
   - pip
 
 #  - pip:

examples/advanced/plot_blockwise_coreg.py

@@ -48,7 +48,7 @@
 
 diff_before = reference_dem - dem_to_be_aligned
 
-diff_before.show(cmap="coolwarm_r", vmin=-10, vmax=10)
+diff_before.plot(cmap="coolwarm_r", vmin=-10, vmax=10)
 plt.show()
 
 # %%
@@ -92,7 +92,7 @@
 
 diff_after = reference_dem - aligned_dem
 
-diff_after.show(cmap="coolwarm_r", vmin=-10, vmax=10)
+diff_after.plot(cmap="coolwarm_r", vmin=-10, vmax=10)
 plt.show()
 
 # %%

examples/advanced/plot_deramp.py

@@ -26,7 +26,7 @@
 # These can be visualized by plotting a change map:
 
 diff_before = reference_dem - dem_to_be_aligned
-diff_before.show(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
+diff_before.plot(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
 
 
 # %%
@@ -41,7 +41,7 @@
 # Then, the new difference can be plotted.
 
 diff_after = reference_dem - corrected_dem
-diff_after.show(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
+diff_after.plot(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
 
 
 # %%

examples/advanced/plot_heterosc_estimation_modelling.py

@@ -269,4 +269,4 @@
 maxc = np.maximum(np.abs(profc), np.abs(planc))
 errors = dh.copy(new_array=dh_err_fun((slope.data, maxc.data)))
 
-errors.show(cmap="Reds", vmin=2, vmax=8, cbar_title=r"Elevation error ($1\sigma$, m)")
+errors.plot(cmap="Reds", vmin=2, vmax=8, cbar_title=r"Elevation error ($1\sigma$, m)")

examples/advanced/plot_norm_regional_hypso.py

@@ -61,7 +61,7 @@
 np.random.seed(42)
 random_nans = (xdem.misc.generate_random_field(dem_1990.shape, corr_size=5) > 0.7) & (glacier_index_map > 0)
 
-random_nans.show()
+random_nans.plot()
 
 # %%
 # The normalized hypsometric signal shows the tendency for elevation change as a function of elevation.

examples/advanced/plot_variogram_estimation_modelling.py

@@ -50,7 +50,7 @@
 # **Does this mean that every pixel has an independent measurement error of** :math:`\pm` **2.5 meters?**
 # Let's plot the elevation differences to visually check the quality of the data.
 plt.figure(figsize=(8, 5))
-dh.show(ax=plt.gca(), cmap="RdYlBu", vmin=-4, vmax=4, cbar_title="Elevation differences (m)")
+dh.plot(ax=plt.gca(), cmap="RdYlBu", vmin=-4, vmax=4, cbar_title="Elevation differences (m)")
 
 # %%
 # We clearly see that the residual elevation differences on stable terrain are not random. The positive and negative
@@ -66,7 +66,7 @@
 # %%
 # We plot the elevation differences after filtering to check that we successively removed glacier signals.
 plt.figure(figsize=(8, 5))
-dh.show(ax=plt.gca(), cmap="RdYlBu", vmin=-4, vmax=4, cbar_title="Elevation differences (m)")
+dh.plot(ax=plt.gca(), cmap="RdYlBu", vmin=-4, vmax=4, cbar_title="Elevation differences (m)")
 
 # %%
 # To quantify the spatial correlation of the data, we sample an empirical variogram.

examples/basic/plot_dem_subtraction.py

@@ -50,7 +50,7 @@
 # It is a new :class:`xdem.DEM` instance, loaded in memory.
 # Let's visualize it:
 
-ddem.show(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
+ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
 
 # %%
 # Let's add some glacier outlines
@@ -61,7 +61,7 @@
 # Need to create a common matplotlib Axes to plot both on the same figure
 # The xlim/ylim commands are necessary only because outlines extend further than the raster extent
 ax = plt.subplot(111)
-ddem.show(ax=ax, cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
+ddem.plot(ax=ax, cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
 glacier_outlines.ds.plot(ax=ax, fc="none", ec="k")
 plt.xlim(ddem.bounds.left, ddem.bounds.right)
 plt.ylim(ddem.bounds.bottom, ddem.bounds.top)

examples/basic/plot_icp_coregistration.py

@@ -22,11 +22,11 @@
 dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
 
 subset_extent = [523000, 8660000, 529000, 8665000]
-dem.crop(subset_extent)
+dem = dem.crop(subset_extent)
 
 # %%
 # Let's plot a hillshade of the mountain for context.
-xdem.terrain.hillshade(dem).show(cmap="gray")
+xdem.terrain.hillshade(dem).plot(cmap="gray")
 
 # %%
 # To try the effects of rotation, we can artificially rotate the DEM using a transformation matrix.
@@ -51,7 +51,7 @@
 # We can plot the difference between the original and rotated DEM.
 # It is now artificially tilting from east down to the west.
 diff_before = dem - rotated_dem
-diff_before.show(cmap="coolwarm_r", vmin=-20, vmax=20)
+diff_before.plot(cmap="coolwarm_r", vmin=-20, vmax=20)
 plt.show()
 
 # %%
@@ -83,7 +83,7 @@
 
     ax = plt.subplot(3, 1, i + 1)
     plt.title(name)
-    diff.show(cmap="coolwarm_r", vmin=-20, vmax=20, ax=ax)
+    diff.plot(cmap="coolwarm_r", vmin=-20, vmax=20, ax=ax)
 
 plt.tight_layout()
 plt.show()

examples/basic/plot_infer_heterosc.py

@@ -35,7 +35,7 @@
 
 # %%
 # The first output corresponds to the error map for the DEM (:math:`\pm` 1\ :math:`\sigma` level):
-errors.show(vmin=2, vmax=7, cmap="Reds", cbar_title=r"Elevation error (1$\sigma$, m)")
+errors.plot(vmin=2, vmax=7, cmap="Reds", cbar_title=r"Elevation error (1$\sigma$, m)")
 
 # %%
 # The second output is the dataframe of 2D binning with slope and maximum curvature:

examples/basic/plot_nuth_kaab.py

@@ -27,7 +27,7 @@
 # These can be visualized by plotting a change map:
 
 diff_before = reference_dem - dem_to_be_aligned
-diff_before.show(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
+diff_before.plot(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
 
 
 # %%
@@ -44,7 +44,7 @@
 # Then, the new difference can be plotted to validate that it improved.
 
 diff_after = reference_dem - aligned_dem
-diff_after.show(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
+diff_after.plot(cmap="coolwarm_r", vmin=-10, vmax=10, cbar_title="Elevation change (m)")
 
 
 # %%

examples/basic/plot_terrain_attributes.py

@@ -34,7 +34,7 @@ def plot_attribute(attribute, cmap, label=None, vlim=None):
     else:
         vlims = {}
 
-    attribute.show(ax=ax, cmap=cmap, add_cbar=add_cbar, cbar_title=label, **vlims)
+    attribute.plot(ax=ax, cmap=cmap, add_cbar=add_cbar, cbar_title=label, **vlims)
 
     plt.xticks([])
     plt.yticks([])

requirements.txt

@@ -1,7 +1,7 @@
 # This file is auto-generated from environment.yml, do not modify.
 # See that file for comments about the need/usage of each dependency.
 
-geopandas>=0.12.0,<0.14
+geopandas>=0.12.0
 numba==0.*
 numpy==1.*
 matplotlib==3.*
@@ -11,7 +11,7 @@ scipy==1.*
 tqdm
 scikit-image==0.*
 scikit-gstat>=1.0
-geoutils==0.0.15
+geoutils==0.1.*
 pip
 setuptools>=64
 setuptools_scm[toml]>=8

tests/test_coreg/test_affine.py

@@ -182,9 +182,9 @@ def test_coreg_example_shift(self, shift_px, coreg_class, points_or_raster, verb
 
         # shift DEM by shift_px
         shifted_ref = self.ref.copy()
-        shifted_ref.shift(shift_px[0] * res, shift_px[1] * res)
+        shifted_ref.shift(shift_px[0] * res, shift_px[1] * res, inplace=True)
 
-        shifted_ref_points = shifted_ref.to_points(as_array=False, subset=subsample, pixel_offset="center").ds
+        shifted_ref_points = shifted_ref.to_points(as_array=False, subsample=subsample, pixel_offset="center").ds
         shifted_ref_points["E"] = shifted_ref_points.geometry.x
         shifted_ref_points["N"] = shifted_ref_points.geometry.y
         shifted_ref_points.rename(columns={"b1": "z"}, inplace=True)

tests/test_coreg/test_base.py

@@ -259,7 +259,7 @@ def test_coreg_raster_and_ndarray_args(self) -> None:
         dem1.data[1, 1] = 100
 
         # Translate the DEM 1 "meter" right and add a vertical shift
-        dem2 = dem1.reproject(dst_bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
+        dem2 = dem1.reproject(bounds=rio.coords.BoundingBox(1, 0, 6, 5), silent=True)
         dem2 += 1
 
         # Create a vertical shift correction for Rasters ("_r") and for arrays ("_a")
@@ -374,7 +374,7 @@ def test_apply_resample(self, inputs: list[Any]) -> None:
                 "dem1.crs",
                 "fit",
                 "warns",
-                "'reference_dem' .* overrides the given 'transform'",
+                "'reference_dem' .* overrides the given *",
             ),
             ("dem1.data", "dem2", "dem1.transform", "None", "fit", "warns", "'dem_to_be_aligned' .* overrides .*"),
             (
@@ -635,7 +635,7 @@ def test_pipeline_pts(self) -> None:
         warnings.simplefilter("ignore")
 
         pipeline = coreg.NuthKaab() + coreg.GradientDescending()
-        ref_points = self.ref.to_points(as_array=False, subset=5000, pixel_offset="center").ds
+        ref_points = self.ref.to_points(as_array=False, subsample=5000, pixel_offset="center").ds
         ref_points["E"] = ref_points.geometry.x
         ref_points["N"] = ref_points.geometry.y
         ref_points.rename(columns={"b1": "z"}, inplace=True)
@@ -783,8 +783,8 @@ def test_blockwise_coreg(self, pipeline: Coreg, subdivision: int) -> None:
     def test_blockwise_coreg_large_gaps(self) -> None:
         """Test BlockwiseCoreg when large gaps are encountered, e.g. around the frame of a rotated DEM."""
         warnings.simplefilter("error")
-        reference_dem = self.ref.reproject(dst_crs="EPSG:3413", dst_res=self.ref.res, resampling="bilinear")
-        dem_to_be_aligned = self.tba.reproject(dst_ref=reference_dem, resampling="bilinear")
+        reference_dem = self.ref.reproject(crs="EPSG:3413", res=self.ref.res, resampling="bilinear")
+        dem_to_be_aligned = self.tba.reproject(ref=reference_dem, resampling="bilinear")
 
         blockwise = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), 64, warn_failures=False)
 

tests/test_coreg/test_biascorr.py

@@ -440,7 +440,7 @@ def test_directionalbias__synthetic(self, angle, nb_freq) -> None:
             synth = self.ref.copy(new_array=synthetic_bias.reshape(np.shape(self.ref.data)))
             import matplotlib.pyplot as plt
 
-            synth.show()
+            synth.plot()
             plt.show()
 
             dirbias = biascorr.DirectionalBias(angle=angle, fit_or_bin="bin", bin_sizes=10000)

tests/test_dem.py

@@ -98,7 +98,7 @@ def test_init__vcrs(self) -> None:
         # Tests 2: instantiation with a file that has a 3D CRS
         # Create such a file
         dem = DEM(fn_img)
-        dem_reproj = dem.reproject(dst_crs=4979)
+        dem_reproj = dem.reproject(crs=4979)
 
         # Save to temporary folder
         temp_dir = tempfile.TemporaryDirectory()
@@ -186,7 +186,7 @@ def test_copy(self) -> None:
         assert isinstance(r2, xdem.dem.DEM)
 
         # Check all immutable attributes are equal
-        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata',
+        # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'bands', 'nodata',
         #                    'res', 'shape', 'transform', 'width']
         # satimg_attrs = ['satellite', 'sensor', 'product', 'version', 'tile_name', 'datetime']
         # dem_attrs = ['vcrs', 'vcrs_grid', 'vcrs_name', 'ccrs']
@@ -263,15 +263,15 @@ def test_to_vcrs(self) -> None:
         dem = DEM(fn_dem)
 
         # Reproject in WGS84 2D
-        dem = dem.reproject(dst_crs=4326)
+        dem = dem.reproject(crs=4326)
         dem_before_trans = dem.copy()
 
         # Set ellipsoid as vertical reference
         dem.set_vcrs(new_vcrs="Ellipsoid")
         ccrs_init = dem.ccrs
         median_before = np.nanmean(dem)
         # Transform to EGM96 geoid
-        dem.to_vcrs(dst_vcrs="EGM96")
+        dem.to_vcrs(vcrs="EGM96")
         median_after = np.nanmean(dem)
 
         # About 32 meters of difference in Svalbard between EGM96 geoid and ellipsoid

xdem/dem.py

@@ -99,8 +99,8 @@ def __init__(
                 warnings.filterwarnings("ignore", message="Parse metadata from file not implemented")
                 super().__init__(filename_or_dataset, silent=silent, **kwargs)
 
-        # Ensure DEM has only one band: self.indexes can be None when data is not loaded through the Raster class
-        if self.indexes is not None and len(self.indexes) > 1:
+        # Ensure DEM has only one band: self.bands can be None when data is not loaded through the Raster class
+        if self.bands is not None and len(self.bands) > 1:
             raise ValueError("DEM rasters should be composed of one band only")
 
         # If the CRS in the raster metadata has a 3rd dimension, could set it as a vertical reference
@@ -233,39 +233,44 @@ def ccrs(self) -> CompoundCRS | CRS | None:
 
     def to_vcrs(
         self,
-        dst_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
-        src_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int | None = None,
+        vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int,
+        force_source_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"]
+        | str
+        | pathlib.Path
+        | VerticalCRS
+        | int
+        | None = None,
     ) -> None:
         """
         Convert the DEM to another vertical coordinate reference system.
 
-        :param dst_vcrs: Destination vertical CRS. Either as a name ("WGS84", "EGM08", "EGM96"),
+        :param vcrs: Destination vertical CRS. Either as a name ("WGS84", "EGM08", "EGM96"),
             an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data)
-        :param src_vcrs: Force a source vertical CRS (uses metadata by default). Same formats as for `dst_vcrs`.
+        :param force_source_vcrs: Force a source vertical CRS (uses metadata by default). Same formats as for `vcrs`.
 
         :return:
         """
 
-        if self.vcrs is None and src_vcrs is None:
+        if self.vcrs is None and force_source_vcrs is None:
             raise ValueError(
                 "The current DEM has no vertical reference, define one with .set_vref() "
                 "or by passing `src_vcrs` to perform a conversion."
             )
 
         # Initial Compound CRS (only exists if vertical CRS is not None, as checked above)
-        if src_vcrs is not None:
+        if force_source_vcrs is not None:
             # Warn if a vertical CRS already existed for that DEM
             if self.vcrs is not None:
                 warnings.warn(
                     category=UserWarning,
                     message="Overriding the vertical CRS of the DEM with the one provided in `src_vcrs`.",
                 )
-            src_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=src_vcrs)
+            src_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=force_source_vcrs)
         else:
             src_ccrs = self.ccrs
 
         # New destination Compound CRS
-        dst_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=_vcrs_from_user_input(vcrs_input=dst_vcrs))
+        dst_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=_vcrs_from_user_input(vcrs_input=vcrs))
 
         # If both compound CCRS are equal, do not run any transform
         if src_ccrs.equals(dst_ccrs):
@@ -284,4 +289,4 @@ def to_vcrs(
         self._data = zz_trans.astype(self.dtypes[0])  # type: ignore
 
         # Update vcrs (which will update ccrs if called)
-        self.set_vcrs(new_vcrs=dst_vcrs)
+        self.set_vcrs(new_vcrs=vcrs)

xdem/spatialstats.py

@@ -2402,7 +2402,7 @@ def number_effective_samples(
         elif isinstance(rasterize_resolution, Raster):
 
             # With a Raster we can get the coordinates directly
-            mask = V.create_mask(rst=rasterize_resolution, as_array=True).squeeze()
+            mask = V.create_mask(raster=rasterize_resolution, as_array=True).squeeze()
             coords = np.array(rasterize_resolution.coords())
             coords_on_mask = coords[:, mask].T