Commit 90586987 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Added output verification for get_overlap_polygon()

geo.vector.topology:
- get_overlap_polygon(): added output verification (automatic hole filling in case of MultiPolygon)

- updated __version__
parent c8af7bd3
......@@ -15,7 +15,7 @@ __all__=[#'compatibility',
'similarity',
'GeoArray']
__version__ = '20170221_01'
__version__ = '20170221_02'
__author__='Daniel Scheffler'
# Validate GDAL version
......
......@@ -24,6 +24,12 @@ def get_overlap_polygon(poly1, poly2, v=False):
"""
overlap_poly = poly1.intersection(poly2)
if not overlap_poly.is_empty:
# check if output is MULTIPOLYGON -> if yes, convert to POLYGON
if overlap_poly.geom_type == 'MultiPolygon':
overlap_poly = fill_holes_within_poly(overlap_poly)
assert overlap_poly.geom_type=='Polygon', \
"get_overlap_polygon() did not return geometry type 'Polygon' but %s." %overlap_poly.geom_type
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)
return {'overlap poly':overlap_poly, 'overlap percentage':overlap_percentage, 'overlap area':overlap_poly.area}
......
......@@ -420,10 +420,9 @@ class GeoArray(object):
def footprint_poly(self):
# FIXME should return polygon in image coordinates if no projection is available
"""Get the footprint polygon of the associated image array (returns an instance of shapely.geometry.Polygon."""
if self._footprint_poly is not None:
return self._footprint_poly
else:
"""Get the footprint polygon of the associated image array (returns an instance of shapely.geometry.Polygon.
"""
if self._footprint_poly is None:
assert self.mask_nodata is not None, 'A nodata mask is needed for calculating the footprint polygon. '
if np.std(self.mask_nodata[:])==0:
# do not run raster2polygon if whole image is filled with data
......@@ -450,8 +449,7 @@ class GeoArray(object):
# assert self.GeoArray.box.mapPoly.contains(Point(XY)) or self.GeoArray.box.mapPoly.touches(Point(XY)), \
# "The corner position '%s' is outside of the %s." % (XY, self.imName)
return self._footprint_poly
return self._footprint_poly
@footprint_poly.setter
......@@ -797,7 +795,8 @@ class GeoArray(object):
driver = gdal.GetDriverByName(fmt)
if driver is None:
raise Exception("'%s' is not a supported GDAL driver." % fmt)
raise Exception("'%s' is not a supported GDAL driver. Refer to www.gdal.org/formats_list.html for full "
"list of GDAL driver codes." % fmt)
if not os.path.isdir(os.path.dirname(out_path)):
os.makedirs(os.path.dirname(out_path))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment