|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
Interface to the GSHHG (coastlines, rivers and borders) database.
The Global Self-consistent Hierarchical High-resolution Geography Database (GSHHG) is a collection of polygons representing land, lakes, rivers and political borders.
If the database is not already available, it will be downloaded automatically on first use.
For more information about GSHHG, see http://www.soest.hawaii.edu/pwessel/gshhg/.
.. note::
**If you use this dataset, please cite:**
Wessel, P., and W. H. F. Smith, A Global Self-consistent, Hierarchical, High-resolution Shoreline Database, J. Geophys. Res., 101, #B4, pp. 8741-8743, 1996. '''
return [(west+360., east+360., south, north)] else: (0., east, south, north)] else:
''' Check if a given bounding box meets the GSHHG conventions.
:param wesn: bounding box as (west, east, south, north) in [deg] '''
w, e, s, n = wesn
return ( w <= e and s <= n and -90.0 <= s <= 90. and -90. <= n <= 90. and -180. <= w < 360. and -180. <= e < 360.)
''' Check if polygon points meet the GSHHG conventions.
:param points: Array of (lat, lon) pairs, shape (N, 2). '''
lats = points[:, 0] lons = points[:, 1]
return ( num.all(-90. <= lats) and num.all(lats <= 90.) and num.all(-180. <= lons) and num.all(lons < 360.))
''' Check which points are contained in a given bounding box.
:param points: Array of (lat lon) pairs, shape (N, 2) [deg]. :param wesn: Region tuple (west, east, south, north) [deg] :param tolerance: increase the size of the test bounding box by *tolerance* [deg] on every side (Some GSHHG polygons have a too tight bounding box).
:returns: Bool array of shape (N,). '''
mask, num.logical_and( num.logical_and( w-tolerance <= points_wrap[:, 1], points_wrap[:, 1] <= e+tolerance), num.logical_and( s-tolerance <= points_wrap[:, 0], points_wrap[:, 0] <= n+tolerance)))
''' Check whether point is contained in a given bounding box.
:param points: Array of (lat lon) pairs, shape (N, 2) [deg]. :param wesn: Region tuple (west, east, south, north) [deg] :param tolerance: increase the size of the test bounding box by *tolerance* [deg] on every side (Some GSHHG polygons have a too tight bounding box).
:rtype: bool '''
and lon <= e+tolerance and s-tolerance <= lat and lat <= n+tolerance):
''' Check whether two bounding boxes intersect.
:param wesn1, wesn2: Region tuples (west, east, south, north) [deg]
:rtype: bool ''' for w1, e1, s1, n1 in split_region_0_360(wesn1): for w2, e2, s2, n2 in split_region_0_360(wesn2): if w2 <= e1 and w1 <= e2 and s2 <= n1 and s1 <= n2: return True
return False
return num.all(points_in_bounding_box(points, wesn, tolerance=tolerance))
lats = points[:, 0] lat_min, lat_max = num.min(lats), num.max(lats)
lons = points[:, 1] lons = lons % 360. lon_min, lon_max = num.min(lons), num.max(lons) if lon_max - lon_min < 180.: return lon_min, lon_max, lat_min, lat_max
lons = (lons - 180.) % 360. - 180. lon_min, lon_max = num.min(lons), num.max(lons) if lon_max - lon_min < 180.: return lon_min, lon_max, lat_min, lat_max
return (-180., 180., lat_min, lat_max)
''' Representation of a GSHHG polygon. '''
'ANTARCTIC_ICE_FRONT', 'ANTARCTIC_GROUNDING_LINE']
''' Initialise a GSHHG polygon
:param gshhg_file: GSHHG binary file :type gshhg_file: str :param offset: This polygon's offset in binary file :type offset: int :param attr: Polygon attributes ``(pid, npoints, _flag, west, east, south, north, area, area_full, container, ancestor)``. See :file:`gshhg.h` for details. :type attr: tuple ''' self.west, self.east, self.south, self.north, self.area, self.area_full, self.container, self.ancestor) = attr
def points(self): ''' Points of the polygon.
Array of (lat, lon) pairs, shape (N, 2).
:rtype: :class:`numpy.ndarray` ''' db, dtype='>i4', count=self.npoints*2)\ .astype(num.float32)\ .reshape(self.npoints, 2)
self._points = self._points[::-1, :]
def lats(self): return self.points[:, 0]
def lons(self): return self.points[:, 1]
return False
''' Check if the polygon is land.
:rtype: bool '''
''' Check if the polygon is a lake.
:rtype: bool ''' return self._is_level(1)
''' Check if the polygon is an island in a lake.
:rtype: bool ''' return self._is_level(2)
''' Check if the polygon is pond on an island in a lake.
:rtype: bool ''' return self._is_level(3)
''' Check if the polygon is antarctic icefront.
:rtype: bool ''' return self._is_level(4)
''' Check if the polygon is antarctic grounding line.
:rtype: bool ''' return self._is_level(5)
''' Check if point lies in polygon.
:param point: (lat, lon) [deg] :type point: tuple :rtype: bool
See :py:func:`pyrocko.orthodrome.contains_points`. ''' self.contains_points(num.asarray(point)[num.newaxis, :])[0])
''' Check if points lie in polygon.
:param points: Array of (lat lon) pairs, shape (N, 2) [deg]. :type points: :class:`numpy.ndarray`
See :py:func:`pyrocko.orthodrome.contains_points`.
:returns: Bool array of shape (N,) ''' self.points, points[mask, :])
return self.level_no < polygon.level_no
def __str__(self): rstr = '''Polygon id: {p.pid} ------------------- Points: {p.npoints} Level: {p.level} Area: {p.area} km**2 Area Full: {p.area_full} km**2 Extent: {p.west} W, {p.east} E, {p.south} S, {p.north} N Source: {p.source} Greenwhich crossed: {p.greenwhich_crossed} Dateline crossed: {p.dateline_crossed} '''.format(p=self) return rstr
''' GSHHG database access.
This class provides methods to select relevant polygons (land, lakes, etc.) for given locations or regions. It also provides robust high-level functions to test if the Earth is dry or wet at given coordinates. '''
''' Initialise the database from GSHHG binary.
:param gshhg_file: Path to file :type gshhg_file: str : '''
% (gshhg_file, time.time()-t0))
self._file, db.tell(), *header)
def _get_database(cls, filename):
path.basename(cls.gshhg_url)) cls.gshhg_url, archive_path, status_callback=get_download_callback( 'Downloading GSHHG database...')) raise util.DownloadError('GSHHG file is corrupted!') else: logger.debug('Using cached %s' % filename)
''' Get all polygons whose bounding boxes contain point.
:param lat: Latitude in [deg] :type lat: float :param lon: Longitude in [deg] :type lon: float :returns: List of :class:`~pyrocko.dataset.gshhg.Polygon` :rtype: list '''
''' Get all polygons whose bounding boxes intersect with a bounding box.
:param west: Western boundary in decimal degree :type west: float :param east: Eastern boundary in decimal degree :type east: float :param north: Northern boundary in decimal degree :type north: float :param south: Southern boundary in decimal degree :type south: float :returns: List of :class:`~pyrocko.dataset.gshhg.Polygon` :rtype: list '''
assert is_valid_bounding_box((west, east, south, north))
rp = [] for p in self.polygons: if bounding_boxes_overlap( p.get_bounding_box(), (west, east, south, north)):
rp.append(p) return rp
''' Check whether a point is on land.
Lakes are considered not land.
:param lat: Latitude in [deg] :type lat: float :param lon: Longitude in [deg] :type lon: float
:rtype: bool '''
p.is_island_in_lake()): elif (p.is_lake() or p.is_antarctic_icefront() or p.is_pond_in_island_in_lake()): if p.contains_point((lat, lon)): land = False
''' Check whether given points are on land.
Lakes are considered not land.
:param points: Array of (lat, lon) pairs, shape (N, 2). :type points: :class:`numpy.ndarray` :return: Boolean land mask :rtype: :class:`numpy.ndarray` of shape (N,) '''
west, east, south, north = bounding_box_covering_points(points)
relevant_polygons = self.get_polygons_within(west, east, south, north) relevant_polygons.sort()
mask = num.zeros(points.shape[0], dtype=num.bool) for p in relevant_polygons: if (p.is_land() or p.is_antarctic_grounding_line() or p.is_island_in_lake()): land = p.contains_points(points) mask[land] = True elif p.is_lake() or p.is_pond_in_island_in_lake(): water = p.contains_points(points) mask[water] = False return mask
def full(cls): ''' Return the full-resolution GSHHG database. ''' return cls(cls._get_database('gshhs_f.b'))
def high(cls): ''' Return the high-resolution GSHHG database. ''' return cls(cls._get_database('gshhs_h.b'))
def intermediate(cls): ''' Return the intermediate-resolution GSHHG database. '''
def low(cls): ''' Return the low-resolution GSHHG database. ''' return cls(cls._get_database('gshhs_l.b'))
def crude(cls): ''' Return the crude-resolution GSHHG database. ''' return cls(cls._get_database('gshhs_c.b')) |