Merge pull request #21 from TUW-GEO/round_bug_fix

Rounding bug fix

CHANGELOG.rst

@@ -2,6 +2,11 @@
 Changelog
 =========
 
+Version 0.2.3
+=============
+
+- fixed rounding of world system coordinates
+
 Version 0.2.2
 =============
 

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
 

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
 

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)
 

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. """