Coverage for py_tools_ds/geo/vector/topology.py: 62%

1# -*- coding: utf-8 -*- 

2 

3# py_tools_ds - A collection of geospatial data analysis tools that simplify standard 

4# operations when handling geospatial raster and vector data as well as projections. 

5# 

6# Copyright (C) 2016-2021 

7# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) 

8# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam, 

9# Germany (https://www.gfz-potsdam.de/) 

10# 

11# This software was developed within the context of the GeoMultiSens project funded 

12# by the German Federal Ministry of Education and Research 

13# (project grant code: 01 IS 14 010 A-C). 

14# 

15# Licensed under the Apache License, Version 2.0 (the "License"); 

16# you may not use this file except in compliance with the License. 

17# You may obtain a copy of the License at 

18# 

19# http://www.apache.org/licenses/LICENSE-2.0 

20# 

21# Unless required by applicable law or agreed to in writing, software 

22# distributed under the License is distributed on an "AS IS" BASIS, 

23# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

24# See the License for the specific language governing permissions and 

25# limitations under the License. 

26 

27import math 

28import warnings 

29import numpy as np 

30from typing import Union # noqa F401 # flake8 issue 

31 

32from shapely.geometry import shape, Polygon, box 

33from shapely.geometry import MultiPolygon # noqa F401 # flake8 issue 

34from ..coord_trafo import mapYX2imYX 

35from ..coord_grid import find_nearest_grid_coord 

36from ..coord_calc import get_corner_coordinates 

37 

38__author__ = "Daniel Scheffler" 

39 

40 

41def get_overlap_polygon(poly1, poly2, v=False): 

42 """Return a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area. 

43 

44 :param poly1: first shapely.Polygon() object 

45 :param poly2: second shapely.Polygon() object 

46 :param v: verbose mode 

47 :return: overlap polygon as shapely.Polygon() object 

48 :return: overlap percentage as float value [%] 

49 :return: area of overlap polygon 

50 """ 

51 # compute overlap polygon 

52 overlap_poly = poly1.intersection(poly2) 

53 if overlap_poly.geom_type == 'GeometryCollection': 

54 overlap_poly = overlap_poly.buffer(0) # converts 'GeometryCollection' to 'MultiPolygon' 

55 

56 if not overlap_poly.is_empty: 

57 # check if output is MultiPolygon or GeometryCollection -> if yes, convert to Polygon 

58 if overlap_poly.geom_type == 'MultiPolygon': 

59 overlap_poly = fill_holes_within_poly(overlap_poly) 

60 assert overlap_poly.geom_type == 'Polygon', \ 

61 "get_overlap_polygon() did not return geometry type 'Polygon' but %s." % overlap_poly.geom_type 

62 

63 overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area 

64 if v: 

65 print('%.2f percent of the image to be shifted is covered by the reference image.' 

66 % overlap_percentage) # pragma: no cover 

67 return {'overlap poly': overlap_poly, 'overlap percentage': overlap_percentage, 

68 'overlap area': overlap_poly.area} 

69 else: 

70 return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0} 

71 

72 

73def get_footprint_polygon(CornerLonLat, fix_invalid=False): 

74 """Convert a list of coordinates into a shapely polygon object. 

75 

76 :param CornerLonLat: a list of coordinate tuples like [[lon,lat], [lon. lat], ..] 

77 in clockwise or counter clockwise order 

78 :param fix_invalid: fix invalid output polygon by returning its convex hull (sometimes this can be different) 

79 :return: a shapely.Polygon() object 

80 """ 

81 if fix_invalid: 

82 with warnings.catch_warnings(): 

83 warnings.simplefilter('ignore') # FIXME not working 

84 outpoly = Polygon(CornerLonLat) 

85 

86 if not outpoly.is_valid: 

87 outpoly = outpoly.convex_hull 

88 else: 

89 outpoly = Polygon(CornerLonLat) 

90 

91 assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \ 

92 'Got coordinates %s.' % CornerLonLat 

93 return outpoly 

94 

95 

96def get_bounds_polygon(gt, rows, cols): 

97 """Get a polygon representing the outer bounds of an image. 

98 

99 :param gt: GDAL geotransform 

100 :param rows: number of rows 

101 :param cols: number of columns 

102 :return: 

103 """ 

104 return get_footprint_polygon(get_corner_coordinates(gt=gt, cols=cols, rows=rows)) 

105 

106 

107def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im, tolerance_ndigits=5): 

108 """Return image coordinates of the smallest box at the given coordinate grid that contains the given map coords box. 

109 

110 :param box_mapYX: input box coordinates as YX-tuples 

111 :param gt_im: geotransform of input box 

112 :param tolerance_ndigits: tolerance to avoid that output image coordinates are rounded to next integer although 

113 they have been very close to an integer before (this avoids float rounding issues) 

114 -> tolerance is given as number of decimal digits of an image coordinate 

115 :return: 

116 """ 

117 xmin, ymin, xmax, ymax = Polygon([(i[1], i[0]) for i in box_mapYX]).bounds # map-bounds box_mapYX 

118 (ymin, xmin), (ymax, xmax) = mapYX2imYX([ymin, xmin], gt_im), mapYX2imYX([ymax, xmax], gt_im) # image coord bounds 

119 

120 # round min coords off and max coords on but tolerate differences below n decimal digits as the integer itself 

121 xmin, ymin, xmax, ymax = np.round([xmin, ymin, xmax, ymax], tolerance_ndigits) 

122 xmin, ymin, xmax, ymax = math.floor(xmin), math.ceil(ymin), math.ceil(xmax), math.floor(ymax) 

123 

124 return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin) # UL_YX,UR_YX,LR_YX,LL_YX 

125 

126 

127def get_largest_onGridPoly_within_poly(outerPoly, gt, rows, cols): 

128 oP_xmin, oP_ymin, oP_xmax, oP_ymax = outerPoly.bounds 

129 xmin, ymax = find_nearest_grid_coord((oP_xmin, oP_ymax), gt, rows, cols, direction='SE') 

130 xmax, ymin = find_nearest_grid_coord((oP_xmax, oP_ymin), gt, rows, cols, direction='NW') 

131 

132 return box(xmin, ymin, xmax, ymax) 

133 

134 

135def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly): 

136 """Return the smallest box that matches the coordinate grid of the given geotransform. 

137 The returned shapely polygon contains image coordinates.""" 

138 xmin, ymin, xmax, ymax = shapelyPoly.bounds # image_coords-bounds 

139 

140 # round min coords off and max coords on 

141 out_poly = box(math.floor(xmin), math.floor(ymin), math.ceil(xmax), math.ceil(ymax)) 

142 

143 return out_poly 

144 

145 

146def find_line_intersection_point(line1, line2): 

147 xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0]) 

148 ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1]) 

149 

150 def det(a, b): return a[0] * b[1] - a[1] * b[0] 

151 div = det(xdiff, ydiff) 

152 if div == 0: 

153 return None, None 

154 d = (det(*line1), det(*line2)) 

155 x = det(d, xdiff) / div 

156 y = det(d, ydiff) / div 

157 

158 return x, y 

159 

160 

161def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0): 

162 """Check if a shapely polygon (inner_poly) contains vertices that are outside of another polygon (outer_poly). 

163 

164 :param inner_poly: the polygon with the vertices to check 

165 :param outer_poly: the polygon where all vertices have to be inside 

166 :param tolerance: tolerance of the decision 

167 """ 

168 return not outer_poly.buffer(tolerance).contains(inner_poly) 

169 

170 

171def fill_holes_within_poly(poly: Union[Polygon, MultiPolygon] 

172 ) -> Polygon: 

173 """Fill the holes within a shapely Polygon or MultiPolygon and return a Polygon with only the outer boundary. 

174 

175 :param poly: <shapely.geometry.Polygon, shapely.geometry.MultiPolygon>, shapely.geometry.GeometryCollection> 

176 :return: 

177 """ 

178 def close_holes(polygon: Polygon) -> Polygon: 

179 if polygon.interiors: 

180 return Polygon(list(polygon.exterior.coords)) 

181 else: 

182 return polygon 

183 

184 if poly.geom_type not in ['Polygon', 'MultiPolygon']: 

185 raise ValueError("Unexpected geometry type %s." % poly.geom_type) 

186 

187 if poly.geom_type == 'Polygon': 

188 # return only the exterior polygon 

189 filled_poly = Polygon(poly.exterior) 

190 

191 else: # 'MultiPolygon' 

192 multipoly_closed = MultiPolygon([close_holes(p) for p in poly.geoms]) 

193 polys_areasorted = list(sorted(multipoly_closed.geoms, key=lambda a: a.area, reverse=True)) 

194 poly_largest = polys_areasorted[0] 

195 

196 polys_disjunct = [p for p in polys_areasorted[1:] if p.disjoint(poly_largest)] 

197 

198 if polys_disjunct: 

199 warnings.warn(RuntimeWarning('The given MultiPolygon contains %d disjunct polygon(s) outside of the ' 

200 'largest polygon. fill_holes_within_poly() will only return the largest ' 

201 'polygon as a filled version.' % len(polys_disjunct))) 

202 

203 filled_poly = poly_largest 

204 

205 return filled_poly