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

3
4
# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
# operations when handling geospatial raster and vector data as well as projections.
5
#
6
7
8
9
# Copyright (C) 2016-2021
# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
#   Germany (https://www.gfz-potsdam.de/)
10
11
12
13
14
#
# 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).
#
15
16
17
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
26

27
import math
28
import warnings
Daniel Scheffler's avatar
Daniel Scheffler committed
29
import numpy as np
Daniel Scheffler's avatar
Daniel Scheffler committed
30
from typing import Union  # noqa F401  # flake8 issue
31

Daniel Scheffler's avatar
Daniel Scheffler committed
32
from geopandas import GeoDataFrame
33
from shapely.geometry import shape, Polygon, box
Daniel Scheffler's avatar
Daniel Scheffler committed
34
from shapely.geometry import MultiPolygon  # noqa F401  # flake8 issue
35
36
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
37
from ..coord_calc import get_corner_coordinates
38

39
40
__author__ = "Daniel Scheffler"

41

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

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

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

64
        overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area
65
        if v:
66
67
            print('%.2f percent of the image to be shifted is covered by the reference image.'
                  % overlap_percentage)  # pragma: no cover
68
69
        return {'overlap poly': overlap_poly, 'overlap percentage': overlap_percentage,
                'overlap area': overlap_poly.area}
70
    else:
71
        return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0}
72
73


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

77
78
    :param CornerLonLat:    a list of coordinate tuples like [[lon,lat], [lon. lat], ..]
                            in clockwise or counter clockwise order
79
    :param fix_invalid:     fix invalid output polygon by returning its convex hull (sometimes this can be different)
80
    :return:                a shapely.Polygon() object
81
    """
82
83
    if fix_invalid:
        with warnings.catch_warnings():
84
            warnings.simplefilter('ignore')  # FIXME not working
85
86
87
88
89
90
91
92
            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.' \
93
                             'Got coordinates %s.' % CornerLonLat
94
95
96
    return outpoly


97
98
99
100
101
102
103
104
105
106
107
def get_bounds_polygon(gt, rows, cols):
    """Get a polygon representing the outer bounds of an image.

    :param gt:      GDAL geotransform
    :param rows:    number of rows
    :param cols:    number of columns
    :return:
    """
    return get_footprint_polygon(get_corner_coordinates(gt=gt, cols=cols, rows=rows))


108
def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im, tolerance_ndigits=5):
Daniel Scheffler's avatar
Daniel Scheffler committed
109
    """Return image coordinates of the smallest box at the given coordinate grid that contains the given map coords box.
110

111
    :param box_mapYX:           input box coordinates as YX-tuples
112
113
114
115
116
117
    :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:
    """
118
119
    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
120
121
122

    # 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)
123
    xmin, ymin, xmax, ymax = math.floor(xmin), math.ceil(ymin), math.ceil(xmax), math.floor(ymax)
124

125
    return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin)  # UL_YX,UR_YX,LR_YX,LL_YX
126
127


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

133
134
135
136
    return box(xmin, ymin, xmax, ymax)


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

141
142
143
144
    # round min coords off and max coords on
    out_poly = box(math.floor(xmin), math.floor(ymin), math.ceil(xmax), math.ceil(ymax))

    return out_poly
145

Daniel Scheffler's avatar
Daniel Scheffler committed
146

147
148
149
150
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])

151
    def det(a, b): return a[0] * b[1] - a[1] * b[0]
152
153
154
155
156
157
    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
158

Daniel Scheffler's avatar
Daniel Scheffler committed
159
160
161
    return x, y


162
def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0):
163
    """Check if a shapely polygon (inner_poly) contains vertices that are outside of another polygon (outer_poly).
Daniel Scheffler's avatar
Daniel Scheffler committed
164
165
166

    :param inner_poly: the polygon with the vertices to check
    :param outer_poly: the polygon where all vertices have to be inside
167
    :param tolerance:  tolerance of the decision
Daniel Scheffler's avatar
Daniel Scheffler committed
168
    """
169
    return not outer_poly.buffer(tolerance).contains(inner_poly)
170
171
172


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

176
    :param poly:  <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection>
177
178
    :return:
    """
179
180
181
    if poly.geom_type not in ['Polygon', 'MultiPolygon']:
        raise ValueError("Unexpected geometry type %s." % poly.geom_type)

182
    if poly.geom_type == 'Polygon':
183
        # return only the exterior polygon
184
        filled_poly = Polygon(poly.exterior)
185

186
    else:  # 'MultiPolygon'
187
        gdf = GeoDataFrame(columns=['geometry'])
188
        gdf['geometry'] = poly.geoms
189
190

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

195
        largest_poly_filled = gdf.loc[gdf['area_filled'].astype(np.float64).idxmax()]['geometry']
196

197
198
199
        # noinspection PyTypeChecker
        gdf['within_or_equal'] = gdf.apply(
            lambda GDF_row:
Daniel Scheffler's avatar
Daniel Scheffler committed
200
            GDF_row.geometry.within(largest_poly_filled.buffer(1e-5)) or
201
202
203
            GDF_row.geometry.equals(largest_poly_filled),
            axis=1)

204
        if False in gdf.within_or_equal.values:
Daniel Scheffler's avatar
Daniel Scheffler committed
205
206
207
208
            n_disjunct_polys = int(np.sum(~gdf.within_or_equal))
            warnings.warn(RuntimeWarning('The given MultiPolygon contains %d disjunct polygone(s) outside of the '
                                         'largest polygone. fill_holes_within_poly() will only return the largest '
                                         'polygone as a filled version.' % n_disjunct_polys))
209

210
        # return the outer boundary of the largest polygon
211
212
213
        filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1))

    return filled_poly