Coverage for py_tools_ds/geo/vector/geometry.py: 98%

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 

27from typing import Tuple 

28 

29from shapely.geometry import Polygon, box 

30import numpy as np 

31 

32from .conversion import shapelyImPoly_to_shapelyMapPoly, get_boxImXY_from_shapelyPoly, \ 

33 shapelyBox2BoxYX, round_shapelyPoly_coords 

34from ..coord_calc import corner_coord_to_minmax 

35from ..coord_trafo import imYX2mapYX, transform_coordArray 

36from ..projection import prj_equal 

37 

38__author__ = "Daniel Scheffler" 

39 

40 

41class boxObj(object): 

42 def __init__(self, **kwargs): 

43 """Create a dynamic/self-updating box object that represents a rectangular or quadratic coordinate box 

44 according to the given keyword arguments. 

45 

46 Note: Either mapPoly+gt or imPoly+gt or wp+ws or boxMapYX+gt or boxImYX+gt must be passed. 

47 

48 :Keyword Arguments: 

49 - gt (tuple): GDAL geotransform (default: (0, 1, 0, 0, 0, -1)) 

50 - prj (str): projection as WKT string 

51 - mapPoly (shapely.Polygon): Polygon with map unit vertices 

52 - imPoly (shapely.Polygon): Polygon with image unit vertices 

53 - wp (tuple): window position in map units (x,y) 

54 - ws (tuple): window size in map units (x,y) 

55 - boxMapYX (list): box map coordinates like [(ULy,ULx), (URy,URx), (LRy,LRx), (LLy,LLx)] 

56 - boxImYX (list): box image coordinates like [(ULy,ULx), (URy,URx), (LRy,LRx), (LLy,LLx)] 

57 """ 

58 # FIXME self.prj is not used 

59 # TODO allow boxObj to be instanced with gt, prj + rows/cols 

60 self.gt = kwargs.get('gt', (0, 1, 0, 0, 0, -1)) 

61 self.prj = kwargs.get('prj', '') 

62 self._mapPoly = kwargs.get('mapPoly', None) 

63 self._imPoly = kwargs.get('imPoly', None) 

64 self.wp = kwargs.get('wp', None) 

65 self._ws = kwargs.get('ws', None) 

66 self._boxMapYX = kwargs.get('boxMapYX', None) 

67 self._boxImYX = kwargs.get('boxImYX', None) 

68 

69 if self._mapPoly: 

70 if self.gt == (0, 1, 0, 0, 0, -1): 

71 raise ValueError(self.gt, "A geotransform must be passed if mapPoly is given.") 

72 else: 

73 # populate self._mapPoly 

74 if self._imPoly: 

75 self._mapPoly = shapelyImPoly_to_shapelyMapPoly(self._imPoly, self.gt) 

76 elif self._boxMapYX: 

77 self.boxMapYX = self._boxMapYX 

78 elif self._boxImYX: 

79 self.boxImYX = self._boxImYX 

80 elif self.wp or self._ws: # asserts wp/ws in map units 

81 assert self.wp and self._ws, \ 

82 "'wp' and 'ws' must be passed together. Got wp=%s and ws=%s." % (self.wp, self._ws) 

83 (wpx, wpy), (wsx, wsy) = self.wp, self._ws 

84 self._mapPoly = box(wpx - wsx / 2, wpy - wsy / 2, wpx + wsx / 2, wpy + wsy / 2) 

85 else: 

86 raise ValueError("No proper set of arguments received.") 

87 

88 # all getters and setters synchronize using self._mapPoly 

89 @property 

90 def mapPoly(self): 

91 imPoly = Polygon(get_boxImXY_from_shapelyPoly(self._mapPoly, self.gt)) 

92 return shapelyImPoly_to_shapelyMapPoly(imPoly, self.gt) 

93 

94 @mapPoly.setter 

95 def mapPoly(self, shapelyPoly): 

96 self._mapPoly = shapelyPoly 

97 

98 @property 

99 def imPoly(self): 

100 return Polygon(get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt)) 

101 

102 @imPoly.setter 

103 def imPoly(self, shapelyImPoly): 

104 self.mapPoly = shapelyImPoly_to_shapelyMapPoly(shapelyImPoly, self.gt) 

105 

106 @property 

107 def boxMapYX(self): 

108 """Return a list of YX coordinate tuples for all corners in the order UL_YX, UR_YX, LR_YX, LL_YX. 

109 

110 :return: UL_YX, UR_YX, LR_YX, LL_YX 

111 """ 

112 return shapelyBox2BoxYX(self.mapPoly, coord_type='map') 

113 

114 @boxMapYX.setter 

115 def boxMapYX(self, mapBoxYX): 

116 mapBoxXY = [(i[1], i[0]) for i in mapBoxYX] 

117 xmin, xmax, ymin, ymax = corner_coord_to_minmax(mapBoxXY) 

118 self.mapPoly = box(xmin, ymin, xmax, ymax) 

119 

120 @property 

121 def boxMapXY(self): 

122 """Return a list of XY coordinate tuples for all corners in the order UL_XY, UR_XY, LR_XY, LL_XY. 

123 

124 :return: UL_XY, UR_XY, LR_XY, LL_XY 

125 """ 

126 return tuple((x, y) for y, x in self.boxMapYX) 

127 

128 @boxMapXY.setter 

129 def boxMapXY(self, mapBoxXY): 

130 xmin, xmax, ymin, ymax = corner_coord_to_minmax(mapBoxXY) 

131 self.mapPoly = box(xmin, ymin, xmax, ymax) 

132 

133 @property 

134 def boxImYX(self): 

135 temp_imPoly = round_shapelyPoly_coords(self.imPoly, precision=0) 

136 floatImBoxYX = shapelyBox2BoxYX(temp_imPoly, coord_type='image') 

137 return [[int(i[0]), int(i[1])] for i in floatImBoxYX] 

138 

139 @boxImYX.setter 

140 def boxImYX(self, imBoxYX): 

141 imBoxXY = [(i[1], i[0]) for i in imBoxYX] 

142 xmin, xmax, ymin, ymax = corner_coord_to_minmax(imBoxXY) 

143 self.imPoly = box(xmin, ymin, xmax, ymax) 

144 

145 @property 

146 def boxImXY(self): 

147 return tuple((x, y) for y, x in self.boxImYX) 

148 

149 @boxImXY.setter 

150 def boxImXY(self, imBoxXY): 

151 xmin, xmax, ymin, ymax = corner_coord_to_minmax(imBoxXY) 

152 self.imPoly = box(xmin, ymin, xmax, ymax) 

153 

154 @property 

155 def boundsMap(self): 

156 """Return xmin,xmax,ymin,ymax in map coordinates.""" 

157 boxMapYX = shapelyBox2BoxYX(self.mapPoly, coord_type='image') 

158 boxMapXY = [(i[1], i[0]) for i in boxMapYX] 

159 return corner_coord_to_minmax(boxMapXY) 

160 

161 @property 

162 def boundsIm(self): 

163 """Return xmin,xmax,ymin,ymax in image coordinates.""" 

164 boxImXY = [(round(x, 0), round(y, 0)) for x, y in get_boxImXY_from_shapelyPoly(self.mapPoly, self.gt)] 

165 return corner_coord_to_minmax(boxImXY) # xmin,xmax,ymin,ymax 

166 

167 @property 

168 def imDimsYX(self): 

169 xmin, xmax, ymin, ymax = self.boundsIm 

170 return (ymax - ymin), (xmax - xmin) 

171 

172 @property 

173 def imDimsXY(self): 

174 return self.imDimsYX[1], self.imDimsYX[0] 

175 

176 @property 

177 def mapDimsYX(self): 

178 xmin, xmax, ymin, ymax = self.boundsMap 

179 return (ymax - ymin), (xmax - xmin) 

180 

181 @property 

182 def mapDimsXY(self): 

183 return self.mapDimsYX[1], self.mapDimsYX[0] 

184 

185 def buffer_imXY(self, buffImX=0, buffImY=0): 

186 # type: (float, float) -> None 

187 """Buffer the box in X- and/or Y-direction. 

188 

189 :param buffImX: <float> buffer value in x-direction as IMAGE UNITS (pixels) 

190 :param buffImY: <float> buffer value in y-direction as IMAGE UNITS (pixels) 

191 """ 

192 xmin, xmax, ymin, ymax = self.boundsIm 

193 xmin, xmax, ymin, ymax = xmin - buffImX, xmax + buffImX, ymin - buffImY, ymax + buffImY 

194 self.imPoly = box(xmin, ymin, xmax, ymax) 

195 

196 def buffer_mapXY(self, buffMapX=0, buffMapY=0): 

197 # type: (float, float) -> None 

198 """Buffer the box in X- and/or Y-direction. 

199 

200 :param buffMapX: <float> buffer value in x-direction as MAP UNITS 

201 :param buffMapY: <float> buffer value in y-direction as MAP UNITS 

202 """ 

203 xmin, xmax, ymin, ymax = self.boundsMap 

204 xmin, xmax, ymin, ymax = xmin - buffMapX, xmax + buffMapX, ymin - buffMapY, ymax + buffMapY 

205 self.mapPoly = box(xmin, ymin, xmax, ymax) 

206 

207 def is_larger_DimXY(self, boundsIm2test): 

208 # type: (tuple) -> (bool,bool) 

209 """Check if the boxObj is larger than a given set of bounding image coordinates (in X- and/or Y-direction). 

210 

211 :param boundsIm2test: <tuple> (xmin,xmax,ymin,ymax) as image coordinates 

212 """ 

213 b2t_xmin, b2t_xmax, b2t_ymin, b2t_ymax = boundsIm2test 

214 xmin, xmax, ymin, ymax = self.boundsIm 

215 x_is_larger = xmin < b2t_xmin or xmax > b2t_xmax 

216 y_is_larger = ymin < b2t_ymin or ymax > b2t_ymax 

217 return x_is_larger, y_is_larger 

218 

219 def get_coordArray_MapXY(self, prj=None): 

220 # type: (str) -> Tuple[np.ndarray] 

221 """Return two coordinate arrays for X and Y coordinates in the given projection. 

222 

223 NOTE: If no projection is given, <boxObj>.prj is used. 

224 

225 :param prj: GDAL projection as WKT string 

226 :return: 

227 """ 

228 if not (self.gt and self.prj): 

229 raise ValueError('geotransform and projection must both be available for computing coordinate array.' 

230 'boxObj.gt=%s; boxobj.prj=%s' % (self.gt, self.prj)) 

231 

232 xmin, xmax, ymin, ymax = self.boundsMap 

233 Xarr, Yarr = np.meshgrid(np.arange(xmin, xmax, abs(self.gt[1])), 

234 np.arange(ymax, ymin, -abs(self.gt[5]))) 

235 if prj and not prj_equal(self.prj, prj): 

236 Xarr, Yarr = transform_coordArray(self.prj, prj, Xarr, Yarr) 

237 

238 return Xarr, Yarr 

239 

240 

241def get_winPoly(wp_imYX, ws, gt, match_grid=0): 

242 # type: (tuple, tuple, list, bool) -> (Polygon, tuple, tuple) 

243 """Create a shapely polygon from a given set of image cordinates, window size and geotransform. 

244 

245 :param wp_imYX: 

246 :param ws: <tuple> X/Y window size 

247 :param gt: 

248 :param match_grid: <bool> 

249 """ 

250 ws = (ws, ws) if isinstance(ws, int) else ws 

251 xmin, xmax, ymin, ymax = ( 

252 wp_imYX[1] - ws[0] / 2, wp_imYX[1] + ws[0] / 2, wp_imYX[0] + ws[1] / 2, wp_imYX[0] - ws[1] / 2) 

253 if match_grid: 

254 xmin, xmax, ymin, ymax = [int(i) for i in [xmin, xmax, ymin, ymax]] 

255 box_YX = (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin) # UL,UR,LR,LL 

256 UL, UR, LR, LL = [imYX2mapYX(imYX, gt) for imYX in box_YX] 

257 box_mapYX = UL, UR, LR, LL 

258 return Polygon([(UL[1], UL[0]), (UR[1], UR[0]), (LR[1], LR[0]), (LL[1], LR[0])]), box_mapYX, box_YX