Skip to content

Small adjustment for 3D CRS awareness #879

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rhugonnet wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Small adjustment for 3D CRS awareness #879

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4259361
Small changes for 3D CRS management directly in GeoUtils
rhugonnet Mar 16, 2026
8561506
Add warning category
rhugonnet Mar 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 18 additions & 4 deletions geoutils/raster/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
import xarray as xr
from affine import Affine
from packaging.version import Version
from rasterio.crs import CRS
from pyproj import CRS
from rasterio.enums import Resampling

from geoutils import profiler
Expand Down Expand Up @@ -636,7 +636,7 @@ def info(self, stats: bool = False, verbose: bool = True) -> None | str:
f"Grid size: {self.width}, {self.height}\n",
f"Number of bands: {self.count:d}\n",
f"Data types: {self.dtype}\n",
f"Coordinate system: {[self.crs.to_string() if self.crs is not None else None]}\n",
f"Coordinate system: {[CRS(self.crs).name if self.crs is not None else None]}\n",
f"Nodata value: {self.nodata}\n",
f"Pixel interpretation: {self.area_or_point}\n",
"Pixel size: {}, {}\n".format(*self.res),
Expand Down Expand Up @@ -942,11 +942,12 @@ def raster_allclose(
rtol=rtol,
)

def georeferenced_grid_equal(self: RasterType, other: RasterType) -> bool:
def georeferenced_grid_equal(self: RasterType, other: RasterType, warn_3d_crs: bool = True) -> bool:
"""
Check that raster shape, geotransform and CRS are equal.

:param other: Another raster.
:param warn_3d_crs: Whether to warn if 3D CRS differs.

:return: Whether the two objects have the same georeferenced grid.
"""
Expand All @@ -955,13 +956,26 @@ def georeferenced_grid_equal(self: RasterType, other: RasterType) -> bool:
if isinstance(other, xr.DataArray):
transform = other.rst.transform
shape = other.rst.shape
crs2d = CRS(other.rst.crs).to_2d()
crs = other.rst.crs
else:
transform = other.transform
shape = other.shape
crs2d = CRS(other.crs).to_2d()
crs = other.crs

return all([self.shape == shape, self.transform == transform, self.crs == crs])
# 2D CRS for self
self_crs2d = CRS(self.crs).to_2d()

gge = all([self.shape == shape, self.transform == transform, self_crs2d == crs2d])

# The only way they differ is if the vertical CRS is different
if gge and crs != self.crs:
warnings.warn("The two rasters have the same 2D CRS but a different vertical CRS: "
f"{CRS(self.crs).name} and {CRS(crs).name}.",
category=UserWarning)

return gge

def get_bounds_projected(self, out_crs: CRS, densify_points: int = 5000) -> rio.coords.BoundingBox:
"""
Expand Down
Loading