Commit ad3f4d04 authored by Daniel Scheffler's avatar Daniel Scheffler Committed by Daniel Scheffler
Browse files

Fix issue #5 (get_overlap_polygon() did not return geometry type 'Polygon' but...

Fix issue #5 (get_overlap_polygon() did not return geometry type 'Polygon' but GeometryCollection.").
parent ebed4b23
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
import math
import warnings
import numpy as np
from typing import TypeVar, Union
from geopandas import GeoDataFrame
from shapely.geometry import shape, Polygon, box, Point
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
from shapely.geometry import shape, Polygon, box, Point, MultiPolygon
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
def get_overlap_polygon(poly1, poly2, v=False):
......@@ -22,19 +22,25 @@ def get_overlap_polygon(poly1, poly2, v=False):
:return: overlap percentage as float value [%]
:return: area of overlap polygon
"""
# compute overlap polygon
overlap_poly = poly1.intersection(poly2)
if overlap_poly.geom_type == 'GeometryCollection':
overlap_poly = overlap_poly.buffer(0) # converts 'GeometryCollection' to 'MultiPolygon'
if not overlap_poly.is_empty:
# check if output is MULTIPOLYGON -> if yes, convert to POLYGON
# check if output is MultiPolygon or GeometryCollection -> 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
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}
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}
else:
return {'overlap poly':None, 'overlap percentage':0, 'overlap area':0}
return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0}
def get_footprint_polygon(CornerLonLat, fix_invalid=False):
......@@ -47,7 +53,7 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False):
if fix_invalid:
with warnings.catch_warnings():
warnings.simplefilter('ignore') # FIXME not working
warnings.simplefilter('ignore') # FIXME not working
outpoly = Polygon(CornerLonLat)
if not outpoly.is_valid:
......@@ -56,30 +62,31 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False):
outpoly = Polygon(CornerLonLat)
assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \
'Got coordinates %s.' %CornerLonLat
'Got coordinates %s.' % CornerLonLat
return outpoly
def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im):
xmin,ymin,xmax,ymax = Polygon([(i[1],i[0]) for i in box_mapYX]).bounds # map-bounds box_mapYX
(ymin,xmin),(ymax,xmax) = mapYX2imYX([ymin,xmin],gt_im), mapYX2imYX([ymax,xmax],gt_im) # image coord bounds
xmin,ymin,xmax,ymax = int(xmin),math.ceil(ymin), math.ceil(xmax), int(ymax) # round min coords off and max coords on
return (ymax,xmin),(ymax,xmax),(ymin,xmax),(ymin,xmin) # UL_YX,UR_YX,LR_YX,LL_YX
xmin, ymin, xmax, ymax = Polygon([(i[1], i[0]) for i in box_mapYX]).bounds # map-bounds box_mapYX
(ymin, xmin), (ymax, xmax) = mapYX2imYX([ymin, xmin], gt_im), mapYX2imYX([ymax, xmax], gt_im) # image coord bounds
xmin, ymin, xmax, ymax = int(xmin), math.ceil(ymin), math.ceil(xmax), int(
ymax) # round min coords off and max coords on
return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin) # UL_YX,UR_YX,LR_YX,LL_YX
def get_largest_onGridPoly_within_poly(outerPoly,gt,rows,cols):
oP_xmin,oP_ymin,oP_xmax,oP_ymax = outerPoly.bounds
xmin,ymax = find_nearest_grid_coord((oP_xmin,oP_ymax),gt,rows,cols,direction='SE')
xmax,ymin = find_nearest_grid_coord((oP_xmax,oP_ymin),gt,rows,cols,direction='NW')
def get_largest_onGridPoly_within_poly(outerPoly, gt, rows, cols):
oP_xmin, oP_ymin, oP_xmax, oP_ymax = outerPoly.bounds
xmin, ymax = find_nearest_grid_coord((oP_xmin, oP_ymax), gt, rows, cols, direction='SE')
xmax, ymin = find_nearest_grid_coord((oP_xmax, oP_ymin), gt, rows, cols, direction='NW')
return box(xmin, ymin, xmax, ymax)
def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly):
"""Returns the smallest box that matches the coordinate grid of the given geotransform.
The returned shapely polygon contains image coordinates."""
xmin,ymin,xmax,ymax = shapelyPoly.bounds # image_coords-bounds
xmin, ymin, xmax, ymax = shapelyPoly.bounds # image_coords-bounds
return box(int(xmin),int(ymin), math.ceil(xmax), math.ceil(ymax)) # round min coords off and max coords on
return box(int(xmin), int(ymin), math.ceil(xmax), math.ceil(ymax)) # round min coords off and max coords on
def find_line_intersection_point(line1, line2):
......@@ -116,16 +123,20 @@ def polyVertices_outside_poly(inner_poly, outer_poly):
def fill_holes_within_poly(poly):
# type: (Union[TypeVar(Polygon), TypeVar(MultiPolygon)]) -> TypeVar(Polygon)
"""Fills the holes within a shapely Polygon or MultiPolygon and returns a Polygon with only the outer boundary.
:param poly: <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>
:param poly: <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection>
:return:
"""
if poly.geom_type not in ['Polygon', 'MultiPolygon']:
raise ValueError("Unexpected geometry type %s." % poly.geom_type)
if poly.geom_type == 'Polygon':
# return only the exterior polygon
return Polygon(poly.exterior)
filled_poly = Polygon(poly.exterior)
elif poly.geom_type == 'MultiPolygon':
else: # 'MultiPolygon'
gdf = GeoDataFrame(columns=['geometry'])
gdf['geometry'] = poly
......@@ -133,7 +144,9 @@ def fill_holes_within_poly(poly):
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.ix[gdf['area_filled'].idxmax()]['geometry']
largest_poly_filled = gdf.loc[gdf['area_filled'].idxmax()]['geometry']
# return the outer boundary of the largest polygon
return Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))
\ No newline at end of file
filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))
return filled_poly
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