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-04 09:52 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Interface to the `GSHHG <https://www.soest.hawaii.edu/pwessel/gshhg/>`_
8(coastlines, rivers and borders) database.
10The Global Self-consistent Hierarchical High-resolution Geography Database
11(GSHHG) is a collection of polygons representing land, lakes, rivers and
12political borders.
14If the database is not already available, it will be downloaded
15automatically on first use.
17For more information about GSHHG, see
18http://www.soest.hawaii.edu/pwessel/gshhg/.
20.. rubric:: Citation
22If you use this dataset, please cite
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.
28'''
30import logging
31import io
32import struct
33import time
34import numpy as num
36from os import path
38from pyrocko import config, orthodrome
39from .util import get_download_callback
42logger = logging.getLogger('pyrocko.dataset.gshhg')
43config = config.config()
45km = 1e3
46micro_deg = 1e-6
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]
61def is_valid_bounding_box(wesn):
62 '''
63 Check if a given bounding box meets the GSHHG conventions.
65 :param wesn: bounding box as (west, east, south, north) in [deg]
66 '''
68 w, e, s, n = wesn
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.)
79def is_valid_polygon(points):
80 '''
81 Check if polygon points meet the GSHHG conventions.
83 :param points: Array of (lat, lon) pairs, shape (N, 2).
84 '''
86 lats = points[:, 0]
87 lons = points[:, 1]
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.))
96def points_in_bounding_box(points, wesn, tolerance=0.1):
97 '''
98 Check which points are contained in a given bounding box.
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).
106 :returns: Bool array of shape (N,).
107 '''
108 points_wrap = points.copy()
109 points_wrap[:, 1] %= 360.
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)))
123 return mask
126def point_in_bounding_box(point, wesn, tolerance=0.1):
127 '''
128 Check whether point is contained in a given bounding box.
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).
136 :rtype: bool
137 '''
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):
147 return True
149 return False
152def bounding_boxes_overlap(wesn1, wesn2):
153 '''
154 Check whether two bounding boxes intersect.
156 :param wesn1:
157 Region tuple ``(west, east, south, north)`` [deg]
158 :type wesn1:
159 :py:class:`tuple` of 4 :py:class:`float`
161 :param wesn2:
162 Region tuple ``(west, east, south, north)`` [deg]
163 :type wesn2:
164 :py:class:`tuple` of 4 :py:class:`float`
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
173 return False
176def is_polygon_in_bounding_box(points, wesn, tolerance=0.1):
177 return num.all(points_in_bounding_box(points, wesn, tolerance=tolerance))
180def bounding_box_covering_points(points):
181 lats = points[:, 0]
182 lat_min, lat_max = num.min(lats), num.max(lats)
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
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
195 return (-180., 180., lat_min, lat_max)
198class Polygon(object):
199 '''
200 Representation of a GSHHG polygon.
201 '''
203 RIVER_NOT_SET = 0
205 LEVELS = ['LAND', 'LAKE', 'ISLAND_IN_LAKE', 'POND_IN_ISLAND_IN_LAKE',
206 'ANTARCTIC_ICE_FRONT', 'ANTARCTIC_GROUNDING_LINE']
208 SOURCE = ['CIA_WDBII', 'WVS', 'AC']
210 def __init__(self, gshhg_file, offset, *attr):
211 '''
212 Initialise a GSHHG polygon
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
228 self.west *= micro_deg
229 self.east *= micro_deg
230 self.south *= micro_deg
231 self.north *= micro_deg
233 self.level_no = (self._flag & 255)
234 self.level = self.LEVELS[self.level_no - 1]
235 self.version = (self._flag >> 8) & 255
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
241 self.source = self.SOURCE[(self._flag >> 24) & 1]
242 if self.level_no >= 5:
243 self.source = self.SOURCE[2]
245 self.river = (self._flag >> 25) & 1
247 scale = 10.**(self._flag >> 26)
248 self.area /= scale
249 self.area_full /= scale
251 self._points = None
252 self._file = gshhg_file
253 self._offset = offset
255 @property
256 def points(self):
257 '''
258 Points of the polygon.
260 Array of (lat, lon) pairs, shape (N, 2).
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)
272 self._points = num.fliplr(self._points)
273 if self.level_no in (2, 4):
274 self._points = self._points[::-1, :]
276 self._points *= micro_deg
277 return self._points
279 @property
280 def lats(self):
281 return self.points[:, 0]
283 @property
284 def lons(self):
285 return self.points[:, 1]
287 def _is_level(self, level):
288 if self.level is self.LEVELS[level]:
289 return True
290 return False
292 def is_land(self):
293 '''
294 Check if the polygon is land.
296 :rtype: bool
297 '''
298 return self._is_level(0)
300 def is_lake(self):
301 '''
302 Check if the polygon is a lake.
304 :rtype: bool
305 '''
306 return self._is_level(1)
308 def is_island_in_lake(self):
309 '''
310 Check if the polygon is an island in a lake.
312 :rtype: bool
313 '''
314 return self._is_level(2)
316 def is_pond_in_island_in_lake(self):
317 '''
318 Check if the polygon is pond on an island in a lake.
320 :rtype: bool
321 '''
322 return self._is_level(3)
324 def is_antarctic_icefront(self):
325 '''
326 Check if the polygon is antarctic icefront.
328 :rtype: bool
329 '''
330 return self._is_level(4)
332 def is_antarctic_grounding_line(self):
333 '''
334 Check if the polygon is antarctic grounding line.
336 :rtype: bool
337 '''
338 return self._is_level(5)
340 def contains_point(self, point):
341 '''
342 Check if point lies in polygon.
344 :param point: (lat, lon) [deg]
345 :type point: tuple
346 :rtype: bool
348 See :py:func:`pyrocko.orthodrome.contains_points`.
349 '''
350 return bool(
351 self.contains_points(num.asarray(point)[num.newaxis, :])[0])
353 def contains_points(self, points):
354 '''
355 Check if points lie in polygon.
357 :param points: Array of (lat lon) pairs, shape (N, 2) [deg].
358 :type points: :class:`numpy.ndarray`
360 See :py:func:`pyrocko.orthodrome.contains_points`.
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, :])
369 return mask
371 def get_bounding_box(self):
372 return (self.west, self.east, self.south, self.north)
374 def __lt__(self, polygon):
375 return self.level_no < polygon.level_no
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
392class GSHHG(object):
393 '''
394 GSHHG database access.
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 '''
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')
404 def __init__(self, gshhg_file):
405 ''' Initialise the database from GSHHG binary.
407 :param gshhg_file: Path to file
408 :type gshhg_file: str
409 :
410 '''
411 t0 = time.time()
412 self._file = gshhg_file
414 self.polygons = []
415 self._read_database()
416 logger.debug('Initialised GSHHG database from %s in [%.4f s]'
417 % (gshhg_file, time.time()-t0))
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)
432 offset = 8 * header[1]
433 db.seek(offset, io.SEEK_CUR)
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
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
457 def get_polygons_at(self, lat, lon):
458 '''
459 Get all polygons whose bounding boxes contain point.
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
474 def get_polygons_within(self, west, east, south, north):
475 '''
476 Get all polygons whose bounding boxes intersect with a bounding box.
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 '''
490 assert is_valid_bounding_box((west, east, south, north))
492 rp = []
493 for p in self.polygons:
494 if bounding_boxes_overlap(
495 p.get_bounding_box(), (west, east, south, north)):
497 rp.append(p)
498 return rp
500 def is_point_on_land(self, lat, lon):
501 '''
502 Check whether a point is on land.
504 Lakes are considered not land.
506 :param lat: Latitude in [deg]
507 :type lat: float
508 :param lon: Longitude in [deg]
509 :type lon: float
511 :rtype: bool
512 '''
514 relevant_polygons = self.get_polygons_at(lat, lon)
515 relevant_polygons.sort()
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
529 def get_land_mask(self, points):
530 '''
531 Check whether given points are on land.
533 Lakes are considered not land.
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 '''
541 west, east, south, north = bounding_box_covering_points(points)
543 relevant_polygons = self.get_polygons_within(west, east, south, north)
544 relevant_polygons.sort()
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
557 @classmethod
558 def full(cls):
559 '''
560 Return the full-resolution GSHHG database.
561 '''
562 return cls(cls._get_database('gshhs_f.b'))
564 @classmethod
565 def high(cls):
566 '''
567 Return the high-resolution GSHHG database.
568 '''
569 return cls(cls._get_database('gshhs_h.b'))
571 @classmethod
572 def intermediate(cls):
573 '''
574 Return the intermediate-resolution GSHHG database.
575 '''
576 return cls(cls._get_database('gshhs_i.b'))
578 @classmethod
579 def low(cls):
580 '''
581 Return the low-resolution GSHHG database.
582 '''
583 return cls(cls._get_database('gshhs_l.b'))
585 @classmethod
586 def crude(cls):
587 '''
588 Return the crude-resolution GSHHG database.
589 '''
590 return cls(cls._get_database('gshhs_c.b'))