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

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 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/>.

24
import math
25
import warnings
Daniel Scheffler's avatar
Daniel Scheffler committed
26
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
27
from typing import Union  # noqa F401  # flake8 issue
28

Daniel Scheffler's avatar
Daniel Scheffler committed
29
from geopandas import GeoDataFrame
Daniel Scheffler's avatar
Daniel Scheffler committed
30
31
from shapely.geometry import shape, Polygon, box, Point
from shapely.geometry import MultiPolygon  # noqa F401  # flake8 issue
32
33
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
34

35
36
__author__ = "Daniel Scheffler"

37

Daniel Scheffler's avatar
Daniel Scheffler committed
38
def get_overlap_polygon(poly1, poly2, v=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
39
    """ Return a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area.
40
41
42
43
44
45
46
47

    :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
    """
48
    # compute overlap polygon
49
    overlap_poly = poly1.intersection(poly2)
50
51
52
    if overlap_poly.geom_type == 'GeometryCollection':
        overlap_poly = overlap_poly.buffer(0)  # converts 'GeometryCollection' to 'MultiPolygon'

53
    if not overlap_poly.is_empty:
54
        # check if output is MultiPolygon or GeometryCollection -> if yes, convert to Polygon
55
56
        if overlap_poly.geom_type == 'MultiPolygon':
            overlap_poly = fill_holes_within_poly(overlap_poly)
57
58
        assert overlap_poly.geom_type == 'Polygon', \
            "get_overlap_polygon() did not return geometry type 'Polygon' but %s." % overlap_poly.geom_type
59

60
        overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
61
62
63
64
        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}
65
    else:
66
        return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0}
67
68


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

72
73
    :param CornerLonLat:    a list of coordinate tuples like [[lon,lat], [lon. lat], ..]
                            in clockwise or counter clockwise order
74
75
    :param fix_invalid:     fix invalid output polygon by returning its convex hull (somtimes this can be different)
    :return:                a shapely.Polygon() object
76
    """
77
78
79

    if fix_invalid:
        with warnings.catch_warnings():
80
            warnings.simplefilter('ignore')  # FIXME not working
81
82
83
84
85
86
87
88
            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.' \
89
                             'Got coordinates %s.' % CornerLonLat
90
91
92
    return outpoly


93
def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im, tolerance_ndigits=5):
Daniel Scheffler's avatar
Daniel Scheffler committed
94
    """Return image coordinates of the smallest box at the given coordinate grid that contains the given map coords box.
95
96
97
98
99
100
101
102

    :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:
    """
103
104
    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
105
106
107
108
109

    # 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)

110
    return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin)  # UL_YX,UR_YX,LR_YX,LL_YX
111
112


113
114
115
116
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')
117
118
119
120
    return box(xmin, ymin, xmax, ymax)


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

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

Daniel Scheffler's avatar
Daniel Scheffler committed
127

128
129
130
131
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])

132
    def det(a, b): return a[0] * b[1] - a[1] * b[0]
133
134
135
136
137
138
    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
139
140
141
142
    return x, y


def polyVertices_outside_poly(inner_poly, outer_poly):
Daniel Scheffler's avatar
Daniel Scheffler committed
143
    """Check if a shapely polygon (inner_poly) contains vertices that do not intersect another polygon (outer_poly)
Daniel Scheffler's avatar
Daniel Scheffler committed
144
145
146
147
148
149
150
151
152
153
154
155
156
157

    :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
158
159
160
161
        return True


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

165
    :param poly:  <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection>
166
167
    :return:
    """
168
169
170
    if poly.geom_type not in ['Polygon', 'MultiPolygon']:
        raise ValueError("Unexpected geometry type %s." % poly.geom_type)

171
    if poly.geom_type == 'Polygon':
172
        # return only the exterior polygon
173
        filled_poly = Polygon(poly.exterior)
174

175
    else:  # 'MultiPolygon'
176
177
178
179
180
181
        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)
182

183
        largest_poly_filled = gdf.loc[gdf['area_filled'].idxmax()]['geometry']
184
185

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

    return filled_poly