Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/gshhg.py: 96%

219 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Interface to the `GSHHG <https://www.soest.hawaii.edu/pwessel/gshhg/>`_ 

8(coastlines, rivers and borders) database. 

9 

10The Global Self-consistent Hierarchical High-resolution Geography Database 

11(GSHHG) is a collection of polygons representing land, lakes, rivers and 

12political borders. 

13 

14If the database is not already available, it will be downloaded 

15automatically on first use. 

16 

17For more information about GSHHG, see 

18http://www.soest.hawaii.edu/pwessel/gshhg/. 

19 

20.. rubric:: Citation 

21 

22If you use this dataset, please cite 

23 

24Wessel, P., and W. H. F. Smith, A Global Self-consistent, Hierarchical, 

25High-resolution Shoreline Database, J. Geophys. Res., 101, #B4, pp. 8741-8743, 

261996. 

27 

28''' 

29 

30import logging 

31import io 

32import struct 

33import time 

34import numpy as num 

35 

36from os import path 

37 

38from pyrocko import config, orthodrome 

39from .util import get_download_callback 

40 

41 

42logger = logging.getLogger('pyrocko.dataset.gshhg') 

43config = config.config() 

44 

45km = 1e3 

46micro_deg = 1e-6 

47 

48 

49def split_region_0_360(wesn): 

50 west, east, south, north = wesn 

51 if west < 0.: 

52 if east <= 0: 

53 return [(west+360., east+360., south, north)] 

54 else: 

55 return [(west+360., 360., south, north), 

56 (0., east, south, north)] 

57 else: 

58 return [wesn] 

59 

60 

61def is_valid_bounding_box(wesn): 

62 ''' 

63 Check if a given bounding box meets the GSHHG conventions. 

64 

65 :param wesn: bounding box as (west, east, south, north) in [deg] 

66 ''' 

67 

68 w, e, s, n = wesn 

69 

70 return ( 

71 w <= e 

72 and s <= n 

73 and -90.0 <= s <= 90. 

74 and -90. <= n <= 90. 

75 and -180. <= w < 360. 

76 and -180. <= e < 360.) 

77 

78 

79def is_valid_polygon(points): 

80 ''' 

81 Check if polygon points meet the GSHHG conventions. 

82 

83 :param points: Array of (lat, lon) pairs, shape (N, 2). 

84 ''' 

85 

86 lats = points[:, 0] 

87 lons = points[:, 1] 

88 

89 return ( 

90 num.all(-90. <= lats) 

91 and num.all(lats <= 90.) 

92 and num.all(-180. <= lons) 

93 and num.all(lons < 360.)) 

94 

95 

96def points_in_bounding_box(points, wesn, tolerance=0.1): 

97 ''' 

98 Check which points are contained in a given bounding box. 

99 

100 :param points: Array of (lat lon) pairs, shape (N, 2) [deg]. 

101 :param wesn: Region tuple (west, east, south, north) [deg] 

102 :param tolerance: increase the size of the test bounding box by 

103 *tolerance* [deg] on every side (Some GSHHG polygons have a too tight 

104 bounding box). 

105 

106 :returns: Bool array of shape (N,). 

107 ''' 

108 points_wrap = points.copy() 

109 points_wrap[:, 1] %= 360. 

110 

111 mask = num.zeros(points_wrap.shape[0], dtype=bool) 

112 for w, e, s, n in split_region_0_360(wesn): 

113 mask = num.logical_or( 

114 mask, 

115 num.logical_and( 

116 num.logical_and( 

117 w-tolerance <= points_wrap[:, 1], 

118 points_wrap[:, 1] <= e+tolerance), 

119 num.logical_and( 

120 s-tolerance <= points_wrap[:, 0], 

121 points_wrap[:, 0] <= n+tolerance))) 

122 

123 return mask 

124 

125 

126def point_in_bounding_box(point, wesn, tolerance=0.1): 

127 ''' 

128 Check whether point is contained in a given bounding box. 

129 

130 :param points: Array of (lat lon) pairs, shape (N, 2) [deg]. 

131 :param wesn: Region tuple (west, east, south, north) [deg] 

132 :param tolerance: increase the size of the test bounding box by 

133 *tolerance* [deg] on every side (Some GSHHG polygons have a too tight 

134 bounding box). 

135 

136 :rtype: bool 

137 ''' 

138 

139 lat, lon = point 

140 lon %= 360. 

141 for w, e, s, n in split_region_0_360(wesn): 

142 if (w-tolerance <= lon 

143 and lon <= e+tolerance 

144 and s-tolerance <= lat 

145 and lat <= n+tolerance): 

146 

147 return True 

148 

149 return False 

150 

151 

152def bounding_boxes_overlap(wesn1, wesn2): 

153 ''' 

154 Check whether two bounding boxes intersect. 

155 

156 :param wesn1: 

157 Region tuple ``(west, east, south, north)`` [deg] 

158 :type wesn1: 

159 :py:class:`tuple` of 4 :py:class:`float` 

160 

161 :param wesn2: 

162 Region tuple ``(west, east, south, north)`` [deg] 

163 :type wesn2: 

164 :py:class:`tuple` of 4 :py:class:`float` 

165 

166 :rtype: bool 

167 ''' 

168 for w1, e1, s1, n1 in split_region_0_360(wesn1): 

169 for w2, e2, s2, n2 in split_region_0_360(wesn2): 

170 if w2 <= e1 and w1 <= e2 and s2 <= n1 and s1 <= n2: 

171 return True 

172 

173 return False 

174 

175 

176def is_polygon_in_bounding_box(points, wesn, tolerance=0.1): 

177 return num.all(points_in_bounding_box(points, wesn, tolerance=tolerance)) 

178 

179 

180def bounding_box_covering_points(points): 

181 lats = points[:, 0] 

182 lat_min, lat_max = num.min(lats), num.max(lats) 

183 

184 lons = points[:, 1] 

185 lons = lons % 360. 

186 lon_min, lon_max = num.min(lons), num.max(lons) 

187 if lon_max - lon_min < 180.: 

188 return lon_min, lon_max, lat_min, lat_max 

189 

190 lons = (lons - 180.) % 360. - 180. 

191 lon_min, lon_max = num.min(lons), num.max(lons) 

192 if lon_max - lon_min < 180.: 

193 return lon_min, lon_max, lat_min, lat_max 

194 

195 return (-180., 180., lat_min, lat_max) 

196 

197 

198class Polygon(object): 

199 ''' 

200 Representation of a GSHHG polygon. 

201 ''' 

202 

203 RIVER_NOT_SET = 0 

204 

205 LEVELS = ['LAND', 'LAKE', 'ISLAND_IN_LAKE', 'POND_IN_ISLAND_IN_LAKE', 

206 'ANTARCTIC_ICE_FRONT', 'ANTARCTIC_GROUNDING_LINE'] 

207 

208 SOURCE = ['CIA_WDBII', 'WVS', 'AC'] 

209 

210 def __init__(self, gshhg_file, offset, *attr): 

211 ''' 

212 Initialise a GSHHG polygon 

213 

214 :param gshhg_file: GSHHG binary file 

215 :type gshhg_file: str 

216 :param offset: This polygon's offset in binary file 

217 :type offset: int 

218 :param attr: Polygon attributes 

219 ``(pid, npoints, _flag, west, east, south, north, 

220 area, area_full, container, ancestor)``. 

221 See :file:`gshhg.h` for details. 

222 :type attr: tuple 

223 ''' 

224 (self.pid, self.npoints, self._flag, 

225 self.west, self.east, self.south, self.north, 

226 self.area, self.area_full, self.container, self.ancestor) = attr 

227 

228 self.west *= micro_deg 

229 self.east *= micro_deg 

230 self.south *= micro_deg 

231 self.north *= micro_deg 

232 

233 self.level_no = (self._flag & 255) 

234 self.level = self.LEVELS[self.level_no - 1] 

235 self.version = (self._flag >> 8) & 255 

236 

237 cross = (self._flag >> 16) & 3 

238 self.greenwhich_crossed = True if cross == 1 or cross == 3 else False 

239 self.dateline_crossed = True if cross == 2 or cross == 3 else False 

240 

241 self.source = self.SOURCE[(self._flag >> 24) & 1] 

242 if self.level_no >= 5: 

243 self.source = self.SOURCE[2] 

244 

245 self.river = (self._flag >> 25) & 1 

246 

247 scale = 10.**(self._flag >> 26) 

248 self.area /= scale 

249 self.area_full /= scale 

250 

251 self._points = None 

252 self._file = gshhg_file 

253 self._offset = offset 

254 

255 @property 

256 def points(self): 

257 ''' 

258 Points of the polygon. 

259 

260 Array of (lat, lon) pairs, shape (N, 2). 

261 

262 :rtype: :class:`numpy.ndarray` 

263 ''' 

264 if self._points is None: 

265 with open(self._file) as db: 

266 db.seek(self._offset) 

267 self._points = num.fromfile( 

268 db, dtype='>i4', count=self.npoints*2)\ 

269 .astype(num.float32)\ 

270 .reshape(self.npoints, 2) 

271 

272 self._points = num.fliplr(self._points) 

273 if self.level_no in (2, 4): 

274 self._points = self._points[::-1, :] 

275 

276 self._points *= micro_deg 

277 return self._points 

278 

279 @property 

280 def lats(self): 

281 return self.points[:, 0] 

282 

283 @property 

284 def lons(self): 

285 return self.points[:, 1] 

286 

287 def _is_level(self, level): 

288 if self.level is self.LEVELS[level]: 

289 return True 

290 return False 

291 

292 def is_land(self): 

293 ''' 

294 Check if the polygon is land. 

295 

296 :rtype: bool 

297 ''' 

298 return self._is_level(0) 

299 

300 def is_lake(self): 

301 ''' 

302 Check if the polygon is a lake. 

303 

304 :rtype: bool 

305 ''' 

306 return self._is_level(1) 

307 

308 def is_island_in_lake(self): 

309 ''' 

310 Check if the polygon is an island in a lake. 

311 

312 :rtype: bool 

313 ''' 

314 return self._is_level(2) 

315 

316 def is_pond_in_island_in_lake(self): 

317 ''' 

318 Check if the polygon is pond on an island in a lake. 

319 

320 :rtype: bool 

321 ''' 

322 return self._is_level(3) 

323 

324 def is_antarctic_icefront(self): 

325 ''' 

326 Check if the polygon is antarctic icefront. 

327 

328 :rtype: bool 

329 ''' 

330 return self._is_level(4) 

331 

332 def is_antarctic_grounding_line(self): 

333 ''' 

334 Check if the polygon is antarctic grounding line. 

335 

336 :rtype: bool 

337 ''' 

338 return self._is_level(5) 

339 

340 def contains_point(self, point): 

341 ''' 

342 Check if point lies in polygon. 

343 

344 :param point: (lat, lon) [deg] 

345 :type point: tuple 

346 :rtype: bool 

347 

348 See :py:func:`pyrocko.orthodrome.contains_points`. 

349 ''' 

350 return bool( 

351 self.contains_points(num.asarray(point)[num.newaxis, :])[0]) 

352 

353 def contains_points(self, points): 

354 ''' 

355 Check if points lie in polygon. 

356 

357 :param points: Array of (lat lon) pairs, shape (N, 2) [deg]. 

358 :type points: :class:`numpy.ndarray` 

359 

360 See :py:func:`pyrocko.orthodrome.contains_points`. 

361 

362 :returns: Bool array of shape (N,) 

363 ''' 

364 mask = points_in_bounding_box(points, self.get_bounding_box()) 

365 if num.any(mask): 

366 mask[mask] = orthodrome.contains_points( 

367 self.points, points[mask, :]) 

368 

369 return mask 

370 

371 def get_bounding_box(self): 

372 return (self.west, self.east, self.south, self.north) 

373 

374 def __lt__(self, polygon): 

375 return self.level_no < polygon.level_no 

376 

377 def __str__(self): 

378 rstr = '''Polygon id: {p.pid} 

379------------------- 

380Points: {p.npoints} 

381Level: {p.level} 

382Area: {p.area} km**2 

383Area Full: {p.area_full} km**2 

384Extent: {p.west} W, {p.east} E, {p.south} S, {p.north} N 

385Source: {p.source} 

386Greenwhich crossed: {p.greenwhich_crossed} 

387Dateline crossed: {p.dateline_crossed} 

388 '''.format(p=self) 

389 return rstr 

390 

391 

392class GSHHG(object): 

393 ''' 

394 GSHHG database access. 

395 

396 This class provides methods to select relevant polygons (land, lakes, etc.) 

397 for given locations or regions. It also provides robust high-level 

398 functions to test if the Earth is dry or wet at given coordinates. 

399 ''' 

400 

401 gshhg_url = 'https://mirror.pyrocko.org/www.soest.hawaii.edu/pwessel/gshhg/gshhg-bin-2.3.7.zip' # noqa 

402 _header_struct = struct.Struct('>IIIiiiiIIii') 

403 

404 def __init__(self, gshhg_file): 

405 ''' Initialise the database from GSHHG binary. 

406 

407 :param gshhg_file: Path to file 

408 :type gshhg_file: str 

409 : 

410 ''' 

411 t0 = time.time() 

412 self._file = gshhg_file 

413 

414 self.polygons = [] 

415 self._read_database() 

416 logger.debug('Initialised GSHHG database from %s in [%.4f s]' 

417 % (gshhg_file, time.time()-t0)) 

418 

419 def _read_database(self): 

420 with open(self._file, mode='rb') as db: 

421 while db: 

422 buf = db.read(self._header_struct.size) 

423 if not buf: 

424 break 

425 header = self._header_struct.unpack_from(buf) 

426 p = Polygon( 

427 self._file, 

428 db.tell(), 

429 *header) 

430 self.polygons.append(p) 

431 

432 offset = 8 * header[1] 

433 db.seek(offset, io.SEEK_CUR) 

434 

435 @classmethod 

436 def _get_database(cls, filename): 

437 file = path.join(config.gshhg_dir, filename) 

438 if not path.exists(file): 

439 from pyrocko import util 

440 import zipfile 

441 

442 archive_path = path.join(config.gshhg_dir, 

443 path.basename(cls.gshhg_url)) 

444 util.download_file( 

445 cls.gshhg_url, archive_path, 

446 status_callback=get_download_callback( 

447 'Downloading GSHHG database...')) 

448 if not zipfile.is_zipfile(archive_path): 

449 raise util.DownloadError('GSHHG file is corrupted!') 

450 logger.info('Unzipping GSHHG database...') 

451 zipf = zipfile.ZipFile(archive_path) 

452 zipf.extractall(config.gshhg_dir) 

453 else: 

454 logger.debug('Using cached %s' % filename) 

455 return file 

456 

457 def get_polygons_at(self, lat, lon): 

458 ''' 

459 Get all polygons whose bounding boxes contain point. 

460 

461 :param lat: Latitude in [deg] 

462 :type lat: float 

463 :param lon: Longitude in [deg] 

464 :type lon: float 

465 :returns: List of :class:`~pyrocko.dataset.gshhg.Polygon` 

466 :rtype: list 

467 ''' 

468 rp = [] 

469 for p in self.polygons: 

470 if point_in_bounding_box((lat, lon), p.get_bounding_box()): 

471 rp.append(p) 

472 return rp 

473 

474 def get_polygons_within(self, west, east, south, north): 

475 ''' 

476 Get all polygons whose bounding boxes intersect with a bounding box. 

477 

478 :param west: Western boundary in decimal degree 

479 :type west: float 

480 :param east: Eastern boundary in decimal degree 

481 :type east: float 

482 :param north: Northern boundary in decimal degree 

483 :type north: float 

484 :param south: Southern boundary in decimal degree 

485 :type south: float 

486 :returns: List of :class:`~pyrocko.dataset.gshhg.Polygon` 

487 :rtype: list 

488 ''' 

489 

490 assert is_valid_bounding_box((west, east, south, north)) 

491 

492 rp = [] 

493 for p in self.polygons: 

494 if bounding_boxes_overlap( 

495 p.get_bounding_box(), (west, east, south, north)): 

496 

497 rp.append(p) 

498 return rp 

499 

500 def is_point_on_land(self, lat, lon): 

501 ''' 

502 Check whether a point is on land. 

503 

504 Lakes are considered not land. 

505 

506 :param lat: Latitude in [deg] 

507 :type lat: float 

508 :param lon: Longitude in [deg] 

509 :type lon: float 

510 

511 :rtype: bool 

512 ''' 

513 

514 relevant_polygons = self.get_polygons_at(lat, lon) 

515 relevant_polygons.sort() 

516 

517 land = False 

518 for p in relevant_polygons: 

519 if (p.is_land() or p.is_antarctic_grounding_line() or 

520 p.is_island_in_lake()): 

521 if p.contains_point((lat, lon)): 

522 land = True 

523 elif (p.is_lake() or p.is_antarctic_icefront() or 

524 p.is_pond_in_island_in_lake()): 

525 if p.contains_point((lat, lon)): 

526 land = False 

527 return land 

528 

529 def get_land_mask(self, points): 

530 ''' 

531 Check whether given points are on land. 

532 

533 Lakes are considered not land. 

534 

535 :param points: Array of (lat, lon) pairs, shape (N, 2). 

536 :type points: :class:`numpy.ndarray` 

537 :return: Boolean land mask 

538 :rtype: :class:`numpy.ndarray` of shape (N,) 

539 ''' 

540 

541 west, east, south, north = bounding_box_covering_points(points) 

542 

543 relevant_polygons = self.get_polygons_within(west, east, south, north) 

544 relevant_polygons.sort() 

545 

546 mask = num.zeros(points.shape[0], dtype=bool) 

547 for p in relevant_polygons: 

548 if (p.is_land() or p.is_antarctic_grounding_line() or 

549 p.is_island_in_lake()): 

550 land = p.contains_points(points) 

551 mask[land] = True 

552 elif p.is_lake() or p.is_pond_in_island_in_lake(): 

553 water = p.contains_points(points) 

554 mask[water] = False 

555 return mask 

556 

557 @classmethod 

558 def full(cls): 

559 ''' 

560 Return the full-resolution GSHHG database. 

561 ''' 

562 return cls(cls._get_database('gshhs_f.b')) 

563 

564 @classmethod 

565 def high(cls): 

566 ''' 

567 Return the high-resolution GSHHG database. 

568 ''' 

569 return cls(cls._get_database('gshhs_h.b')) 

570 

571 @classmethod 

572 def intermediate(cls): 

573 ''' 

574 Return the intermediate-resolution GSHHG database. 

575 ''' 

576 return cls(cls._get_database('gshhs_i.b')) 

577 

578 @classmethod 

579 def low(cls): 

580 ''' 

581 Return the low-resolution GSHHG database. 

582 ''' 

583 return cls(cls._get_database('gshhs_l.b')) 

584 

585 @classmethod 

586 def crude(cls): 

587 ''' 

588 Return the crude-resolution GSHHG database. 

589 ''' 

590 return cls(cls._get_database('gshhs_c.b'))