Fixed issue #13 (TypeError when a MultiPolygon instance is passed to...

Fixed issue #13 (TypeError when a MultiPolygon instance is passed to fill_holes_within_poly()). Added tests for fill_holes_within_poly() and get_overlap_polygon(). Tests are now also executable on Windows. Updated version info and HISTORY.rst.
Signed-off-by: Daniel Scheffler <danschef@gfz-potsdam.de>

HISTORY.rst

@@ -2,6 +2,14 @@
 History
 =======
 
+0.15.9 (2020-10-30)
+-------------------
+
+* Fixed issue #13 (TypeError when a MultiPolygon instance is passed to fill_holes_within_poly()).
+* Added tests for fill_holes_within_poly() and get_overlap_polygon()
+* Tests are now also executable on Windows.
+
+
 0.15.8 (2020-10-26)
 -------------------
 

py_tools_ds/geo/vector/topology.py

@@ -59,7 +59,8 @@ def get_overlap_polygon(poly1, poly2, v=False):
 
         overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
         if v:
-            print('%.2f percent of the image to be shifted is covered by the reference image.' % overlap_percentage)
+            print('%.2f percent of the image to be shifted is covered by the reference image.'
+                  % overlap_percentage)  # pragma: no cover
         return {'overlap poly': overlap_poly, 'overlap percentage': overlap_percentage,
                 'overlap area': overlap_poly.area}
     else:
@@ -158,6 +159,7 @@ def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0):
     elif inner_poly.intersects(outer_poly.buffer(tolerance)):
         # check if all vertices intersect with outer poly
         GDF = GeoDataFrame(np.swapaxes(np.array(inner_poly.exterior.coords.xy), 0, 1), columns=['X', 'Y'])
+        # noinspection PyTypeChecker
         return False in GDF.apply(lambda GDF_row: Point(GDF_row.X, GDF_row.Y).intersects(outer_poly), axis=1).values
     else:
         # inner_poly does not intersect out_poly -> all vertices are outside
@@ -166,7 +168,7 @@ def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0):
 
 def fill_holes_within_poly(poly):
     # type: (Union[Polygon, MultiPolygon]) -> Polygon
-    """Fill the holes within a shapely Polygon or MultiPolygon and returns a Polygon with only the outer boundary.
+    """Fill the holes within a shapely Polygon or MultiPolygon and return a Polygon with only the outer boundary.
 
     :param poly:  <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection>
     :return:
@@ -180,14 +182,27 @@ def fill_holes_within_poly(poly):
 
     else:  # 'MultiPolygon'
         gdf = GeoDataFrame(columns=['geometry'])
-        gdf['geometry'] = poly
+        gdf['geometry'] = poly.geoms
 
         # get the area of each polygon of the multipolygon EXCLUDING the gaps in it
+        # noinspection PyTypeChecker
         gdf['area_filled'] = gdf.apply(
             lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1)
 
         largest_poly_filled = gdf.loc[gdf['area_filled'].idxmax()]['geometry']
 
+        # noinspection PyTypeChecker
+        gdf['within_or_equal'] = gdf.apply(
+            lambda GDF_row:
+            GDF_row.geometry.within(largest_poly_filled) or
+            GDF_row.geometry.equals(largest_poly_filled),
+            axis=1)
+
+        if False in gdf.within_or_equal:
+            warnings.warn(RuntimeWarning('The given MultiPolygon contains disjunct polygones outside of the largest '
+                                         'polygone. fill_holes_within_poly() will only return the largest polygone as '
+                                         'a filled version.'))
+
         # return the outer boundary of the largest polygon
         filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))
 

py_tools_ds/version.py

@@ -19,5 +19,5 @@
 # You should have received a copy of the GNU Lesser General Public License along
 # with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-__version__ = '0.15.8'
-__versionalias__ = '20201026_01'
+__version__ = '0.15.9'
+__versionalias__ = '20201030_01'

setup.py

@@ -47,8 +47,8 @@ requirements = [
     'six',
     'spectral'
 ]
-setup_requirements = []  # TODO(danschef): put setup requirements (distutils extensions, etc.) here
-test_requirements = requirements + ["coverage", "nose", "nose2", "nose-htmloutput", "rednose"]
+setup_requirements = ['setuptools']
+test_requirements = requirements + ["coverage", "nose", "nose2", "nose-htmloutput", "rednose", "shapely"]
 
 setup(
     name='py_tools_ds',

tests/test_geo/test_coord_grid.py

@@ -47,3 +47,7 @@ class Test_move_shapelyPoly_to_image_grid(unittest.TestCase):
         poly_on_grid = move_shapelyPoly_to_image_grid(poly_local, (0, 1, 0, 0, 0, -1), rows=6281, cols=11162)
         self.assertTrue(isinstance(poly_on_grid, Polygon))
         self.assertEqual(str(poly_on_grid), 'POLYGON ((5708 -3262, 5708 -3006, 5452 -3006, 5452 -3262, 5708 -3262))')
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_geo/test_coord_trafo.py

@@ -30,7 +30,7 @@ test_coord_trafo
 Tests for `py_tools_ds.geo.coord_trafo` module.
 """
 
-from unittest import TestCase
+from unittest import TestCase, main
 from shapely.geometry import Polygon
 import numpy as np
 
@@ -58,3 +58,7 @@ class Test_transform_any_prj(TestCase):
         self.assertTrue(isinstance(out_y, float))
         self.assertTrue(np.all(np.isclose(np.array([out_x, out_y]),
                                           np.array([utm_x, utm_y]))))
+
+
+if __name__ == '__main__':
+    main()

tests/test_geo/test_map_info.py

@@ -155,3 +155,7 @@ class Test_mapinfo2geotransform(unittest.TestCase):
         gt = mapinfo2geotransform(map_info_local_rotated)
         self.assertTrue(isinstance(gt, (tuple, list)))
         self.assertEqual(gt, geotransform_local_rotated)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_geo/test_projection.py

@@ -81,3 +81,7 @@ class Test_EPSG2WKT(unittest.TestCase):
         wkt = EPSG2WKT(self.epsg_utm)
         self.assertTrue(isinstance(wkt, str), "EPSG2WKT returned a %s object instead of a string!" % type(wkt))
         self.assertNotEqual(wkt, "", msg="EPSG2WKT returned an empty WKT string!")
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_geo/test_vector/test_geometry.py

@@ -29,7 +29,7 @@ test_geometry
 Tests for `py_tools_ds.geo.vector.geometry` module.
 """
 
-from unittest import TestCase
+from unittest import TestCase, main
 
 from py_tools_ds.geo.vector.geometry import boxObj, get_winPoly
 
@@ -219,3 +219,7 @@ class Test_get_winPoly(TestCase):
 
     def test_get_winPoly_match_grid(self):
         get_winPoly(wp_imYX=self.wp_imYX, ws=self.ws, gt=self.gt, match_grid=True)
+
+
+if __name__ == '__main__':
+    main()

tests/test_geo/test_vector/test_topology.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# py_tools_ds
+#
+# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
+#
+# This software was developed within the context of the GeoMultiSens project funded
+# by the German Federal Ministry of Education and Research
+# (project grant code: 01 IS 14 010 A-C).
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License along
+# with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+test_topology
+-------------
+
+Tests for `py_tools_ds.geo.vector.topology` module.
+"""
+
+from unittest import TestCase, main
+from shapely.geometry import Polygon, MultiPolygon, Point
+
+from py_tools_ds.geo.vector.topology import get_overlap_polygon, fill_holes_within_poly
+
+
+class Test_get_overlap_polygon(TestCase):
+    def setUp(self):
+        self.polygon_1 = Polygon([(0, 0), (0, 5), (5, 5), (5, 0)])
+        self.polygon_2 = Polygon([(1, 0), (1, 5), (6, 5), (6, 0)])
+        self.polygon_3 = Polygon([(-1, 0), (-1, -5), (-6, -5), (-6, 0)])
+
+    def test_with_overlap(self):
+        out = get_overlap_polygon(self.polygon_1, self.polygon_2)
+        self.assertIsInstance(out, dict)
+
+        self.assertTrue(out['overlap poly'].equals(Polygon([(1, 0), (1, 5), (5, 5), (5, 0)])))
+        self.assertEqual(out['overlap percentage'], 80)
+        self.assertEqual(out['overlap area'], 20)
+
+    def test_no_overlap(self):
+        out = get_overlap_polygon(self.polygon_1, self.polygon_3)
+        self.assertIsInstance(out, dict)
+
+        self.assertTrue(out['overlap poly'] is None)
+        self.assertEqual(out['overlap percentage'], 0)
+        self.assertEqual(out['overlap area'], 0)
+
+    # TODO test the case where the output polygon is a multipolygon
+
+
+class Test_fill_holes_within_poly(TestCase):
+    def setUp(self):
+        self.polygon_1 = Polygon([(0, 0), (5, 0), (5, 5), (0, 5)])  # large poly
+        self.polygon_2 = Polygon([(1, 1), (1, 2), (2, 2), (2, 1)])  # small poly within large poly
+        self.polygon_3 = Polygon([(3, 3), (3, 4), (5, 4), (4, 3)])  # small poly within large poly
+        self.multipolygon_1 = MultiPolygon([self.polygon_1, self.polygon_2])
+        self.multipolygon_2 = MultiPolygon([self.polygon_2, self.polygon_3])
+
+    def test_polygon(self):
+        filled_poly = fill_holes_within_poly(self.polygon_1)
+        self.assertIsInstance(filled_poly, Polygon)
+        self.assertTrue(filled_poly == self.polygon_1)
+
+    def test_multipolygon_with_smaller_poly_inside(self):
+        filled_poly = fill_holes_within_poly(self.multipolygon_1)
+        self.assertIsInstance(filled_poly, Polygon)
+        self.assertTrue(filled_poly == self.polygon_1)
+
+    def test_multipolygon_is_disjunct(self):
+        with self.assertWarns(RuntimeWarning):
+            filled_poly = fill_holes_within_poly(self.multipolygon_2)
+
+        self.assertIsInstance(filled_poly, Polygon)
+        self.assertTrue(filled_poly == self.polygon_3)
+
+    def test_empty_multipolygon(self):
+        with self.assertRaises(ValueError):
+            fill_holes_within_poly(MultiPolygon())
+
+    def test_wrong_input(self):
+        with self.assertRaises(ValueError):
+            # noinspection PyTypeChecker
+            fill_holes_within_poly(Point((0, 1)))
+
+
+if __name__ == '__main__':
+    main()

tests/test_io/test_raster/test_gdal.py

@@ -30,7 +30,7 @@ test_gdal
 Tests for `py_tools_ds.io.raster.gdal` module.
 """
 
-from unittest import TestCase
+from unittest import TestCase, main
 import numpy as np
 from gdal import Dataset
 from pandas import DataFrame
@@ -95,3 +95,7 @@ class Test_get_GDAL_driverList(TestCase):
     def test_check_output(self):
         drvList = get_GDAL_driverList()
         self.assertIsInstance(drvList, DataFrame)
+
+
+if __name__ == '__main__':
+    main()

tests/test_numeric/test_array.py

@@ -58,3 +58,7 @@ class Test_get_array_tilebounds(unittest.TestCase):
 
     def test_pixel(self):
         self._run_and_validate(tile_shape=(1, 1))
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_numeric/test_vector.py

@@ -29,7 +29,7 @@ test_vector
 Tests for `py_tools_ds.numeric.vector` module.
 """
 
-from unittest import TestCase
+from unittest import TestCase, main
 
 from py_tools_ds.numeric.vector import find_nearest
 
@@ -89,3 +89,7 @@ class Test_find_nearest(TestCase):
         # value above
         out = find_nearest(self.array, 330060.7, roundAlg='on', tolerance=0.0001)
         self.assertEqual(out, 330060.69999999984)  # would return 330070.7999999998 without tolerance
+
+
+if __name__ == '__main__':
+    main()