Coverage for gws-app/gws/lib/shape/__init__.py: 56%

1"""Shapes. 

2 

3A shape (``gws.Shape``) is a geo-referenced geometry: a Shapely geometry 

4together with a ``gws.Crs``. Shapes are used for feature geometries, search 

5geometries and map extents throughout the application. 

6 

7The package is a single module. It provides: 

8 

9- constructors that create a ``Shape`` from WKT and EWKT, WKB and EWKB (binary or 

10 hex), SQLAlchemy/GeoAlchemy WKB elements, GeoJSON geometries, shape props or 

11 dicts, extents, ``gws.Bounds`` and x/y coordinates, 

12- the ``Shape`` class, which implements the ``gws.Shape`` interface: conversion to 

13 WKB, WKT, GeoJSON and props, spatial predicates, union and intersection, type 

14 conversions, buffering with a tolerance and transformation to other CRS, 

15- the ``Props`` class for shapes sent to and from the client. 

16 

17Constructors raise ``Error`` if the input cannot be parsed or has no CRS. 

18EWKT and EWKB inputs carry their own SRID; for plain WKT and WKB a default CRS 

19must be given. GeoJSON inputs and extents are expected in the axis order of the 

20CRS, unless ``always_xy`` is set. Circles (``{"type": "Circle", "center": ..., 

21"radius": ...}``), as sent by the client, are converted to polygons. 

22 

23Binary predicates and set operations transform the other shape to the CRS of 

24this shape first. ``to_geojson`` transforms to WGS84 unless asked to keep the CRS. 

25 

26Example:: 

27 

28 import gws.lib.shape 

29 import gws.lib.crs 

30 

31 shape = gws.lib.shape.from_wkt('POINT(10 20)', gws.lib.crs.WGS84) 

32 area = shape.tolerance_polygon(5).transformed_to(gws.lib.crs.WEBMERCATOR) 

33 ewkt = area.to_ewkt() 

34 

35 other = gws.lib.shape.from_wkt('SRID=4326;POLYGON((0 0,30 0,30 30,0 30,0 0))') 

36 print(other.contains(shape)) 

37""" 

38 

39# @TODO support for SQL/MM extensions 

40 

41import struct 

42import re 

43import shapely.errors 

44import shapely.geometry 

45import shapely.ops 

46import shapely.wkb 

47import shapely.wkt 

48 

49import gws 

50import gws.lib.crs 

51import gws.lib.sa as sa 

52 

53_TOLERANCE_QUAD_SEGS = 6 

54_MIN_TOLERANCE_RADIUS = 0.01 

55 

56 

57class Error(gws.Error): 

58 """Invalid geometry or CRS.""" 

59 

60 pass 

61 

62 

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

64 """Create a shape from a WKT or EWKT string. 

65 

66 Args: 

67 wkt: A WKT or EWKT string. 

68 default_crs: CRS to use if the string has no SRID. 

69 

70 Returns: 

71 A Shape object. 

72 

73 Raises: 

74 ``Error``: If the string is invalid or there is no CRS. 

75 """ 

76 

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

78 # EWKT 

79 c = wkt.index(';') 

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

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

82 wkt = wkt[c + 1 :] 

83 elif default_crs: 

84 crs = default_crs 

85 else: 

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

87 

88 try: 

89 geom = shapely.wkt.loads(wkt) 

90 except shapely.errors.ShapelyError as exc: 

91 raise Error('invalid WKT') from exc 

92 return Shape(geom, crs) 

93 

94 

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

96 """Create a shape from a WKB or EWKB byte string. 

97 

98 Args: 

99 wkb: A WKB or EWKB byte string. 

100 default_crs: CRS to use if the data has no SRID. 

101 

102 Returns: 

103 A Shape object. 

104 

105 Raises: 

106 ``Error``: If the data is invalid or there is no CRS. 

107 """ 

108 

109 return _from_wkb(wkb, default_crs) 

110 

111 

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

113 """Create a shape from a hex-encoded WKB or EWKB string. 

114 

115 Args: 

116 wkb: A hex-encoded WKB or EWKB string. 

117 default_crs: CRS to use if the data has no SRID. 

118 

119 Returns: 

120 A Shape object. 

121 

122 Raises: 

123 ``Error``: If the data is invalid or there is no CRS. 

124 """ 

125 

126 try: 

127 b = bytes.fromhex(wkb) 

128 except ValueError as exc: 

129 raise Error('invalid WKB hex') from exc 

130 return _from_wkb(b, default_crs) 

131 

132 

133def _from_wkb(wkb: bytes, default_crs): 

134 """Create a shape from WKB or EWKB bytes, reading the SRID from the EWKB header.""" 

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

136 

137 try: 

138 byte_order = wkb[0] 

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

140 except (IndexError, struct.error) as exc: 

141 raise Error('invalid WKB') from exc 

142 

143 if header[1] & 0x20000000: 

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

145 elif default_crs: 

146 crs = default_crs 

147 else: 

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

149 

150 try: 

151 geom = shapely.wkb.loads(wkb) 

152 except shapely.errors.ShapelyError as exc: 

153 raise Error('invalid WKB') from exc 

154 return Shape(geom, crs) 

155 

156 

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

158 """Create a shape from a GeoAlchemy WKB element. 

159 

160 The CRS is taken from the EWKB data, then from the SRID of the element, then 

161 from ``default_crs``. 

162 

163 Args: 

164 element: A WKB element, with binary or hex-encoded data. 

165 default_crs: CRS to use if neither the data nor the element has a valid SRID. 

166 

167 Returns: 

168 A Shape object. 

169 

170 Raises: 

171 ``Error``: If the data is invalid or there is no CRS. 

172 """ 

173 data = element.data 

174 if isinstance(data, str): 

175 wkb = bytes.fromhex(data) 

176 else: 

177 wkb = bytes(data) 

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

179 return _from_wkb(wkb, crs or default_crs) 

180 

181 

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

183 """Create a shape from a GeoJSON geometry dict. 

184 

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

186 A ``Circle`` geometry with ``center`` and ``radius`` is converted to a polygon. 

187 The coordinates are assumed to be in the axis order of the CRS, unless ``always_xy`` is ``True``. 

188 

189 Args: 

190 geojson: A GeoJSON geometry dict. 

191 crs: A Crs object. 

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

193 

194 Returns: 

195 A Shape object. 

196 

197 Raises: 

198 ``Error``: If the geometry is invalid. 

199 """ 

200 

201 geom = _shapely_shape(geojson) 

202 if crs.isYX and not always_xy: 

203 geom = _swap_xy(geom) 

204 return Shape(geom, crs) 

205 

206 

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

208 """Create a shape from a properties object. 

209 

210 Args: 

211 props: A properties object with ``crs`` and ``geometry`` (a GeoJSON geometry dict). 

212 

213 Returns: 

214 A Shape object. 

215 

216 Raises: 

217 ``Error``: If the CRS or the geometry is invalid. 

218 """ 

219 

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

221 if not crs: 

222 raise Error('missing or invalid crs') 

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

224 return Shape(geom, crs) 

225 

226 

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

228 """Create a shape from a dictionary. 

229 

230 Args: 

231 d: A dictionary with the keys ``crs`` and ``geometry`` (a GeoJSON geometry dict). 

232 

233 Returns: 

234 A Shape object. 

235 

236 Raises: 

237 ``Error``: If the CRS or the geometry is invalid. 

238 """ 

239 

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

241 if not crs: 

242 raise Error('missing or invalid crs') 

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

244 return Shape(geom, crs) 

245 

246 

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

248 """Create a polygon shape from an extent. 

249 

250 Args: 

251 extent: An extent. 

252 crs: A Crs object. 

253 always_xy: If ``True``, the extent is assumed to be in the XY (lon/lat) order, 

254 otherwise in the axis order of the CRS. 

255 

256 Returns: 

257 A Shape object. 

258 """ 

259 

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

261 if crs.isYX and not always_xy: 

262 geom = _swap_xy(geom) 

263 return Shape(geom, crs) 

264 

265 

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

267 """Create a polygon shape from a Bounds object. 

268 

269 Args: 

270 bounds: A Bounds object. 

271 

272 Returns: 

273 A Shape object. 

274 """ 

275 

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

277 

278 

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

280 """Create a point shape from coordinates. 

281 

282 Args: 

283 x: X coordinate (lon/easting). 

284 y: Y coordinate (lat/northing). 

285 crs: A Crs object. 

286 

287 Returns: 

288 A Shape object. 

289 """ 

290 

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

292 

293 

294def _swap_xy(geom): 

295 """Return a copy of a Shapely geometry with x and y swapped.""" 

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

297 return y, x 

298 

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

300 

301 

302_CIRCLE_RESOLUTION = 64 

303 

304 

305def _shapely_shape(d): 

306 """Create a Shapely geometry from a GeoJSON dict, raising ``Error`` if it is invalid.""" 

307 try: 

308 return _shapely_shape2(d) 

309 except (shapely.errors.ShapelyError, AttributeError, TypeError, ValueError) as exc: 

310 raise Error('invalid geometry') from exc 

311 

312 

313def _shapely_shape2(d): 

314 """Create a Shapely geometry from a GeoJSON dict, converting circles to polygons.""" 

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

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

317 return geom.buffer( 

318 d.get('radius'), 

319 resolution=_CIRCLE_RESOLUTION, 

320 cap_style=shapely.geometry.CAP_STYLE.round, 

321 join_style=shapely.geometry.JOIN_STYLE.round, 

322 ) 

323 

324 return shapely.geometry.shape(d) 

325 

326 

327## 

328 

329 

330class Props(gws.Props): 

331 """Shape properties object.""" 

332 

333 crs: str 

334 geometry: dict 

335 

336 

337## 

338 

339 

340class Shape(gws.Shape): 

341 """Shape implemented with a Shapely geometry.""" 

342 

343 geom: shapely.geometry.base.BaseGeometry 

344 """Shapely geometry.""" 

345 

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

347 """Create a shape. 

348 

349 Args: 

350 geom: Shapely geometry. 

351 crs: CRS of the geometry. 

352 """ 

353 super().__init__() 

354 self.geom = geom 

355 self.crs = crs 

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

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

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

359 

360 def __str__(self): 

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

362 

363 def area(self): 

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

365 

366 def bounds(self): 

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

368 

369 def centroid(self): 

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

371 

372 def center(self): 

373 c = self.geom.centroid 

374 return c.x, c.y 

375 

376 def to_wkb(self): 

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

378 

379 def to_wkb_hex(self): 

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

381 

382 def to_ewkb(self): 

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

384 

385 def to_ewkb_hex(self): 

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

387 

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

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

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

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

392 return s 

393 

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

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

396 

397 def to_geojson(self, keep_crs=False): 

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

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

400 # coords order is always XY 

401 

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

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

404 

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

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

407 return shapely.geometry.mapping(new_geom) 

408 

409 def to_precision(self, prec: int): 

410 geom = shapely.set_precision(self.geom, 10 ** -prec) 

411 return Shape(geom, self.crs) 

412 

413 def to_props(self): 

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

415 

416 def is_empty(self): 

417 return self.geom.is_empty 

418 

419 def is_ring(self): 

420 return self.geom.is_ring 

421 

422 def is_simple(self): 

423 return self.geom.is_simple 

424 

425 def is_valid(self): 

426 return self.geom.is_valid 

427 

428 def equals(self, other): 

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

430 

431 def contains(self, other): 

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

433 

434 def covers(self, other): 

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

436 

437 def covered_by(self, other): 

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

439 

440 def crosses(self, other): 

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

442 

443 def disjoint(self, other): 

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

445 

446 def intersects(self, other): 

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

448 

449 def overlaps(self, other): 

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

451 

452 def touches(self, other): 

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

454 

455 def within(self, other): 

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

457 

458 def _binary_predicate(self, other, op): 

459 """Apply a Shapely predicate to this shape and another one, transformed to this CRS.""" 

460 s = other.transformed_to(self.crs) 

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

462 

463 def union(self, others): 

464 if not others: 

465 return self 

466 

467 geoms = [self.geom] 

468 for s in others: 

469 s = s.transformed_to(self.crs) 

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

471 

472 geom = shapely.ops.unary_union(geoms) 

473 return Shape(geom, self.crs) 

474 

475 def intersection(self, *others): 

476 if not others: 

477 return self 

478 

479 geom = self.geom 

480 for s in others: 

481 s = s.transformed_to(self.crs) 

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

483 

484 return Shape(geom, self.crs) 

485 

486 def to_multi(self): 

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

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

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

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

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

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

493 return self 

494 

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

496 if new_type == self.type: 

497 return self 

498 if new_type == gws.GeometryType.geometry: 

499 return self 

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

501 return self.to_multi() 

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

503 return self.to_multi() 

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

505 return self.to_multi() 

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

507 

508 def to_2d(self): 

509 geom = shapely.force_2d(self.geom) 

510 if geom is self.geom: 

511 return self 

512 return Shape(geom, self.crs) 

513 

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

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

516 

517 if not tolerance and is_poly: 

518 return self 

519 

520 # we need a polygon even if tolerance = 0 

521 tolerance = tolerance or _MIN_TOLERANCE_RADIUS 

522 quad_segs = quad_segs or _TOLERANCE_QUAD_SEGS 

523 

524 geom = self.geom 

525 if self.crs.isGeographic: 

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

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

528 

529 if is_poly: 

530 cs = shapely.geometry.CAP_STYLE.flat 

531 js = shapely.geometry.JOIN_STYLE.mitre 

532 else: 

533 cs = shapely.geometry.CAP_STYLE.round 

534 js = shapely.geometry.JOIN_STYLE.round 

535 

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

537 if self.crs.isGeographic: 

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

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

540 

541 return Shape(geom, self.crs) 

542 

543 def transformed_to(self, crs): 

544 if crs == self.crs: 

545 return self 

546 tr = self.crs.transformer(crs) 

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

548 return Shape(dg, crs)