From f8ea8ac673f9dc6f712364ccf726ef2d7c625be8 Mon Sep 17 00:00:00 2001 From: Amaury Dehecq <3285905+adehecq@users.noreply.github.com> Date: Wed, 14 Feb 2024 19:42:39 +0100 Subject: [PATCH] Mirror geoutils API changes (#462) Co-authored-by: Romain Hugonnet --- .gitignore | 1 + dev-environment.yml | 6 ++--- environment.yml | 4 +-- examples/advanced/plot_blockwise_coreg.py | 4 +-- examples/advanced/plot_deramp.py | 4 +-- .../plot_heterosc_estimation_modelling.py | 2 +- examples/advanced/plot_norm_regional_hypso.py | 2 +- .../plot_variogram_estimation_modelling.py | 4 +-- examples/basic/plot_dem_subtraction.py | 4 +-- examples/basic/plot_icp_coregistration.py | 8 +++--- examples/basic/plot_infer_heterosc.py | 2 +- examples/basic/plot_nuth_kaab.py | 4 +-- examples/basic/plot_terrain_attributes.py | 2 +- requirements.txt | 4 +-- tests/test_coreg/test_affine.py | 4 +-- tests/test_coreg/test_base.py | 10 +++---- tests/test_coreg/test_biascorr.py | 2 +- tests/test_dem.py | 8 +++--- xdem/dem.py | 27 +++++++++++-------- xdem/spatialstats.py | 2 +- 20 files changed, 55 insertions(+), 49 deletions(-) diff --git a/.gitignore b/.gitignore index 8730d8dc..0882d928 100644 --- a/.gitignore +++ b/.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 diff --git a/dev-environment.yml b/dev-environment.yml index 63a046ea..7b3add3f 100644 --- a/dev-environment.yml +++ b/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 diff --git a/environment.yml b/environment.yml index 190afb5e..f74242b2 100644 --- a/environment.yml +++ b/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: diff --git a/examples/advanced/plot_blockwise_coreg.py b/examples/advanced/plot_blockwise_coreg.py index 47dfc65c..873d3716 100644 --- a/examples/advanced/plot_blockwise_coreg.py +++ b/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() # %% diff --git a/examples/advanced/plot_deramp.py b/examples/advanced/plot_deramp.py index 218c737b..20c4f205 100644 --- a/examples/advanced/plot_deramp.py +++ b/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)") # %% diff --git a/examples/advanced/plot_heterosc_estimation_modelling.py b/examples/advanced/plot_heterosc_estimation_modelling.py index ede8a95c..a979fbae 100644 --- a/examples/advanced/plot_heterosc_estimation_modelling.py +++ b/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)") diff --git a/examples/advanced/plot_norm_regional_hypso.py b/examples/advanced/plot_norm_regional_hypso.py index bf6ff0b3..c6eacada 100644 --- a/examples/advanced/plot_norm_regional_hypso.py +++ b/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. diff --git a/examples/advanced/plot_variogram_estimation_modelling.py b/examples/advanced/plot_variogram_estimation_modelling.py index 230471ec..a32eddfa 100644 --- a/examples/advanced/plot_variogram_estimation_modelling.py +++ b/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. diff --git a/examples/basic/plot_dem_subtraction.py b/examples/basic/plot_dem_subtraction.py index f27e4bdf..f5def25f 100644 --- a/examples/basic/plot_dem_subtraction.py +++ b/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) diff --git a/examples/basic/plot_icp_coregistration.py b/examples/basic/plot_icp_coregistration.py index 679a6176..2fbf28e6 100644 --- a/examples/basic/plot_icp_coregistration.py +++ b/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() diff --git a/examples/basic/plot_infer_heterosc.py b/examples/basic/plot_infer_heterosc.py index 2d1fe528..691d49ec 100644 --- a/examples/basic/plot_infer_heterosc.py +++ b/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: diff --git a/examples/basic/plot_nuth_kaab.py b/examples/basic/plot_nuth_kaab.py index 61bf86b8..e390f7af 100644 --- a/examples/basic/plot_nuth_kaab.py +++ b/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)") # %% diff --git a/examples/basic/plot_terrain_attributes.py b/examples/basic/plot_terrain_attributes.py index 50a132e3..d0b79d3d 100644 --- a/examples/basic/plot_terrain_attributes.py +++ b/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([]) diff --git a/requirements.txt b/requirements.txt index 97d485c0..82f06cdc 100644 --- a/requirements.txt +++ b/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 diff --git a/tests/test_coreg/test_affine.py b/tests/test_coreg/test_affine.py index baf5e54e..251d68c8 100644 --- a/tests/test_coreg/test_affine.py +++ b/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) diff --git a/tests/test_coreg/test_base.py b/tests/test_coreg/test_base.py index d7ff64ef..aab04f54 100644 --- a/tests/test_coreg/test_base.py +++ b/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) diff --git a/tests/test_coreg/test_biascorr.py b/tests/test_coreg/test_biascorr.py index b7a7e6b8..04d8f943 100644 --- a/tests/test_coreg/test_biascorr.py +++ b/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) diff --git a/tests/test_dem.py b/tests/test_dem.py index 0179d06d..0e4224cd 100644 --- a/tests/test_dem.py +++ b/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,7 +263,7 @@ 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 @@ -271,7 +271,7 @@ def test_to_vcrs(self) -> None: 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 diff --git a/xdem/dem.py b/xdem/dem.py index d1dca259..70fb5544 100644 --- a/xdem/dem.py +++ b/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) diff --git a/xdem/spatialstats.py b/xdem/spatialstats.py index 79f18d17..94a729b8 100644 --- a/xdem/spatialstats.py +++ b/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