Geographical datasets

Pyrocko offers access to commonly used geographical datasets, such as


The pyrocko.dataset.topo subpackage offers quick access to some popular global topography datasets.

The following example draws gridded topography for Mount Vesuvius from SRTMGL3 ([1], [2], [3], [4]) in local cartesian coordinates.

import numpy as num
from pyrocko import plot, orthodrome as od
from pyrocko.dataset import topo

lon_min, lon_max, lat_min, lat_max = 14.34, 14.50, 40.77, 40.87
dem_name = 'SRTMGL3'

# extract gridded topography (possibly downloading first)
tile = topo.get(dem_name, (lon_min, lon_max, lat_min, lat_max))

# geographic to local cartesian coordinates
lons = tile.x()
lats = tile.y()
lons2 = num.tile(lons, lats.size)
lats2 = num.repeat(lats, lons.size)
norths, easts = od.latlon_to_ne_numpy(lats[0], lons[0], lats2, lons2)
norths = norths.reshape((lats.size, lons.size))
easts = easts.reshape((lats.size, lons.size))

# plot it
plt = plot.mpl_init(fontsize=10.)
fig = plt.figure(figsize=plot.mpl_papersize('a5', 'landscape'))
axes = fig.add_subplot(1, 1, 1, aspect=1.0)
cbar = axes.pcolormesh(easts, norths,,
                       cmap='gray', shading='gouraud')
fig.colorbar(cbar).set_label('Altitude [m]')
axes.set_xlim(easts.min(), easts.max())
axes.set_ylim(norths.min(), norths.max())
axes.set_xlabel('Easting [m]')
axes.set_ylabel('Northing [m]')


GSHHG coastal database

The GSHHG database is a high-resolution geography data set. We implement functions to extract coordinates of landmasks.

Classes covered in this example:
import numpy as num
from pyrocko.dataset.gshhg import GSHHG
from matplotlib import pyplot as plt

gshhg = GSHHG.intermediate()
# gshhg = GSHHG.full()
# gshhg = GSHHG.low()

gshhg.is_point_on_land(lat=32.1, lon=44.2)

# Create a landmask for a regular grid
lats = num.linspace(30., 50., 500)
lons = num.linspace(2., 10., 500)

lon_grid, lat_grid = num.meshgrid(lons, lats)
coordinates = num.array([lat_grid.ravel(), lon_grid.ravel()]).T

land_mask = gshhg.get_land_mask(coordinates).reshape(*lat_grid.shape)

plt.pcolormesh(lons, lats, land_mask, cmap='Greys', shading='nearest')



Tectonic plates and boundaries (PB2003)

An updated digital model of plate boundaries [6] offers a global set of present plate boundaries on the Earth. This database used in pyrocko.plot.automap.

Classes covered in this example:
from pyrocko.dataset.tectonics import PeterBird2003

poi_lat = 12.4
poi_lon = 133.5

PB = PeterBird2003()
plates = PB.get_plates()

for plate in plates:
    if plate.contains_point((poi_lat, poi_lon)):

>>> PS
>>> Philippine Sea Plate



Global strain rate (GSRM1)

Access to the global strain rate data set from Kreemer et al. (2003) [7].

Classes to be covered in this example:


To be documented by example!