topology.py 7.31 KB
Newer Older
1
2
3
# -*- coding: utf-8 -*-

import math
4
import warnings
Daniel Scheffler's avatar
Daniel Scheffler committed
5
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
6
from typing import Union  # noqa F401  # flake8 issue
7

Daniel Scheffler's avatar
Daniel Scheffler committed
8
from geopandas import GeoDataFrame
Daniel Scheffler's avatar
Daniel Scheffler committed
9
10
from shapely.geometry import shape, Polygon, box, Point
from shapely.geometry import MultiPolygon  # noqa F401  # flake8 issue
11
12
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
13

14
15
__author__ = "Daniel Scheffler"

16

Daniel Scheffler's avatar
Daniel Scheffler committed
17
def get_overlap_polygon(poly1, poly2, v=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
18
    """ Return a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area.
19
20
21
22
23
24
25
26

    :param poly1:   first shapely.Polygon() object
    :param poly2:   second shapely.Polygon() object
    :param v:       verbose mode
    :return:        overlap polygon as shapely.Polygon() object
    :return:        overlap percentage as float value [%]
    :return:        area of overlap polygon
    """
27
    # compute overlap polygon
28
    overlap_poly = poly1.intersection(poly2)
29
30
31
    if overlap_poly.geom_type == 'GeometryCollection':
        overlap_poly = overlap_poly.buffer(0)  # converts 'GeometryCollection' to 'MultiPolygon'

32
    if not overlap_poly.is_empty:
33
        # check if output is MultiPolygon or GeometryCollection -> if yes, convert to Polygon
34
35
        if overlap_poly.geom_type == 'MultiPolygon':
            overlap_poly = fill_holes_within_poly(overlap_poly)
36
37
        assert overlap_poly.geom_type == 'Polygon', \
            "get_overlap_polygon() did not return geometry type 'Polygon' but %s." % overlap_poly.geom_type
38

39
        overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
40
41
42
43
        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}
44
    else:
45
        return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0}
46
47


48
def get_footprint_polygon(CornerLonLat, fix_invalid=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
49
50
    """ Convert a list of coordinates into a shapely polygon object.

51
52
    :param CornerLonLat:    a list of coordinate tuples like [[lon,lat], [lon. lat], ..]
                            in clockwise or counter clockwise order
53
54
    :param fix_invalid:     fix invalid output polygon by returning its convex hull (somtimes this can be different)
    :return:                a shapely.Polygon() object
55
    """
56
57
58

    if fix_invalid:
        with warnings.catch_warnings():
59
            warnings.simplefilter('ignore')  # FIXME not working
60
61
62
63
64
65
66
67
            outpoly = Polygon(CornerLonLat)

            if not outpoly.is_valid:
                outpoly = outpoly.convex_hull
    else:
        outpoly = Polygon(CornerLonLat)

    assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \
68
                             'Got coordinates %s.' % CornerLonLat
69
70
71
    return outpoly


72
def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im, tolerance_ndigits=5):
Daniel Scheffler's avatar
Daniel Scheffler committed
73
    """Return image coordinates of the smallest box at the given coordinate grid that contains the given map coords box.
74
75
76
77
78
79
80
81

    :param box_mapYX:           input box ccordinates as YX-tuples
    :param gt_im:               geotransform of input box
    :param tolerance_ndigits:   tolerance to avoid that output image coordinates are rounded to next integer although
                                they have been very close to an integer before (this avoids float rounding issues)
                                -> tolerance is given as number of decimal digits of an image coordinate
    :return:
    """
82
83
    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
84
85
86
87
88

    # round min coords off and max coords on but tolerate differences below n decimal digits as the integer itself
    xmin, ymin, xmax, ymax = np.round([xmin, ymin, xmax, ymax], tolerance_ndigits)
    xmin, ymin, xmax, ymax = int(xmin), math.ceil(ymin), math.ceil(xmax), int(ymax)

89
    return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin)  # UL_YX,UR_YX,LR_YX,LL_YX
90
91


92
93
94
95
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')
96
97
98
99
    return box(xmin, ymin, xmax, ymax)


def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly):
Daniel Scheffler's avatar
Daniel Scheffler committed
100
    """Return the smallest box that matches the coordinate grid of the given geotransform.
101
    The returned shapely polygon contains image coordinates."""
102
    xmin, ymin, xmax, ymax = shapelyPoly.bounds  # image_coords-bounds
103

104
    return box(int(xmin), int(ymin), math.ceil(xmax), math.ceil(ymax))  # round min coords off and max coords on
105

Daniel Scheffler's avatar
Daniel Scheffler committed
106

107
108
109
110
def find_line_intersection_point(line1, line2):
    xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
    ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])

111
    def det(a, b): return a[0] * b[1] - a[1] * b[0]
112
113
114
115
116
117
    div = det(xdiff, ydiff)
    if div == 0:
        return None, None
    d = (det(*line1), det(*line2))
    x = det(d, xdiff) / div
    y = det(d, ydiff) / div
Daniel Scheffler's avatar
Daniel Scheffler committed
118
119
120
121
    return x, y


def polyVertices_outside_poly(inner_poly, outer_poly):
Daniel Scheffler's avatar
Daniel Scheffler committed
122
    """Check if a shapely polygon (inner_poly) contains vertices that do not intersect another polygon (outer_poly)
Daniel Scheffler's avatar
Daniel Scheffler committed
123
124
125
126
127
128
129
130
131
132
133
134
135
136

    :param inner_poly: the polygon with the vertices to check
    :param outer_poly: the polygon where all vertices have to be inside
    """

    if inner_poly.within(outer_poly):
        # all vertices are inside outer_poly
        return False
    elif inner_poly.intersects(outer_poly):
        # check if all vertices intersect with outer poly
        GDF = GeoDataFrame(np.swapaxes(np.array(inner_poly.exterior.coords.xy), 0, 1), columns=['X', 'Y'])
        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
137
138
139
140
        return True


def fill_holes_within_poly(poly):
Daniel Scheffler's avatar
Daniel Scheffler committed
141
    # type: (Union[Polygon, MultiPolygon]) -> Polygon
Daniel Scheffler's avatar
Daniel Scheffler committed
142
    """Fill the holes within a shapely Polygon or MultiPolygon and returns a Polygon with only the outer boundary.
143

144
    :param poly:  <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection>
145
146
    :return:
    """
147
148
149
    if poly.geom_type not in ['Polygon', 'MultiPolygon']:
        raise ValueError("Unexpected geometry type %s." % poly.geom_type)

150
    if poly.geom_type == 'Polygon':
151
        # return only the exterior polygon
152
        filled_poly = Polygon(poly.exterior)
153

154
    else:  # 'MultiPolygon'
155
156
157
158
159
160
        gdf = GeoDataFrame(columns=['geometry'])
        gdf['geometry'] = poly

        # get the area of each polygon of the multipolygon EXCLUDING the gaps in it
        gdf['area_filled'] = gdf.apply(
            lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1)
161

162
        largest_poly_filled = gdf.loc[gdf['area_filled'].idxmax()]['geometry']
163
164

        # return the outer boundary of the largest polygon
165
166
167
        filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))

    return filled_poly