# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
'''Interface to use CRUST2.0 model by Laske, Masters and Reif.
All functions defined in this module return SI units (m, m/s, kg/m^3).
.. note::
Please refer to the REM web site if you use this model:
http://igppweb.ucsd.edu/~gabi/rem.html
or
Bassin, C., Laske, G. and Masters, G., The Current Limits of Resolution for
Surface Wave Tomography in North America, EOS Trans AGU, 81, F897, 2000. A
description of CRUST 5.1 can be found in: Mooney, Laske and Masters, Crust
5.1: a global crustal model at 5x5 degrees, JGR, 103, 727-747, 1998.
Usage
-----
::
>>> from pyrocko import crust2x2
>>> p = crust2x2.get_profile(10., 20.)
>>> print p
type, name: G2, Archean 0.5 km seds.
elevation: 529
crustal thickness: 38500
average vp, vs, rho: 6460.7 3665.1 2867.5
mantle ave. vp, vs, rho: 8200 4700 3400
0 3810 1940 920 ice
0 1500 0 1020 water
500 2500 1200 2100 soft sed.
0 4000 2100 2400 hard sed.
12500 6200 3600 2800 upper crust
13000 6400 3600 2850 middle crust
13000 6800 3800 2950 lower crust
>>> print p.get_weeded()
[[ 0. 500. 500. 13000. 13000. 26000. 26000. 39000. 39000.]
[ 2500. 2500. 6200. 6200. 6400. 6400. 6800. 6800. 8200.]
[ 1200. 1200. 3600. 3600. 3600. 3600. 3800. 3800. 4700.]
[ 2100. 2100. 2800. 2800. 2850. 2850. 2950. 2950. 3400.]]
Constants
---------
============== ==============
Layer id Layer name
============== ==============
LICE ice
LWATER water
LSOFTSED soft sediments
LHARDSED hard sediments
LUPPERCRUST upper crust
LMIDDLECRUST middle crust
LLOWERCRUST lower crust
LBELOWCRUST below crust
============== ==============
Contents
--------
'''
from __future__ import absolute_import, division
import os
import copy
import math
from io import StringIO
import numpy as num
LICE, LWATER, LSOFTSED, LHARDSED, LUPPERCRUST, LMIDDLECRUST, \
LLOWERCRUST, LBELOWCRUST = list(range(8))
[docs]class Crust2Profile(object):
'''Representation of a CRUST2.0 key profile.'''
layer_names = (
'ice', 'water', 'soft sed.', 'hard sed.', 'upper crust',
'middle crust', 'lower crust', 'mantle')
def __init__(self, ident, name, vp, vs, rho, thickness, elevation):
self._ident = ident
self._name = name
self._vp = vp
self._vs = vs
self._rho = rho
self._thickness = thickness
self._crustal_thickness = None
self._elevation = elevation
[docs] def get_weeded(self, include_waterlayer=False):
'''Get layers used in the profile.
:param include_waterlayer: include water layer if ``True``. Default is
``False``
:returns: NumPy array with rows ``depth``, ``vp``, ``vs``, ``density``
'''
depth = 0.
layers = []
for ilayer, thickness, vp, vs, rho in zip(
range(8),
self._thickness,
self._vp[:-1],
self._vs[:-1],
self._rho[:-1]):
if thickness == 0.0:
continue
if not include_waterlayer and ilayer == LWATER:
continue
layers.append([depth, vp, vs, rho])
layers.append([depth+thickness, vp, vs, rho])
depth += thickness
layers.append([
depth,
self._vp[LBELOWCRUST],
self._vs[LBELOWCRUST],
self._rho[LBELOWCRUST]])
return num.array(layers).T
[docs] def get_layer(self, ilayer):
'''Get parameters for a layer.
:param ilayer: id of layer
:returns: thickness, vp, vs, density
'''
if ilayer == LBELOWCRUST:
thickness = num.Inf
else:
thickness = self._thickness[ilayer]
return thickness, self._vp[ilayer], self._vs[ilayer], self._rho[ilayer]
[docs] def set_elevation(self, elevation):
self._elevation = elevation
[docs] def set_layer_thickness(self, ilayer, thickness):
self._thickness[ilayer] = thickness
[docs] def elevation(self):
return self._elevation
def __str__(self):
vvp, vvs, vrho, vthi = self.averages()
return '''type, name: %s, %s
elevation: %15.5g
crustal thickness: %15.5g
average vp, vs, rho: %15.5g %15.5g %15.5g
mantle ave. vp, vs, rho: %15.5g %15.5g %15.5g
%s''' % (self._ident, self._name, self._elevation, vthi, vvp, vvs, vrho,
self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST],
'\n'.join([
'%15.5g %15.5g %15.5g %15.5g %s' % x
for x in zip(
self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1],
Crust2Profile.layer_names)]))
[docs] def crustal_thickness(self):
'''Get total crustal thickness
Takes into account ice layer.
Does not take into account water layer.
'''
return num.sum(self._thickness[3:]) + self._thickness[LICE]
[docs] def averages(self):
'''Get crustal averages for vp, vs and density and total crustal thickness,
Takes into account ice layer.
Does not take into account water layer.
'''
vthi = self.crustal_thickness()
vvp = num.sum(self._thickness[3:] / self._vp[3:-1]) + \
self._thickness[LICE] / self._vp[LICE]
vvs = num.sum(self._thickness[3:] / self._vs[3:-1]) + \
self._thickness[LICE] / self._vs[LICE]
vrho = num.sum(self._thickness[3:] * self._rho[3:-1]) + \
self._thickness[LICE] * self._rho[LICE]
vvp = vthi / vvp
vvs = vthi / vvs
vrho = vrho / vthi
return vvp, vvs, vrho, vthi
def _sa2arr(sa):
return num.array([float(x) for x in sa], dtype=num.float)
def _wrap(x, mi, ma):
if mi <= x and x <= ma:
return x
return x - math.floor((x-mi)/(ma-mi)) * (ma-mi)
def _clip(x, mi, ma):
return min(max(mi, x), ma)
[docs]class Crust2(object):
'''Access CRUST2.0 model.
:param directory: Directory with the data files which contain the
CRUST2.0 model data. If this is set to ``None``, builtin CRUST2.0
files are used.
'''
fn_keys = 'CNtype2_key.txt'
fn_elevation = 'CNelevatio2.txt'
fn_map = 'CNtype2.txt'
nlo = 180
nla = 90
_instance = None
def __init__(self, directory=None):
self.profile_keys = []
self._directory = directory
self._typemap = None
self._load_crustal_model()
[docs] def get_profile(self, *args, **kwargs):
'''
Get crustal profile at a specific location or raw profile for given
key.
Get profile for location ``(lat, lon)``, or raw profile for given
string key.
:rtype: instance of :py:class:`Crust2Profile`
'''
lat = kwargs.pop('lat', None)
lon = kwargs.pop('lon', None)
if len(args) == 2:
lat, lon = args
if lat is not None and lon is not None:
return self._typemap[self._indices(float(lat), float(lon))]
else:
return self._raw_profiles[args[0]]
def _indices(self, lat, lon):
lat = _clip(lat, -90., 90.)
lon = _wrap(lon, -180., 180.)
dlo = 360./Crust2.nlo
dla = 180./Crust2.nla
cola = 90.-lat
ilat = _clip(int(cola/dla), 0, Crust2.nla-1)
ilon = int((lon+180.)/dlo) % Crust2.nlo
return ilat, ilon
def _load_crustal_model(self):
if self._directory is not None:
path_keys = os.path.join(self._directory, Crust2.fn_keys)
f = open(path_keys, 'r')
else:
from .crust2x2_data import decode, type2_key, type2, elevation
f = StringIO(decode(type2_key))
# skip header
for i in range(5):
f.readline()
profiles = {}
while True:
line = f.readline()
if not line:
break
ident, name = line.split(None, 1)
line = f.readline()
vp = _sa2arr(line.split()) * 1000.
line = f.readline()
vs = _sa2arr(line.split()) * 1000.
line = f.readline()
rho = _sa2arr(line.split()) * 1000.
line = f.readline()
toks = line.split()
thickness = _sa2arr(toks[:-2]) * 1000.
assert ident not in profiles
profiles[ident] = Crust2Profile(
ident.strip(), name.strip(), vp, vs, rho, thickness, 0.0)
f.close()
self._raw_profiles = profiles
self.profile_keys = sorted(profiles.keys())
if self._directory is not None:
path_map = os.path.join(self._directory, Crust2.fn_map)
f = open(path_map, 'r')
else:
f = StringIO(decode(type2))
f.readline() # header
amap = {}
for ila, line in enumerate(f):
keys = line.split()[1:]
for ilo, key in enumerate(keys):
amap[ila, ilo] = copy.deepcopy(profiles[key])
f.close()
if self._directory is not None:
path_elevation = os.path.join(self._directory, Crust2.fn_elevation)
f = open(path_elevation, 'r')
else:
f = StringIO(decode(elevation))
f.readline()
for ila, line in enumerate(f):
for ilo, s in enumerate(line.split()[1:]):
p = amap[ila, ilo]
p.set_elevation(float(s))
if p.elevation() < 0.:
p.set_layer_thickness(LWATER, -p.elevation())
f.close()
self._typemap = amap
[docs] @staticmethod
def instance():
'''Get the global default Crust2 instance.'''
if Crust2._instance is None:
Crust2._instance = Crust2()
return Crust2._instance
[docs]def get_profile_keys():
'''Get list of all profile keys.'''
crust2 = Crust2.instance()
return list(crust2.profile_keys)
[docs]def get_profile(*args, **kwargs):
'''Get Crust2x2 profile for given location or profile key.
Get profile for (lat,lon) or raw profile for given string key.
'''
crust2 = Crust2.instance()
return crust2.get_profile(*args, **kwargs)
[docs]def plot_crustal_thickness(crust2=None, filename='crustal_thickness.pdf'):
'''
Create a quick and dirty plot of the crustal thicknesses defined in
CRUST2.0.
'''
if crust2 is None:
crust2 = Crust2.instance()
def func(lat, lon):
return crust2.get_profile(lat, lon).crustal_thickness(),
plot(func, filename, zscaled_unit='km', zscaled_unit_factor=0.001)
[docs]def plot_vp_belowcrust(crust2=None, filename='vp_below_crust.pdf'):
'''
Create a quick and dirty plot of vp below the crust, as defined in
CRUST2.0.
'''
if crust2 is None:
crust2 = Crust2.instance()
def func(lat, lon):
return crust2.get_profile(lat, lon).get_layer(LBELOWCRUST)[1]
plot(func, filename, zscaled_unit='km/s', zscaled_unit_factor=0.001)
[docs]def plot(func, filename, **kwargs):
nlats, nlons = 91, 181
lats = num.linspace(-90., 90., nlats)
lons = num.linspace(-180., 180., nlons)
vecfunc = num.vectorize(func, [num.float])
latss, lonss = num.meshgrid(lats, lons)
thickness = vecfunc(latss, lonss)
from pyrocko.plot import gmtpy
cm = gmtpy.cm
marg = (1.5*cm, 2.5*cm, 1.5*cm, 1.5*cm)
p = gmtpy.Simple(
width=20*cm, height=10*cm, margins=marg,
with_palette=True, **kwargs)
p.density_plot(gmtpy.tabledata(lons, lats, thickness.T))
p.save(filename)
