From 931c9dc540f21394e708dc2ec1fd7b2c2d39d649 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Thu, 19 Jan 2023 21:27:13 -0900 Subject: [PATCH] Add `proximity` and `buffer_metric` functionalities with `metric` option based on finding local UTM (#336) * Add latlon_to_utm and utm_to_epsg functions * Add raster proximity function * Add vector proximity, buffer_metric functions, and add metric argument to buffer_without_overlap * Add test for utm_to_epsg * Fix circular import * Add tests for latlon_to_utm * Add test for buffer_metric * Update tests for buffer_without_overlap * Linting * Make proximity function more modular * Add more parameters to proximity * Incremental commit on proximity tests * Finalize proximity tests * Linting * Fix description of buffer metric * Incremental commit on eriks comments * Finalize accounting for comments * Fix code broken when accounting for comments * Linting * Linting --- geoutils/georaster/raster.py | 127 +++++++++++++++++++++- geoutils/geovector.py | 129 +++++++++++++++++++++-- geoutils/projtools.py | 47 +++++++++ tests/test_georaster.py | 198 +++++++++++++++++++++++++++++++++++ tests/test_geovector.py | 99 +++++++++++++++++- tests/test_projtools.py | 70 +++++++++++++ 6 files changed, 661 insertions(+), 9 deletions(-) diff --git a/geoutils/georaster/raster.py b/geoutils/georaster/raster.py index fc787af7..10aac780 100644 --- a/geoutils/georaster/raster.py +++ b/geoutils/georaster/raster.py @@ -28,7 +28,7 @@ from rasterio.features import shapes from rasterio.plot import show as rshow from rasterio.warp import Resampling -from scipy.ndimage import map_coordinates +from scipy.ndimage import distance_transform_edt, map_coordinates import geoutils.geovector as gv from geoutils._typing import AnyNumber, ArrayLike, DTypeLike @@ -1146,6 +1146,17 @@ def copy(self: RasterType, new_array: np.ndarray | None = None) -> RasterType: return cp + def equal_georeferenced_grid(self: RasterType, raster: RasterType) -> bool: + """ + Check that grid shape, geotransform and CRS are equal. + + :param raster: Another Raster object + + :return: Whether the two objects have the same georeferenced grid + """ + + return all([self.shape == raster.shape, self.transform == raster.transform, self.crs == raster.crs]) + @overload def get_nanarray(self, return_mask: Literal[False] = False) -> np.ndarray: ... @@ -2510,3 +2521,117 @@ def polygonize( gdf.set_crs(self.crs, inplace=True) return gv.Vector(gdf) + + def proximity( + self, + vector: Vector | None = None, + target_values: list[float] | None = None, + geometry_type: str = "boundary", + in_or_out: Literal["in"] | Literal["out"] | Literal["both"] = "both", + distance_unit: Literal["pixel"] | Literal["georeferenced"] = "georeferenced", + ) -> Raster: + """ + Proximity to this Raster's target pixels, or to a Vector's geometry, computed for each cell of this Raster's + grid. + + When passing a Vector, by default, the boundary of the geometry will be used. The full geometry can be used by + passing "geometry", or any lower dimensional geometry attribute such as "centroid", "envelope" or "convex_hull". + See all geometry attributes in the Shapely documentation at https://shapely.readthedocs.io/. + + :param vector: Vector for which to compute the proximity to geometry, + if not provided computed on this Raster target pixels. + :param target_values: (Only with Raster) List of target values to use for the proximity, + defaults to all non-zero values. + :param geometry_type: (Only with a Vector) Type of geometry to use for the proximity, defaults to 'boundary'. + :param in_or_out: (Only with a Vector) Compute proximity only 'in' or 'out'-side the geometry, or 'both'. + :param distance_unit: Distance unit, either 'georeferenced' or 'pixel'. + + :return: Proximity raster. + """ + + proximity = proximity_from_vector_or_raster( + raster=self, + vector=vector, + target_values=target_values, + geometry_type=geometry_type, + in_or_out=in_or_out, + distance_unit=distance_unit, + ) + + return self.copy(new_array=proximity) + + +# ----------------------------------------- +# Additional stand-alone utility functions +# ----------------------------------------- + + +def proximity_from_vector_or_raster( + raster: Raster, + vector: Vector | None = None, + target_values: list[float] | None = None, + geometry_type: str = "boundary", + in_or_out: Literal["in"] | Literal["out"] | Literal["both"] = "both", + distance_unit: Literal["pixel"] | Literal["georeferenced"] = "georeferenced", +) -> np.ndarray: + """ + (This function is defined here as mostly raster-based, but used in a class method for both Raster and Vector) + Proximity to a Raster's target values if no Vector is provided, otherwise to a Vector's geometry type + rasterized on the Raster. + + :param raster: Raster to burn the proximity grid on. + :param vector: Vector for which to compute the proximity to geometry, + if not provided computed on the Raster target pixels. + :param target_values: (Only with a Raster) List of target values to use for the proximity, + defaults to all non-zero values. + :param geometry_type: (Only with a Vector) Type of geometry to use for the proximity, defaults to 'boundary'. + :param in_or_out: (Only with a Vector) Compute proximity only 'in' or 'out'-side the geometry, or 'both'. + :param distance_unit: Distance unit, either 'georeferenced' or 'pixel'. + """ + + # 1/ First, if there is a vector input, we rasterize the geometry type + # (works with .boundary that is a LineString (.exterior exists, but is a LinearRing) + if vector is not None: + # We create a geodataframe with the geometry type + boundary_shp = gpd.GeoDataFrame(geometry=vector.ds.__getattr__(geometry_type), crs=vector.crs) + # We mask the pixels that make up the geometry type + mask_boundary = Vector(boundary_shp).create_mask(raster).squeeze() + + else: + # We mask target pixels + if target_values is not None: + mask_boundary = np.logical_or.reduce([raster.get_nanarray() == target_val for target_val in target_values]) + # Otherwise, all non-zero values are considered targets + else: + mask_boundary = raster.get_nanarray().astype(bool) + + # 2/ Now, we compute the distance matrix relative to the masked geometry type + if distance_unit.lower() == "georeferenced": + sampling: int | tuple[float | int, float | int] = raster.res + elif distance_unit.lower() == "pixel": + sampling = 1 + else: + raise ValueError('Distance unit must be either "georeferenced" or "pixel".') + + # If not all pixels are targets, then we compute the distance + non_targets = np.count_nonzero(mask_boundary) + if non_targets > 0: + proximity = distance_transform_edt(~mask_boundary, sampling=sampling) + # Otherwise, pass an array full of nodata + else: + proximity = np.ones(np.shape(mask_boundary)) * np.nan + + # 3/ If there was a vector input, apply the in_and_out argument to optionally mask inside/outside + if vector is not None: + if in_or_out == "both": + pass + elif in_or_out in ["in", "out"]: + mask_polygon = Vector(vector.ds).create_mask(raster).squeeze() + if in_or_out == "in": + proximity[~mask_polygon] = 0 + else: + proximity[mask_polygon] = 0 + else: + raise ValueError('The type of proximity must be one of "in", "out" or "both".') + + return proximity diff --git a/geoutils/geovector.py b/geoutils/geovector.py index 5e2b9830..4ac04e2b 100644 --- a/geoutils/geovector.py +++ b/geoutils/geovector.py @@ -364,7 +364,7 @@ def query(self: VectorType, expression: str, inplace: bool = False) -> VectorTyp """ Query the Vector dataset with a valid Pandas expression. - :param expression: A python-like expression to evaluate. Example: "col1 > col2" + :param expression: A python-like expression to evaluate. Example: "col1 > col2". :param inplace: Whether the query should modify the data in place or return a modified copy. :returns: Vector resulting from the provided query expression or itself if inplace=True. @@ -379,12 +379,105 @@ def query(self: VectorType, expression: str, inplace: bool = False) -> VectorTyp new_vector.__init__(self.ds.query(expression)) # type: ignore return new_vector # type: ignore + def proximity( + self, + raster: gu.Raster | None = None, + grid_size: tuple[int, int] = (1000, 1000), + geometry_type: str = "boundary", + in_or_out: Literal["in"] | Literal["out"] | Literal["both"] = "both", + distance_unit: Literal["pixel"] | Literal["georeferenced"] = "georeferenced", + ) -> gu.Raster: + """ + Proximity to this Vector's geometry computed for each cell of a Raster grid, or for a grid size with + this Vector's extent. + + If passed, the Raster georeferenced grid will be used to burn the proximity. If not, a grid size can be + passed to create a georeferenced grid with the bounds and CRS of this Vector. + + By default, the boundary of the Vector's geometry will be used. The full geometry can be used by + passing "geometry", + or any lower dimensional geometry attribute such as "centroid", "envelope" or "convex_hull". + See all geometry attributes in the Shapely documentation at https://shapely.readthedocs.io/. + + :param raster: Raster to burn the proximity grid on. + :param grid_size: If no Raster is provided, grid size to use with this Vector's extent and CRS + (defaults to 1000 x 1000). + :param geometry_type: Type of geometry to use for the proximity, defaults to 'boundary'. + :param in_or_out: Compute proximity only 'in' or 'out'-side the polygon, or 'both'. + :param distance_unit: Distance unit, either 'georeferenced' or 'pixel'. + + :return: Proximity raster. + """ + + # 0/ If no Raster is passed, create one on the Vector bounds of size 1000 x 1000 + if raster is None: + + # TODO: this bit of code is common in several vector functions (rasterize, etc): move out as common code? + # By default, use self's bounds + if self.bounds is None: + raise ValueError("To automatically rasterize on the vector, bounds need to be defined.") + + # Calculate raster shape + left, bottom, right, top = self.bounds + + # Calculate raster transform + transform = rio.transform.from_bounds(left, bottom, right, top, grid_size[0], grid_size[1]) + + raster = gu.Raster.from_array(data=np.zeros((1000, 1000)), transform=transform, crs=self.crs) + + proximity = gu.georaster.raster.proximity_from_vector_or_raster( + raster=raster, vector=self, geometry_type=geometry_type, in_or_out=in_or_out, distance_unit=distance_unit + ) + + return raster.copy(new_array=proximity) + + def buffer_metric(self, buffer_size: float) -> Vector: + """ + Buffer the vector in a metric system (UTM only). + + The outlines are projected to a local UTM, then reverted to the original projection after buffering. + + :param buffer_size: Buffering distance in meters + + :return: Buffered shapefile + """ + + from geoutils.projtools import latlon_to_utm, utm_to_epsg + + # Get a rough centroid in geographic coordinates (ignore the warning that it is not the most precise): + with warnings.catch_warnings(): + warnings.simplefilter(action="ignore", category=UserWarning) + shp_wgs84 = self.ds.to_crs(epsg=4326) + lat, lon = shp_wgs84.centroid.y.values[0], shp_wgs84.centroid.x.values[0] + del shp_wgs84 + + # Get the EPSG code of the local UTM + utm = latlon_to_utm(lat, lon) + epsg = utm_to_epsg(utm) + + # Reproject the shapefile in the local UTM + ds_utm = self.ds.to_crs(epsg=epsg) + + # Buffer the shapefile + ds_buffered = ds_utm.buffer(distance=buffer_size) + del ds_utm + + # Revert-project the shapefile in the original CRS + ds_buffered_origproj = ds_buffered.to_crs(crs=self.ds.crs) + del ds_buffered + + # Return a Vector object of the buffered GeoDataFrame + # TODO: Clarify what is conserved in the GeoSeries and what to pass the GeoDataFrame to not lose any attributes + vector_buffered = Vector(gpd.GeoDataFrame(geometry=ds_buffered_origproj.geometry, crs=self.ds.crs)) + + return vector_buffered + def get_bounds_projected(self, out_crs: CRS, densify_pts: int = 5000) -> rio.coords.BoundingBox: """ Return self's bounds in the given CRS. :param out_crs: Output CRS - :param densify_pts: Maximum points to be added between image corners to account for non linear edges. + :param densify_pts: Maximum points to be added between image corners to account for nonlinear edges. Reduce if time computation is really critical (ms) or increase if extent is not accurate enough. """ @@ -393,14 +486,14 @@ def get_bounds_projected(self, out_crs: CRS, densify_pts: int = 5000) -> rio.coo return new_bounds - def buffer_without_overlap(self, buffer_size: int | float, plot: bool = False) -> Vector: + def buffer_without_overlap(self, buffer_size: int | float, metric: bool = True, plot: bool = False) -> Vector: """ Returns a Vector object containing self's geometries extended by a buffer, without overlapping each other. The algorithm is based upon this tutorial: https://statnmap.com/2020-07-31-buffer-area-for-nearest-neighbour/. - The buffered polygons are created using Voronoi polygons in order to delineate the "area of influence" \ + The buffered polygons are created using Voronoi polygons in order to delineate the "area of influence" of each geometry. - The buffer is slightly inaccurate where two geometries touch, due to the nature of the Voronoi polygons,\ + The buffer is slightly inaccurate where two geometries touch, due to the nature of the Voronoi polygons, hence one geometry "steps" slightly on the neighbor buffer in some cases. The algorithm may also yield unexpected results on very simple geometries. @@ -416,12 +509,32 @@ def buffer_without_overlap(self, buffer_size: int | float, plot: bool = False) - >>> plt.show() # doctest: +SKIP :param buffer_size: Buffer size in self's coordinate system units. + :param metric: Whether to perform the buffering in a local metric system (default: True). :param plot: Set to True to show intermediate plots, useful for understanding or debugging. :returns: A Vector containing the buffered geometries. """ + + from geoutils.projtools import latlon_to_utm, utm_to_epsg + + # Project in local UTM if metric is True + if metric: + # Get a rough centroid in geographic coordinates (ignore the warning that it is not the most precise): + with warnings.catch_warnings(): + warnings.simplefilter(action="ignore", category=UserWarning) + shp_wgs84 = self.ds.to_crs(epsg=4326) + lat, lon = shp_wgs84.centroid.y.values[0], shp_wgs84.centroid.x.values[0] + del shp_wgs84 + + # Get the EPSG code of the local UTM + utm = latlon_to_utm(lat, lon) + epsg = utm_to_epsg(utm) + + gdf = self.ds.to_crs(epsg=epsg) + else: + gdf = self.ds + # Dissolve all geometries into one - gdf = self.ds merged = gdf.dissolve() # Add buffer around geometries @@ -478,6 +591,10 @@ def buffer_without_overlap(self, buffer_size: int | float, plot: bool = False) - ax4.set_title("Final buffer") plt.show() + # Reverse-project to the original CRS if metric is True + if metric: + merged_voronoi = merged_voronoi.to_crs(crs=self.crs) + return Vector(merged_voronoi) diff --git a/geoutils/projtools.py b/geoutils/projtools.py index ef36faa9..41184885 100644 --- a/geoutils/projtools.py +++ b/geoutils/projtools.py @@ -16,6 +16,53 @@ from shapely.geometry.polygon import Polygon +def latlon_to_utm(lat: float, lon: float) -> str: + """ + Get UTM zone for a given latitude and longitude coordinates. + + :param lat: Latitude coordinate. + :param lon: Longitude coordinate. + + :returns: UTM zone. + """ + + if not ( + isinstance(lat, (float, np.floating, int, np.integer)) + and isinstance(lon, (float, np.floating, int, np.integer)) + ): + raise TypeError("Latitude and longitude must be floats or integers.") + + if not -180 <= lon < 180: + raise ValueError("Longitude value is out of range [-180, 180[.") + if not -90 <= lat < 90: + raise ValueError("Latitude value is out of range [-90, 90[.") + + # Get UTM zone from name string of crs info + utm_zone = pyproj.database.query_utm_crs_info( + "WGS 84", area_of_interest=pyproj.aoi.AreaOfInterest(lon, lat, lon, lat) + )[0].name.split(" ")[-1] + + return str(utm_zone) + + +def utm_to_epsg(utm: str) -> int: + """ + Get EPSG code of UTM zone. + + :param utm: UTM zone. + + :return: EPSG of UTM zone. + """ + + # Whether UTM is passed as single or double digits, homogenize to single-digit + utm = str(int(utm[:-1])) + utm[-1].upper() + + # Get corresponding EPSG + epsg = pyproj.CRS(f"WGS 84 / UTM Zone {utm}").to_epsg() + + return int(epsg) + + def bounds2poly( boundsGeom: list[float] | rio.io.DatasetReader, in_crs: CRS | None = None, diff --git a/tests/test_georaster.py b/tests/test_georaster.py index 74c24698..113a7f52 100644 --- a/tests/test_georaster.py +++ b/tests/test_georaster.py @@ -26,6 +26,43 @@ DO_PLOT = False +def run_gdal_proximity( + input_raster: gu.Raster, target_values: list[float] | None, distunits: str = "GEO" +) -> np.ndarray: + """Run GDAL's ComputeProximity and return the read numpy array.""" + # Rasterio strongly recommends against importing gdal along rio, so this is done here instead. + from osgeo import gdal, gdalconst + + # Initiate empty GDAL raster for proximity output + drv = gdal.GetDriverByName("MEM") + proxy_ds = drv.Create("", input_raster.shape[1], input_raster.shape[0], 1, gdal.GetDataTypeByName("Float32")) + proxy_ds.GetRasterBand(1).SetNoDataValue(-9999) + + # Save input in temporary file to read with GDAL + # (avoids the nightmare of setting nodata, transform, crs in GDAL format...) + with tempfile.TemporaryDirectory() as temp_dir: + temp_path = os.path.join(temp_dir, "input.tif") + input_raster.save(temp_path) + ds_raster_in = gdal.Open(temp_path, gdalconst.GA_ReadOnly) + + # Define GDAL options + proximity_options = ["DISTUNITS=" + distunits] + if target_values is not None: + proximity_options.insert(0, "VALUES=" + ",".join([str(tgt) for tgt in target_values])) + + # Compute proximity + gdal.ComputeProximity(ds_raster_in.GetRasterBand(1), proxy_ds.GetRasterBand(1), proximity_options) + # Save array + proxy_array = proxy_ds.GetRasterBand(1).ReadAsArray().astype("float32") + proxy_array[proxy_array == -9999] = np.nan + + # Close GDAL datasets + proxy_ds = None + ds_raster_in = None + + return proxy_array + + class TestRaster: landsat_b4_path = examples.get_path("everest_landsat_b4") @@ -1886,6 +1923,115 @@ def test_polygonize(self, example: str) -> None: mask = img > value mask.polygonize(in_value=1) + # Test all options, with both an artificial Raster (that has all target values) and a real Raster + @pytest.mark.parametrize("distunits", ["GEO", "PIXEL"]) # type: ignore + # 0 and 1,2,3 are especially useful for the artificial Raster, and 112 for the real Raster + @pytest.mark.parametrize("target_values", [[1, 2, 3], [0], [112], None]) # type: ignore + @pytest.mark.parametrize( + "raster", + [ + gu.Raster(landsat_b4_path), + gu.Raster.from_array( + np.arange(25, dtype="int32").reshape(5, 5), transform=rio.transform.from_origin(0, 5, 1, 1), crs=4326 + ), + ], + ) # type: ignore + def test_proximity_against_gdal(self, distunits: str, target_values: list[float] | None, raster: gu.Raster) -> None: + """Test that proximity matches the results of GDAL for any parameter.""" + + # We generate proximity with GDAL and GeoUtils + gdal_proximity = run_gdal_proximity(raster, target_values=target_values, distunits=distunits) + # We translate distunits GDAL option into its GeoUtils equivalent + if distunits == "GEO": + distance_unit = "georeferenced" + else: + distance_unit = "pixel" + geoutils_proximity = ( + raster.proximity(distance_unit=distance_unit, target_values=target_values) + .data.data.squeeze() + .astype("float32") + ) + + # The results should be the same in all cases + try: + # In some cases, the proximity differs slightly (generally <1%) for complex settings + # (Landsat Raster with target of 112) + # It looks like GDAL might not have the right value, + # so this particular case is treated differently in tests + if target_values is not None and target_values[0] == 112 and raster.filename is not None: + + # Get index and number of not almost equal point (tolerance of 10-4) + ind_not_almost_equal = np.abs(gdal_proximity - geoutils_proximity) > 1e-04 + nb_not_almost_equal = np.count_nonzero(ind_not_almost_equal) + # Check that this is a minority of points (less than 0.5%) + assert nb_not_almost_equal < 0.005 * raster.width * raster.height + + # Replace these exceptions by zero in both + gdal_proximity[ind_not_almost_equal] = 0.0 + geoutils_proximity[ind_not_almost_equal] = 0.0 + # Check that all the rest is almost equal + assert np.allclose(gdal_proximity, geoutils_proximity, atol=1e-04, equal_nan=True) + + # Otherwise, results are exactly equal + else: + assert np.array_equal(gdal_proximity, geoutils_proximity, equal_nan=True) + + # For debugging + except Exception as exception: + + import matplotlib.pyplot as plt + + # Plotting the xdem and GDAL attributes for comparison (plotting "diff" can also help debug) + plt.subplot(121) + plt.imshow(gdal_proximity) + # plt.imshow(np.abs(gdal_proximity - geoutils_proximity)>0.1) + plt.colorbar() + plt.subplot(122) + plt.imshow(geoutils_proximity) + # plt.imshow(raster.data.data == 112) + plt.colorbar() + plt.show() + + # ind_not_equal = np.abs(gdal_proximity - geoutils_proximity)>0.1 + # print(gdal_proximity[ind_not_equal]) + # print(geoutils_proximity[ind_not_equal]) + + raise exception + + def test_proximity_parameters(self) -> None: + """ + Test that new (different to GDAL's) proximity parameters run. + No need to test the results specifically, as those rely entirely on the previous test with GDAL, + and tests in rasterize and shapely. + #TODO: Maybe add one test with an artificial vector to check it works as intended + """ + + # -- Test 1: with self's Raster alone -- + raster1 = gu.Raster(self.landsat_b4_path) + prox1 = raster1.proximity() + + # The raster should have the same extent, resolution and CRS + assert raster1.equal_georeferenced_grid(prox1) + + # It should change with target values specified + prox2 = raster1.proximity(target_values=[255]) + assert not np.array_equal(prox1.data, prox2.data) + + # -- Test 2: with a vector provided -- + vector = gu.Vector(self.everest_outlines_path) + + # With default options (boundary geometry) + raster1.proximity(vector=vector) + + # With the base geometry + raster1.proximity(vector=vector, geometry_type="geometry") + + # With another geometry option + raster1.proximity(vector=vector, geometry_type="centroid") + + # With only inside proximity + raster1.proximity(vector=vector, in_or_out="in") + def test_to_points(self) -> None: """Test the outputs of the to_points method and that it doesn't load if not needed.""" # Create a small raster to test point sampling on @@ -2048,6 +2194,58 @@ def test_equal(self) -> None: r2.set_nodata(34) assert r1 != r2 + def test_equal_georeferenced_grid(self) -> None: + """ + Test that equal for shape, crs and transform work as expected + """ + + # -- Test 1: based on a copy -- + r1 = self.r1 + r2 = r1.copy() + assert r1.equal_georeferenced_grid(r2) + + # Change data + r2.data += 1 + assert r1.equal_georeferenced_grid(r2) + + # Change mask (False by default) + r2 = r1.copy() + r2.data[0, 0] = np.ma.masked + assert r1.equal_georeferenced_grid(r2) + + # Change fill_value (999999 by default) + r2 = r1.copy() + r2.data.fill_value = 0 + assert r1.equal_georeferenced_grid(r2) + + # Change dtype + r2 = r1.copy() + r2 = r2.astype("float32") + assert r1.equal_georeferenced_grid(r2) + + # Change transform + r2 = r1.copy() + r2.transform = rio.transform.from_bounds(0, 0, 1, 1, self.width + 1, self.height) + assert not r1.equal_georeferenced_grid(r2) + + # Change CRS + r2 = r1.copy() + r2.crs = rio.crs.CRS.from_epsg(4326) + assert not r1.equal_georeferenced_grid(r2) + + # Change nodata + r2 = r1.copy() + r2.set_nodata(34) + assert r1.equal_georeferenced_grid(r2) + + # -- Test 2: based on another Raster with one different georeferenced grid attribute -- + + assert not r1.equal_georeferenced_grid(self.r1_wrong_crs) + + assert not r1.equal_georeferenced_grid(self.r1_wrong_shape) + + assert not r1.equal_georeferenced_grid(self.r1_wrong_transform) + # List of operations with two operands ops_2args = [ "__add__", diff --git a/tests/test_geovector.py b/tests/test_geovector.py index 6b28bf14..91187d35 100644 --- a/tests/test_geovector.py +++ b/tests/test_geovector.py @@ -120,6 +120,38 @@ def test_crop(self, data: list[str]) -> None: # Check that some features were indeed removed assert np.sum(~np.array(intersects_old)) > 0 + def test_proximity(self) -> None: + """ + The core functionality is already tested against GDAL in test_raster: just verify the vector-specific behaviour. + #TODO: add an artificial test as well (mirroring TODO in test_georaster) + """ + + vector = gu.Vector(self.everest_outlines_path) + + # -- Test 1: with a Raster provided -- + raster1 = gu.Raster(self.landsat_b4_crop_path) + prox1 = vector.proximity(raster=raster1) + + # The proximity should have the same extent, resolution and CRS + assert raster1.equal_georeferenced_grid(prox1) + + # With the base geometry + vector.proximity(raster=raster1, geometry_type="geometry") + + # With another geometry option + vector.proximity(raster=raster1, geometry_type="centroid") + + # With only inside proximity + vector.proximity(raster=raster1, in_or_out="in") + + # -- Test 2: with no Raster provided, just grid size -- + + # Default grid size + vector.proximity() + + # With specific grid size + vector.proximity(grid_size=(100, 100)) + class TestSynthetic: @@ -258,6 +290,69 @@ def test_generate_voronoi(self) -> None: with pytest.raises(ValueError, match=expected_message): voronoi = gu.geovector.generate_voronoi_polygons(self.vector.ds) + def test_buffer_metric(self) -> None: + """Check that metric buffering works""" + + # Case with two squares: test that the buffered area is without deformations + # https://epsg.io/32631 + utm31_x_center = 500000 + utm31_y_center = 4649776 + poly1_utm31 = Polygon( + [ + (utm31_x_center, utm31_y_center), + (utm31_x_center + 1, utm31_y_center), + (utm31_x_center + 1, utm31_y_center + 1), + (utm31_x_center, utm31_y_center + 1), + ] + ) + + poly2_utm31 = Polygon( + [ + (utm31_x_center + 10, utm31_y_center + 10), + (utm31_x_center + 11, utm31_y_center + 10), + (utm31_x_center + 11, utm31_y_center + 11), + (utm31_x_center + 10, utm31_y_center + 11), + ] + ) + + # We initiate the squares of size 1x1 in a UTM projection + two_squares = gu.Vector(gpd.GeoDataFrame(geometry=[poly1_utm31, poly2_utm31], crs="EPSG:32631")) + + # Their area should now be 1 for each polygon + assert two_squares.ds.area.values[0] == 1 + assert two_squares.ds.area.values[1] == 1 + + # We buffer them + two_squares_utm_buffered = two_squares.buffer_metric(buffer_size=1.0) + + # Their area should now be 1 (square) + 4 (buffer along the sides) + 4*(pi*1**2 /4) + # (buffer of corners = quarter-disks) + expected_area = 1 + 4 + np.pi + assert two_squares_utm_buffered.ds.area.values[0] == pytest.approx(expected_area, abs=0.01) + assert two_squares_utm_buffered.ds.area.values[1] == pytest.approx(expected_area, abs=0.01) + + # And the new GeoDataFrame should exactly match that of one buffer from the original one + direct_gpd_buffer = gu.Vector( + gpd.GeoDataFrame(geometry=two_squares.ds.buffer(distance=1.0).geometry, crs=two_squares.crs) + ) + assert_geodataframe_equal(direct_gpd_buffer.ds, two_squares_utm_buffered.ds) + + # Now, if we reproject the original vector in a non-metric system + two_squares_geographic = gu.Vector(two_squares.ds.to_crs(epsg=4326)) + # We buffer directly the Vector object in the non-metric system + two_squares_geographic_buffered = two_squares_geographic.buffer_metric(buffer_size=1.0) + # Then, we reproject that vector in the UTM zone + two_squares_geographic_buffered_reproj = gu.Vector( + two_squares_geographic_buffered.ds.to_crs(crs=two_squares.crs) + ) + + # Their area should now be the same as before for each polygon + assert two_squares_geographic_buffered_reproj.ds.area.values[0] == pytest.approx(expected_area, abs=0.01) + assert two_squares_geographic_buffered_reproj.ds.area.values[0] == pytest.approx(expected_area, abs=0.01) + + # And this time, it is the reprojected GeoDataFrame that should almost match (within a tolerance of 10e-06) + assert all(direct_gpd_buffer.ds.geom_almost_equals(two_squares_geographic_buffered_reproj.ds)) + def test_buffer_without_overlap(self) -> None: """ Check that non-overlapping buffer feature works. Does not work on simple geometries, so test on MultiPolygon. @@ -269,7 +364,7 @@ def test_buffer_without_overlap(self) -> None: # Check with buffers that should not overlap # ------------------------------------------ buffer_size = 2 - buffer = two_squares.buffer_without_overlap(buffer_size) + buffer = two_squares.buffer_without_overlap(buffer_size, metric=False) # Output should be of same size as input and same geometry type assert len(buffer.ds) == len(two_squares.ds) @@ -297,7 +392,7 @@ def test_buffer_without_overlap(self) -> None: # Case 2 - Check with buffers that overlap -> this case is actually not the expected result ! # ------------------------------- buffer_size = 5 - buffer = two_squares.buffer_without_overlap(buffer_size) + buffer = two_squares.buffer_without_overlap(buffer_size, metric=False) # Output should be of same size as input and same geometry type assert len(buffer.ds) == len(two_squares.ds) diff --git a/tests/test_projtools.py b/tests/test_projtools.py index f98c3f5f..7572a070 100644 --- a/tests/test_projtools.py +++ b/tests/test_projtools.py @@ -2,6 +2,7 @@ Test projtools """ import numpy as np +import pyproj.exceptions import pytest import geoutils as gu @@ -16,6 +17,75 @@ class TestProjTools: landsat_rgb_path = examples.get_path("everest_landsat_rgb") aster_dem_path = examples.get_path("exploradores_aster_dem") + def test_latlon_to_utm(self) -> None: + + # First: Check errors are raised when format is invalid + + # If format is invalid + with pytest.raises(TypeError): + pt.latlon_to_utm("1", 100) # type: ignore + # If values are outside limits: latitude above or below 90 + with pytest.raises(ValueError): + pt.latlon_to_utm(91, 0) + with pytest.raises(ValueError): + pt.latlon_to_utm(-91, 0) + # Or longitude above or below 180 + with pytest.raises(ValueError): + pt.latlon_to_utm(0, -181) + with pytest.raises(ValueError): + pt.latlon_to_utm(0, 181) + + # Second: check that the UTM zone is correct + # Lower left belongs to the zone, so 0, 0 should be in zone 31N + assert pt.latlon_to_utm(0, 0) == "30N" + # Test extreme zones + assert pt.latlon_to_utm(-79, -179) == "1S" + assert pt.latlon_to_utm(79, -179) == "1N" + assert pt.latlon_to_utm(-79, 179) == "60S" + assert pt.latlon_to_utm(79, 179) == "60N" + # Test some middles zones + assert pt.latlon_to_utm(1, -59) == "21N" + assert pt.latlon_to_utm(1, 61) == "41N" + assert pt.latlon_to_utm(-1, -121) == "10S" + assert pt.latlon_to_utm(-1, 119) == "50S" + + # Third, check that any floating or integer type works + assert pt.latlon_to_utm(0.1, 0.1) + assert pt.latlon_to_utm(np.float32(0.1), np.float32(0.1)) + assert pt.latlon_to_utm(np.int16(0), np.int16(0)) + + def test_utm_to_epsg(self) -> None: + """Check that the EPSG codes derived from UTM zones are correct""" + + # First: Check errors are raised when format is invalid + + # If there isn't 2 digits for the code + with pytest.raises(pyproj.exceptions.CRSError): + pt.utm_to_epsg("100N") + # If type is incorrect + with pytest.raises(TypeError): + pt.utm_to_epsg({"utm": "10N"}) # type: ignore + # If the code digits does not exist + with pytest.raises(pyproj.exceptions.CRSError): + pt.utm_to_epsg("61N") + # If the north-south zone letter is incorrect + with pytest.raises(pyproj.exceptions.CRSError): + pt.utm_to_epsg("61E") + + # Second: Check that the EPSG code is correct + # https://epsg.io/32601 + assert pt.utm_to_epsg("1N") == 32601 + # https://epsg.io/32701 + assert pt.utm_to_epsg("1S") == 32701 + # https://epsg.io/32&660 + assert pt.utm_to_epsg("60N") == 32660 + # https://epsg.io/32760 + assert pt.utm_to_epsg("60S") == 32760 + + # Third: Check that different format work: single digit, lower-case + assert pt.utm_to_epsg("1N") == pt.utm_to_epsg("01N") == pt.utm_to_epsg("01n") + assert pt.utm_to_epsg("08s") == pt.utm_to_epsg("8S") == pt.utm_to_epsg("08S") + @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path]) # type: ignore def test_latlon_reproject(self, example: str) -> None: """