diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 1772dd8..4539976 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -2,6 +2,11 @@ Changelog ========= +Version 0.2.3 +============= + +- fixed rounding of world system coordinates + Version 0.2.2 ============= diff --git a/src/geospade/raster.py b/src/geospade/raster.py index 97dca11..0dddb0c 100644 --- a/src/geospade/raster.py +++ b/src/geospade/raster.py @@ -23,6 +23,7 @@ from geospade.tools import bbox_to_polygon from geospade.tools import rasterise_polygon from geospade.tools import rel_extent +from geospade.tools import _round_polygon_coords from geospade.transform import build_geotransform from geospade.transform import xy2ij from geospade.transform import ij2xy @@ -689,7 +690,7 @@ def touches(self, other, sref=None) -> bool: True if both geometries touch each other, false if not. """ - return self.boundary.Touches(other) + return _round_polygon_coords(self.boundary, DECIMALS).Touches(_round_polygon_coords(other, DECIMALS)) @_align_geom(align=True) def within(self, other, sref=None) -> bool: @@ -771,7 +772,7 @@ def slice_by_geom(self, other, sref=None, snap_to_grid=True, inplace=False, **kw return intersection = self.boundary.Intersection(other) - bbox = np.around(intersection.GetEnvelope(), decimals=DECIMALS) + bbox = intersection.GetEnvelope() if snap_to_grid: new_ll_x, new_ll_y = self.snap_to_grid(bbox[0], bbox[2], px_origin="ll") new_ur_x, new_ur_y = self.snap_to_grid(bbox[1], bbox[3], px_origin="ur") @@ -1286,8 +1287,9 @@ def __eq__(self, other) -> bool: True if both raster geometries are the same, otherwise false. """ - - return self.outer_boundary_corners == other.outer_boundary_corners and \ + this_corners = np.around(np.array(self.outer_boundary_corners), decimals=DECIMALS) + other_corners = np.around(np.array(other.outer_boundary_corners), decimals=DECIMALS) + return np.all(this_corners == other_corners) and \ self.n_rows == other.n_rows and \ self.n_cols == other.n_cols diff --git a/src/geospade/tools.py b/src/geospade/tools.py index 1a98a3c..633c56d 100644 --- a/src/geospade/tools.py +++ b/src/geospade/tools.py @@ -80,8 +80,8 @@ def polar_point(x_ori, y_ori, dist, angle, deg=True) -> Tuple[np.ndarray or floa """ angle = np.radians(angle) if deg else angle - x_pp = np.around(x_ori + dist * np.cos(angle), decimals=DECIMALS) - y_pp = np.around(y_ori + dist * np.sin(angle), decimals=DECIMALS) + x_pp = x_ori + dist * np.cos(angle) + y_pp = y_ori + dist * np.sin(angle) return x_pp, y_pp diff --git a/src/geospade/transform.py b/src/geospade/transform.py index 13d57ea..d3209ed 100644 --- a/src/geospade/transform.py +++ b/src/geospade/transform.py @@ -166,8 +166,6 @@ def transform_geom(geom, other_sref) -> ogr.Geometry: geometry_out = geom.Clone() # transform geometry to new spatial reference system. geometry_out.TransformTo(other_sref.osr_sref) - # round geometry coordinates - geometry_out = _round_geom_coords(geometry_out, DECIMALS) # assign new spatial reference system geometry_out.AssignSpatialReference(other_sref.osr_sref) diff --git a/tests/test_raster.py b/tests/test_raster.py index 4da2c3a..7605888 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -18,6 +18,13 @@ from geospade.raster import RasterGeometry from geospade.raster import MosaicGeometry from geospade.raster import RegularMosaicGeometry +from geospade import DECIMALS + + +def assert_extent(this_extent, other_extent): + this_extent = np.around(np.array(this_extent), decimals=DECIMALS) + other_extent = np.around(np.array(other_extent), decimals=DECIMALS) + assert np.all(this_extent == other_extent) class RasterGeometryTest(unittest.TestCase): @@ -58,14 +65,14 @@ def setUp(self): def test_from_extent(self): """ Tests setting up a raster geometry from a given extent. """ - self.assertTupleEqual(self.raster_geom.outer_boundary_extent, self.extent) + assert_extent(self.raster_geom.outer_boundary_extent, self.extent) def test_from_geom(self): """ Tests setting up a raster geometry from a given geometry. """ raster_geom = RasterGeometry.from_geometry(self.ogr_geom, self.x_pixel_size, self.y_pixel_size) - self.assertTupleEqual(raster_geom.outer_boundary_corners, tuple(self.sh_geom.exterior.coords)[:-1]) + assert_extent(raster_geom.outer_boundary_corners, tuple(self.sh_geom.exterior.coords)[:-1]) def test_from_json(self): """ Tests the creation of a raster geometry from a JSON file containing a raster geometry definition. """ @@ -99,7 +106,7 @@ def test_get_common_geom(self): ur_x, ur_y = raster_geom_c.rc2xy(0, raster_geom_b.n_cols - 1, px_origin="ur") new_extent = (ll_x, ll_y, ur_x, ur_y) - self.assertTupleEqual(raster_geom.outer_boundary_extent, new_extent) + assert_extent(raster_geom.outer_boundary_extent, new_extent) # test if error is raised, when raster geometries with different resolutions are joined try: @@ -156,7 +163,7 @@ def test_vertices(self): (self.extent[2], self.extent[3]), (self.extent[2], self.extent[1])) - self.assertTupleEqual(self.raster_geom.outer_boundary_corners, vertices) + assert_extent(self.raster_geom.outer_boundary_corners, vertices) def test_x_coords(self): """ Tests coordinate retrieval along x dimension. """ @@ -241,7 +248,7 @@ def test_intersection(self): self.x_pixel_size, self.y_pixel_size) raster_geom_intsct = self.raster_geom.slice_by_geom(raster_geom_shifted, inplace=False) - assert raster_geom_intsct.outer_boundary_extent == extent_intsct + assert_extent(raster_geom_intsct.outer_boundary_extent, extent_intsct) # test intersection with no overlap extent_no_ovlp = (self.extent[2] + 1., self.extent[3] + 1., @@ -531,7 +538,7 @@ def test_slice_tiles_by_geom(self): geom = any_geom2ogr_geom(self._get_roi(), sref=self.mosaic_geom.sref) intscted_tiles = self.mosaic_geom.slice_by_geom(geom, inplace=False).tiles outer_boundary_extent = RasterGeometry.from_raster_geometries(intscted_tiles).outer_boundary_extent - self.assertTupleEqual(outer_boundary_extent, self._get_roi()) + assert_extent(outer_boundary_extent, self._get_roi()) def test_slice_by_geom(self): """ Tests sub-setting a mosaic geometry with another geometry. """ @@ -668,7 +675,7 @@ def test_slice_tiles_by_geom(self): geom = any_geom2ogr_geom(self._get_roi(), sref=self.mosaic_geom.sref) intscted_tiles = self.mosaic_geom.slice_by_geom(geom, inplace=False).tiles outer_boundary_extent = RasterGeometry.from_raster_geometries(intscted_tiles).outer_boundary_extent - self.assertTupleEqual(outer_boundary_extent, self._get_roi()) + assert_extent(outer_boundary_extent, self._get_roi()) def test_slice_by_geom(self): """ Tests sub-setting a regular mosaic geometry with another geometry. """