Coverage for gws-app/gws/base/shape/__init__.py: 58%

1"""Shape object. 

2 

3The Shape object represents a geo-referenced geometry. 

4Internally, it holds a pointer to a Shapely geometry object and a Crs object. 

5""" 

6 

7# @TODO support for SQL/MM extensions 

8 

9import struct 

10import re 

11import shapely.geometry 

12import shapely.ops 

13import shapely.wkb 

14import shapely.wkt 

15 

16import gws 

17import gws.lib.crs 

18import gws.lib.sa as sa 

19 

20_TOLERANCE_QUAD_SEGS = 6 

21_MIN_TOLERANCE_RADIUS = 0.01 

22 

23 

24class Error(gws.Error): 

25 pass 

26 

27 

28def from_wkt(wkt: str, default_crs: gws.Crs = None) -> gws.Shape: 

29 """Creates a shape object from a WKT string. 

30 

31 Args: 

32 wkt: A WKT or EWKT string. 

33 default_crs: Default Crs. 

34 

35 Returns: 

36 A Shape object. 

37 """ 

38 

39 if wkt.startswith('SRID='): 

40 # EWKT 

41 c = wkt.index(';') 

42 srid = wkt[len('SRID=') : c] 

43 crs = gws.lib.crs.require(int(srid)) 

44 wkt = wkt[c + 1 :] 

45 elif default_crs: 

46 crs = default_crs 

47 else: 

48 raise Error('missing or invalid crs for WKT') 

49 

50 return Shape(shapely.wkt.loads(wkt), crs) 

51 

52 

53def from_wkb(wkb: bytes, default_crs: gws.Crs = None) -> gws.Shape: 

54 """Creates a shape object from a WKB byte string. 

55 

56 Args: 

57 wkb: A WKB or EWKB byte string. 

58 default_crs: Default Crs. 

59 

60 Returns: 

61 A Shape object. 

62 """ 

63 

64 return _from_wkb(wkb, default_crs) 

65 

66 

67def from_wkb_hex(wkb: str, default_crs: gws.Crs = None) -> gws.Shape: 

68 """Creates a shape object from a hex-encoded WKB string. 

69 

70 Args: 

71 wkb: A hex-encoded WKB or EWKB byte string. 

72 default_crs: Default Crs. 

73 

74 Returns: 

75 A Shape object. 

76 """ 

77 

78 return _from_wkb(bytes.fromhex(wkb), default_crs) 

79 

80 

81def _from_wkb(wkb: bytes, default_crs): 

82 # http://libgeos.org/specifications/wkb/#extended-wkb 

83 

84 byte_order = wkb[0] 

85 header = struct.unpack('<cLL' if byte_order == 1 else '>cLL', wkb[:9]) 

86 

87 if header[1] & 0x20000000: 

88 crs = gws.lib.crs.require(header[2]) 

89 elif default_crs: 

90 crs = default_crs 

91 else: 

92 raise Error('missing or invalid crs for WKB') 

93 

94 geom = shapely.wkb.loads(wkb) 

95 return Shape(geom, crs) 

96 

97 

98def from_wkb_element(element: sa.geo.WKBElement, default_crs: gws.Crs = None): 

99 data = element.data 

100 if isinstance(data, str): 

101 wkb = bytes.fromhex(data) 

102 else: 

103 wkb = bytes(data) 

104 crs = gws.lib.crs.get(element.srid) 

105 return _from_wkb(wkb, crs or default_crs) 

106 

107 

108def from_geojson(geojson: dict, crs: gws.Crs, always_xy=False) -> gws.Shape: 

109 """Creates a shape object from a GeoJSON geometry dict. 

110 

111 Parses a dict as a GeoJSON geometry object (https://www.rfc-editor.org/rfc/rfc7946#section-3.1). 

112 

113 The coordinates are assumed to be in the projection order, unless ``always_xy`` is ``True``. 

114 

115 Args: 

116 geojson: A GeoJSON geometry dict 

117 crs: A Crs object. 

118 always_xy: If ``True``, coordinates are assumed to be in the XY (lon/lat) order 

119 

120 Returns: 

121 A Shape object. 

122 """ 

123 

124 geom = _shapely_shape(geojson) 

125 if crs.isYX and not always_xy: 

126 geom = _swap_xy(geom) 

127 return Shape(geom, crs) 

128 

129 

130def from_props(props: gws.Props) -> gws.Shape: 

131 """Creates a Shape from a properties object. 

132 

133 Args: 

134 props: A properties object. 

135 Returns: 

136 A Shape object. 

137 """ 

138 

139 crs = gws.lib.crs.get(props.get('crs')) 

140 if not crs: 

141 raise Error('missing or invalid crs') 

142 geom = _shapely_shape(props.get('geometry')) 

143 return Shape(geom, crs) 

144 

145 

146def from_dict(d: dict) -> gws.Shape: 

147 """Creates a Shape from a dictionary. 

148 

149 Args: 

150 d: A dictionary with the keys 'crs' and 'geometry'. 

151 Returns: 

152 A Shape object. 

153 """ 

154 

155 crs = gws.lib.crs.get(d.get('crs')) 

156 if not crs: 

157 raise Error('missing or invalid crs') 

158 geom = _shapely_shape(d.get('geometry')) 

159 return Shape(geom, crs) 

160 

161 

162def from_extent(extent: gws.Extent, crs: gws.Crs, always_xy=False) -> gws.Shape: 

163 """Creates a polygon Shape from an extent. 

164 

165 Args: 

166 extent: A hex-encoded WKB byte string. 

167 crs: A Crs object. 

168 always_xy: If ``True``, coordinates are assumed to be in the XY (lon/lat) order 

169 

170 Returns: 

171 A Shape object. 

172 """ 

173 

174 geom = shapely.geometry.box(*extent) 

175 if crs.isYX and not always_xy: 

176 geom = _swap_xy(geom) 

177 return Shape(geom, crs) 

178 

179 

180def from_bounds(bounds: gws.Bounds) -> gws.Shape: 

181 """Creates a polygon Shape from a Bounds object. 

182 

183 Args: 

184 bounds: A Bounds object. 

185 

186 Returns: 

187 A Shape object. 

188 """ 

189 

190 return Shape(shapely.geometry.box(*bounds.extent), bounds.crs) 

191 

192 

193def from_xy(x: float, y: float, crs: gws.Crs) -> gws.Shape: 

194 """Creates a point Shape from coordinates. 

195 

196 Args: 

197 x: X coordinate (lon/easting) 

198 y: Y coordinate (lat/northing) 

199 crs: A Crs object. 

200 

201 Returns: 

202 A Shape object. 

203 """ 

204 

205 return Shape(shapely.geometry.Point(x, y), crs) 

206 

207 

208def _swap_xy(geom): 

209 def f(x: float, y: float, z: float = None) -> tuple[float, float]: 

210 return y, x 

211 

212 return shapely.ops.transform(f, geom) 

213 

214 

215_CIRCLE_RESOLUTION = 64 

216 

217 

218def _shapely_shape(d): 

219 if d.get('type').upper() == 'CIRCLE': 

220 geom = shapely.geometry.Point(d.get('center')) 

221 return geom.buffer( 

222 d.get('radius'), 

223 resolution=_CIRCLE_RESOLUTION, 

224 cap_style=shapely.geometry.CAP_STYLE.round, 

225 join_style=shapely.geometry.JOIN_STYLE.round, 

226 ) 

227 

228 return shapely.geometry.shape(d) 

229 

230 

231## 

232 

233 

234class Props(gws.Props): 

235 """Shape properties object.""" 

236 

237 crs: str 

238 geometry: dict 

239 

240 

241## 

242 

243 

244class Shape(gws.Shape): 

245 geom: shapely.geometry.base.BaseGeometry 

246 

247 def __init__(self, geom, crs: gws.Crs): 

248 super().__init__() 

249 self.geom = geom 

250 self.crs = crs 

251 self.type = self.geom.geom_type.lower() 

252 self.x = getattr(self.geom, 'x', None) 

253 self.y = getattr(self.geom, 'y', None) 

254 

255 def __str__(self): 

256 return '{Geometry:' + self.geom.geom_type.upper() + '}' 

257 

258 def area(self): 

259 return getattr(self.geom, 'area', 0) 

260 

261 def bounds(self): 

262 return gws.Bounds(crs=self.crs, extent=self.geom.bounds) 

263 

264 def centroid(self): 

265 return Shape(self.geom.centroid, self.crs) 

266 

267 def to_wkb(self): 

268 return shapely.wkb.dumps(self.geom) 

269 

270 def to_wkb_hex(self): 

271 return shapely.wkb.dumps(self.geom, hex=True) 

272 

273 def to_ewkb(self): 

274 return shapely.wkb.dumps(self.geom, srid=self.crs.srid) 

275 

276 def to_ewkb_hex(self): 

277 return shapely.wkb.dumps(self.geom, srid=self.crs.srid, hex=True) 

278 

279 def to_wkt(self, trim=False, rounding_precision=-1, output_dimension=3): 

280 s = shapely.wkt.dumps(self.geom, trim=trim, rounding_precision=rounding_precision, output_dimension=output_dimension) 

281 s = re.sub(r'\s*([,()])\s*', r'\1', s) 

282 s = re.sub(r'\s+', ' ', s.strip()) 

283 return s 

284 

285 def to_ewkt(self, trim=False, rounding_precision=-1, output_dimension=3): 

286 return f'SRID={self.crs.srid};' + self.to_wkt(trim=trim, rounding_precision=rounding_precision, output_dimension=output_dimension) 

287 

288 def to_geojson(self, keep_crs=False): 

289 # see https://datatracker.ietf.org/doc/html/rfc7946#section-4 

290 # convert to WGS lon,lat unless keep_crs is true 

291 # coords order is always XY 

292 

293 if keep_crs or self.crs == gws.lib.crs.WGS84: 

294 return shapely.geometry.mapping(self.geom) 

295 

296 tr = self.crs.transformer(gws.lib.crs.WGS84) 

297 new_geom = shapely.ops.transform(tr, self.geom) 

298 return shapely.geometry.mapping(new_geom) 

299 

300 def to_props(self): 

301 return gws.ShapeProps(crs=self.crs.epsg, geometry=shapely.geometry.mapping(self.geom)) 

302 

303 def is_empty(self): 

304 return self.geom.is_empty 

305 

306 def is_ring(self): 

307 return self.geom.is_ring 

308 

309 def is_simple(self): 

310 return self.geom.is_simple 

311 

312 def is_valid(self): 

313 return self.geom.is_valid 

314 

315 def equals(self, other): 

316 return self._binary_predicate(other, 'equals') 

317 

318 def contains(self, other): 

319 return self._binary_predicate(other, 'contains') 

320 

321 def covers(self, other): 

322 return self._binary_predicate(other, 'covers') 

323 

324 def covered_by(self, other): 

325 return self._binary_predicate(other, 'covered_by') 

326 

327 def crosses(self, other): 

328 return self._binary_predicate(other, 'crosses') 

329 

330 def disjoint(self, other): 

331 return self._binary_predicate(other, 'disjoint') 

332 

333 def intersects(self, other): 

334 return self._binary_predicate(other, 'intersects') 

335 

336 def overlaps(self, other): 

337 return self._binary_predicate(other, 'overlaps') 

338 

339 def touches(self, other): 

340 return self._binary_predicate(other, 'touches') 

341 

342 def within(self, other): 

343 return self._binary_predicate(other, 'within') 

344 

345 def _binary_predicate(self, other, op): 

346 s = other.transformed_to(self.crs) 

347 return getattr(self.geom, op)(getattr(s, 'geom')) 

348 

349 def union(self, others): 

350 if not others: 

351 return self 

352 

353 geoms = [self.geom] 

354 for s in others: 

355 s = s.transformed_to(self.crs) 

356 geoms.append(getattr(s, 'geom')) 

357 

358 geom = shapely.ops.unary_union(geoms) 

359 return Shape(geom, self.crs) 

360 

361 def intersection(self, *others): 

362 if not others: 

363 return self 

364 

365 geom = self.geom 

366 for s in others: 

367 s = s.transformed_to(self.crs) 

368 geom = geom.intersection(getattr(s, 'geom')) 

369 

370 return Shape(geom, self.crs) 

371 

372 def to_multi(self): 

373 if self.type == gws.GeometryType.point: 

374 return Shape(shapely.geometry.MultiPoint([self.geom]), self.crs) 

375 if self.type == gws.GeometryType.linestring: 

376 return Shape(shapely.geometry.MultiLineString([self.geom]), self.crs) 

377 if self.type == gws.GeometryType.polygon: 

378 return Shape(shapely.geometry.MultiPolygon([self.geom]), self.crs) 

379 return self 

380 

381 def to_type(self, new_type: gws.GeometryType): 

382 if new_type == self.type: 

383 return self 

384 if new_type == gws.GeometryType.geometry: 

385 return self 

386 if self.type == gws.GeometryType.point and new_type == gws.GeometryType.multipoint: 

387 return self.to_multi() 

388 if self.type == gws.GeometryType.linestring and new_type == gws.GeometryType.multilinestring: 

389 return self.to_multi() 

390 if self.type == gws.GeometryType.polygon and new_type == gws.GeometryType.multipolygon: 

391 return self.to_multi() 

392 raise Error(f'cannot convert {self.type!r} to {new_type!r}') 

393 

394 def to_2d(self): 

395 geom = shapely.force_2d(self.geom) 

396 if geom is self.geom: 

397 return self 

398 return Shape(geom, self.crs) 

399 

400 def tolerance_polygon(self, tolerance=None, quad_segs=None): 

401 is_poly = self.type in (gws.GeometryType.polygon, gws.GeometryType.multipolygon) 

402 

403 if not tolerance and is_poly: 

404 return self 

405 

406 # we need a polygon even if tolerance = 0 

407 tolerance = tolerance or _MIN_TOLERANCE_RADIUS 

408 quad_segs = quad_segs or _TOLERANCE_QUAD_SEGS 

409 

410 geom = self.geom 

411 if self.crs.isGeographic: 

412 tr = self.crs.transformer(gws.lib.crs.WEBMERCATOR) 

413 geom = shapely.ops.transform(tr, self.geom) 

414 

415 if is_poly: 

416 cs = shapely.geometry.CAP_STYLE.flat 

417 js = shapely.geometry.JOIN_STYLE.mitre 

418 else: 

419 cs = shapely.geometry.CAP_STYLE.round 

420 js = shapely.geometry.JOIN_STYLE.round 

421 

422 geom = geom.buffer(tolerance, quad_segs, cap_style=cs, join_style=js) 

423 if self.crs.isGeographic: 

424 tr = gws.lib.crs.WEBMERCATOR.transformer(self.crs) 

425 geom = shapely.ops.transform(tr, geom) 

426 

427 return Shape(geom, self.crs) 

428 

429 def transformed_to(self, crs): 

430 if crs == self.crs: 

431 return self 

432 tr = self.crs.transformer(crs) 

433 dg = shapely.ops.transform(tr, self.geom) 

434 return Shape(dg, crs)