1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Interface to the GSHHG (coastlines, rivers and borders) database. 

8 

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

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

11political borders. 

12 

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

14automatically on first use. 

15 

16For more information about GSHHG, see 

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

18 

19.. note:: 

20 

21 **If you use this dataset, please cite:** 

22 

23 Wessel, P., and W. H. F. 

24 Smith, A Global Self-consistent, Hierarchical, High-resolution 

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

26''' 

27 

28from __future__ import absolute_import 

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, wesn2: Region tuples (west, east, south, north) [deg] 

157 

158 :rtype: bool 

159 ''' 

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

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

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

163 return True 

164 

165 return False 

166 

167 

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

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

170 

171 

172def bounding_box_covering_points(points): 

173 lats = points[:, 0] 

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

175 

176 lons = points[:, 1] 

177 lons = lons % 360. 

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

179 if lon_max - lon_min < 180.: 

180 return lon_min, lon_max, lat_min, lat_max 

181 

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

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

184 if lon_max - lon_min < 180.: 

185 return lon_min, lon_max, lat_min, lat_max 

186 

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

188 

189 

190class Polygon(object): 

191 ''' 

192 Representation of a GSHHG polygon. 

193 ''' 

194 

195 RIVER_NOT_SET = 0 

196 

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

198 'ANTARCTIC_ICE_FRONT', 'ANTARCTIC_GROUNDING_LINE'] 

199 

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

201 

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

203 ''' 

204 Initialise a GSHHG polygon 

205 

206 :param gshhg_file: GSHHG binary file 

207 :type gshhg_file: str 

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

209 :type offset: int 

210 :param attr: Polygon attributes 

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

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

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

214 :type attr: tuple 

215 ''' 

216 (self.pid, self.npoints, self._flag, 

217 self.west, self.east, self.south, self.north, 

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

219 

220 self.west *= micro_deg 

221 self.east *= micro_deg 

222 self.south *= micro_deg 

223 self.north *= micro_deg 

224 

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

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

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

228 

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

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

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

232 

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

234 if self.level_no >= 5: 

235 self.source = self.SOURCE[2] 

236 

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

238 

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

240 self.area /= scale 

241 self.area_full /= scale 

242 

243 self._points = None 

244 self._file = gshhg_file 

245 self._offset = offset 

246 

247 @property 

248 def points(self): 

249 ''' 

250 Points of the polygon. 

251 

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

253 

254 :rtype: :class:`numpy.ndarray` 

255 ''' 

256 if self._points is None: 

257 with open(self._file) as db: 

258 db.seek(self._offset) 

259 self._points = num.fromfile( 

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

261 .astype(num.float32)\ 

262 .reshape(self.npoints, 2) 

263 

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

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

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

267 

268 self._points *= micro_deg 

269 return self._points 

270 

271 @property 

272 def lats(self): 

273 return self.points[:, 0] 

274 

275 @property 

276 def lons(self): 

277 return self.points[:, 1] 

278 

279 def _is_level(self, level): 

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

281 return True 

282 return False 

283 

284 def is_land(self): 

285 ''' 

286 Check if the polygon is land. 

287 

288 :rtype: bool 

289 ''' 

290 return self._is_level(0) 

291 

292 def is_lake(self): 

293 ''' 

294 Check if the polygon is a lake. 

295 

296 :rtype: bool 

297 ''' 

298 return self._is_level(1) 

299 

300 def is_island_in_lake(self): 

301 ''' 

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

303 

304 :rtype: bool 

305 ''' 

306 return self._is_level(2) 

307 

308 def is_pond_in_island_in_lake(self): 

309 ''' 

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

311 

312 :rtype: bool 

313 ''' 

314 return self._is_level(3) 

315 

316 def is_antarctic_icefront(self): 

317 ''' 

318 Check if the polygon is antarctic icefront. 

319 

320 :rtype: bool 

321 ''' 

322 return self._is_level(4) 

323 

324 def is_antarctic_grounding_line(self): 

325 ''' 

326 Check if the polygon is antarctic grounding line. 

327 

328 :rtype: bool 

329 ''' 

330 return self._is_level(5) 

331 

332 def contains_point(self, point): 

333 ''' 

334 Check if point lies in polygon. 

335 

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

337 :type point: tuple 

338 :rtype: bool 

339 

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

341 ''' 

342 return bool( 

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

344 

345 def contains_points(self, points): 

346 ''' 

347 Check if points lie in polygon. 

348 

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

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

351 

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

353 

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

355 ''' 

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

357 if num.any(mask): 

358 mask[mask] = orthodrome.contains_points( 

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

360 

361 return mask 

362 

363 def get_bounding_box(self): 

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

365 

366 def __lt__(self, polygon): 

367 return self.level_no < polygon.level_no 

368 

369 def __str__(self): 

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

371------------------- 

372Points: {p.npoints} 

373Level: {p.level} 

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

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

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

377Source: {p.source} 

378Greenwhich crossed: {p.greenwhich_crossed} 

379Dateline crossed: {p.dateline_crossed} 

380 '''.format(p=self) 

381 return rstr 

382 

383 

384class GSHHG(object): 

385 ''' 

386 GSHHG database access. 

387 

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

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

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

391 ''' 

392 

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

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

395 

396 def __init__(self, gshhg_file): 

397 ''' Initialise the database from GSHHG binary. 

398 

399 :param gshhg_file: Path to file 

400 :type gshhg_file: str 

401 : 

402 ''' 

403 t0 = time.time() 

404 self._file = gshhg_file 

405 

406 self.polygons = [] 

407 self._read_database() 

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

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

410 

411 def _read_database(self): 

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

413 while db: 

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

415 if not buf: 

416 break 

417 header = self._header_struct.unpack_from(buf) 

418 p = Polygon( 

419 self._file, 

420 db.tell(), 

421 *header) 

422 self.polygons.append(p) 

423 

424 offset = 8 * header[1] 

425 db.seek(offset, io.SEEK_CUR) 

426 

427 @classmethod 

428 def _get_database(cls, filename): 

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

430 if not path.exists(file): 

431 from pyrocko import util 

432 import zipfile 

433 

434 archive_path = path.join(config.gshhg_dir, 

435 path.basename(cls.gshhg_url)) 

436 util.download_file( 

437 cls.gshhg_url, archive_path, 

438 status_callback=get_download_callback( 

439 'Downloading GSHHG database...')) 

440 if not zipfile.is_zipfile(archive_path): 

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

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

443 zipf = zipfile.ZipFile(archive_path) 

444 zipf.extractall(config.gshhg_dir) 

445 else: 

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

447 return file 

448 

449 def get_polygons_at(self, lat, lon): 

450 ''' 

451 Get all polygons whose bounding boxes contain point. 

452 

453 :param lat: Latitude in [deg] 

454 :type lat: float 

455 :param lon: Longitude in [deg] 

456 :type lon: float 

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

458 :rtype: list 

459 ''' 

460 rp = [] 

461 for p in self.polygons: 

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

463 rp.append(p) 

464 return rp 

465 

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

467 ''' 

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

469 

470 :param west: Western boundary in decimal degree 

471 :type west: float 

472 :param east: Eastern boundary in decimal degree 

473 :type east: float 

474 :param north: Northern boundary in decimal degree 

475 :type north: float 

476 :param south: Southern boundary in decimal degree 

477 :type south: float 

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

479 :rtype: list 

480 ''' 

481 

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

483 

484 rp = [] 

485 for p in self.polygons: 

486 if bounding_boxes_overlap( 

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

488 

489 rp.append(p) 

490 return rp 

491 

492 def is_point_on_land(self, lat, lon): 

493 ''' 

494 Check whether a point is on land. 

495 

496 Lakes are considered not land. 

497 

498 :param lat: Latitude in [deg] 

499 :type lat: float 

500 :param lon: Longitude in [deg] 

501 :type lon: float 

502 

503 :rtype: bool 

504 ''' 

505 

506 relevant_polygons = self.get_polygons_at(lat, lon) 

507 relevant_polygons.sort() 

508 

509 land = False 

510 for p in relevant_polygons: 

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

512 p.is_island_in_lake()): 

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

514 land = True 

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

516 p.is_pond_in_island_in_lake()): 

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

518 land = False 

519 return land 

520 

521 def get_land_mask(self, points): 

522 ''' 

523 Check whether given points are on land. 

524 

525 Lakes are considered not land. 

526 

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

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

529 :return: Boolean land mask 

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

531 ''' 

532 

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

534 

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

536 relevant_polygons.sort() 

537 

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

539 for p in relevant_polygons: 

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

541 p.is_island_in_lake()): 

542 land = p.contains_points(points) 

543 mask[land] = True 

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

545 water = p.contains_points(points) 

546 mask[water] = False 

547 return mask 

548 

549 @classmethod 

550 def full(cls): 

551 ''' 

552 Return the full-resolution GSHHG database. 

553 ''' 

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

555 

556 @classmethod 

557 def high(cls): 

558 ''' 

559 Return the high-resolution GSHHG database. 

560 ''' 

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

562 

563 @classmethod 

564 def intermediate(cls): 

565 ''' 

566 Return the intermediate-resolution GSHHG database. 

567 ''' 

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

569 

570 @classmethod 

571 def low(cls): 

572 ''' 

573 Return the low-resolution GSHHG database. 

574 ''' 

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

576 

577 @classmethod 

578 def crude(cls): 

579 ''' 

580 Return the crude-resolution GSHHG database. 

581 ''' 

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