Geographical datasets

Pyrocko offers access to commonly used geographical datasets, such as

Topography

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, tile.data,
                       cmap='gray', shading='gouraud')
fig.colorbar(cbar).set_label('Altitude [m]')
axes.set_title(dem_name)
axes.set_xlim(easts.min(), easts.max())
axes.set_ylim(norths.min(), norths.max())
axes.set_xlabel('Easting [m]')
axes.set_ylabel('Northing [m]')
fig.savefig('topo_example.png')
../../_images/topo_example.png

Footnotes

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')
plt.show()

Download gshhg_example.py

Footnotes

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.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)):
        print(plate.name)
        print(PB.full_name(plate.name))

'''
>>> PS
>>> Philippine Sea Plate
'''

Download tectonics_example.py

Footnotes

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:

Warning

To be documented by example!

Footnotes