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

Merge branch 'bugfix/fix_get_overlap_polygon'

parents 063b540d d9c1dc93
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler" __author__ = "Daniel Scheffler"
import math import math
import warnings import warnings
import numpy as np import numpy as np
from typing import TypeVar, Union
from geopandas import GeoDataFrame from geopandas import GeoDataFrame
from shapely.geometry import shape, Polygon, box, Point from shapely.geometry import shape, Polygon, box, Point, MultiPolygon
from ..coord_trafo import mapYX2imYX from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord from ..coord_grid import find_nearest_grid_coord
def get_overlap_polygon(poly1, poly2, v=False): def get_overlap_polygon(poly1, poly2, v=False):
...@@ -22,19 +22,25 @@ 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: overlap percentage as float value [%]
:return: area of overlap polygon :return: area of overlap polygon
""" """
# compute overlap polygon
overlap_poly = poly1.intersection(poly2) 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: 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': if overlap_poly.geom_type == 'MultiPolygon':
overlap_poly = fill_holes_within_poly(overlap_poly) overlap_poly = fill_holes_within_poly(overlap_poly)
assert overlap_poly.geom_type=='Polygon', \ assert overlap_poly.geom_type == 'Polygon', \
"get_overlap_polygon() did not return geometry type 'Polygon' but %s." %overlap_poly.geom_type "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 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) if v:
return {'overlap poly':overlap_poly, 'overlap percentage':overlap_percentage, 'overlap area':overlap_poly.area} 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: 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): def get_footprint_polygon(CornerLonLat, fix_invalid=False):
...@@ -47,7 +53,7 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False): ...@@ -47,7 +53,7 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False):
if fix_invalid: if fix_invalid:
with warnings.catch_warnings(): with warnings.catch_warnings():
warnings.simplefilter('ignore') # FIXME not working warnings.simplefilter('ignore') # FIXME not working
outpoly = Polygon(CornerLonLat) outpoly = Polygon(CornerLonLat)
if not outpoly.is_valid: if not outpoly.is_valid:
...@@ -56,30 +62,31 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False): ...@@ -56,30 +62,31 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False):
outpoly = Polygon(CornerLonLat) outpoly = Polygon(CornerLonLat)
assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \ 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 return outpoly
def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im): 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 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 (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 xmin, ymin, xmax, ymax = int(xmin), math.ceil(ymin), math.ceil(xmax), int(
return (ymax,xmin),(ymax,xmax),(ymin,xmax),(ymin,xmin) # UL_YX,UR_YX,LR_YX,LL_YX 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): def get_largest_onGridPoly_within_poly(outerPoly, gt, rows, cols):
oP_xmin,oP_ymin,oP_xmax,oP_ymax = outerPoly.bounds 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') 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') xmax, ymin = find_nearest_grid_coord((oP_xmax, oP_ymin), gt, rows, cols, direction='NW')
return box(xmin, ymin, xmax, ymax) return box(xmin, ymin, xmax, ymax)
def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly): def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly):
"""Returns the smallest box that matches the coordinate grid of the given geotransform. """Returns the smallest box that matches the coordinate grid of the given geotransform.
The returned shapely polygon contains image coordinates.""" 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): def find_line_intersection_point(line1, line2):
...@@ -116,16 +123,20 @@ def polyVertices_outside_poly(inner_poly, outer_poly): ...@@ -116,16 +123,20 @@ def polyVertices_outside_poly(inner_poly, outer_poly):
def fill_holes_within_poly(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. """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: :return:
""" """
if poly.geom_type not in ['Polygon', 'MultiPolygon']:
raise ValueError("Unexpected geometry type %s." % poly.geom_type)
if poly.geom_type == 'Polygon': if poly.geom_type == 'Polygon':
# return only the exterior 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 = GeoDataFrame(columns=['geometry'])
gdf['geometry'] = poly gdf['geometry'] = poly
...@@ -133,7 +144,9 @@ def fill_holes_within_poly(poly): ...@@ -133,7 +144,9 @@ def fill_holes_within_poly(poly):
gdf['area_filled'] = gdf.apply( gdf['area_filled'] = gdf.apply(
lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1) 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 the outer boundary of the largest polygon
return Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1)) filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))
\ No newline at end of file
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