1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6Interface to use CRUST2.0 model by Laske, Masters and Reif.
8All functions defined in this module return SI units (m, m/s, kg/m^3).
10.. note::
11 Please refer to the REM web site if you use this model:
13 http://igppweb.ucsd.edu/~gabi/rem.html
15 or
17 Bassin, C., Laske, G. and Masters, G., The Current Limits of Resolution for
18 Surface Wave Tomography in North America, EOS Trans AGU, 81, F897, 2000. A
19 description of CRUST 5.1 can be found in: Mooney, Laske and Masters, Crust
20 5.1: a global crustal model at 5x5 degrees, JGR, 103, 727-747, 1998.
22Usage
23-----
25::
27 >>> from pyrocko import crust2x2
28 >>> p = crust2x2.get_profile(10., 20.)
29 >>> print p
30 type, name: G2, Archean 0.5 km seds.
31 elevation: 529
32 crustal thickness: 38500
33 average vp, vs, rho: 6460.7 3665.1 2867.5
34 mantle ave. vp, vs, rho: 8200 4700 3400
36 0 3810 1940 920 ice
37 0 1500 0 1020 water
38 500 2500 1200 2100 soft sed.
39 0 4000 2100 2400 hard sed.
40 12500 6200 3600 2800 upper crust
41 13000 6400 3600 2850 middle crust
42 13000 6800 3800 2950 lower crust
43 >>> print p.get_weeded()
44 [[ 0. 500. 500. 13000. 13000. 26000. 26000. 39000. 39000.]
45 [ 2500. 2500. 6200. 6200. 6400. 6400. 6800. 6800. 8200.]
46 [ 1200. 1200. 3600. 3600. 3600. 3600. 3800. 3800. 4700.]
47 [ 2100. 2100. 2800. 2800. 2850. 2850. 2950. 2950. 3400.]]
50Constants
51---------
53============== ==============
54Layer id Layer name
55============== ==============
56LICE ice
57LWATER water
58LSOFTSED soft sediments
59LHARDSED hard sediments
60LUPPERCRUST upper crust
61LMIDDLECRUST middle crust
62LLOWERCRUST lower crust
63LBELOWCRUST below crust
64============== ==============
66Contents
67--------
68'''
70import os
71import copy
72import math
73from io import StringIO
75import numpy as num
78LICE, LWATER, LSOFTSED, LHARDSED, LUPPERCRUST, LMIDDLECRUST, \
79 LLOWERCRUST, LBELOWCRUST = list(range(8))
82class Crust2Profile(object):
83 '''
84 Representation of a CRUST2.0 key profile.
85 '''
87 layer_names = (
88 'ice', 'water', 'soft sed.', 'hard sed.', 'upper crust',
89 'middle crust', 'lower crust', 'mantle')
91 def __init__(self, ident, name, vp, vs, rho, thickness, elevation):
92 self._ident = ident
93 self._name = name
94 self._vp = vp
95 self._vs = vs
96 self._rho = rho
97 self._thickness = thickness
98 self._crustal_thickness = None
99 self._elevation = elevation
101 def get_weeded(self, include_waterlayer=False):
102 '''
103 Get layers used in the profile.
105 :param include_waterlayer: include water layer if ``True``. Default is
106 ``False``
108 :returns: NumPy array with rows ``depth``, ``vp``, ``vs``, ``density``
109 '''
110 depth = 0.
111 layers = []
112 for ilayer, thickness, vp, vs, rho in zip(
113 range(8),
114 self._thickness,
115 self._vp[:-1],
116 self._vs[:-1],
117 self._rho[:-1]):
119 if thickness == 0.0:
120 continue
122 if not include_waterlayer and ilayer == LWATER:
123 continue
125 layers.append([depth, vp, vs, rho])
126 layers.append([depth+thickness, vp, vs, rho])
127 depth += thickness
129 layers.append([
130 depth,
131 self._vp[LBELOWCRUST],
132 self._vs[LBELOWCRUST],
133 self._rho[LBELOWCRUST]])
135 return num.array(layers).T
137 def get_layer(self, ilayer):
138 '''
139 Get parameters for a layer.
141 :param ilayer: id of layer
142 :returns: thickness, vp, vs, density
143 '''
145 if ilayer == LBELOWCRUST:
146 thickness = num.Inf
147 else:
148 thickness = self._thickness[ilayer]
150 return thickness, self._vp[ilayer], self._vs[ilayer], self._rho[ilayer]
152 def set_elevation(self, elevation):
153 self._elevation = elevation
155 def set_layer_thickness(self, ilayer, thickness):
156 self._thickness[ilayer] = thickness
158 def elevation(self):
159 return self._elevation
161 def __str__(self):
163 vvp, vvs, vrho, vthi = self.averages()
165 return '''type, name: %s, %s
166elevation: %15.5g
167crustal thickness: %15.5g
168average vp, vs, rho: %15.5g %15.5g %15.5g
169mantle ave. vp, vs, rho: %15.5g %15.5g %15.5g
171%s''' % (self._ident, self._name, self._elevation, vthi, vvp, vvs, vrho,
172 self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST],
173 '\n'.join([
174 '%15.5g %15.5g %15.5g %15.5g %s' % x
175 for x in zip(
176 self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1],
177 Crust2Profile.layer_names)]))
179 def crustal_thickness(self):
180 '''
181 Get total crustal thickness.
183 Takes into account ice layer.
184 Does not take into account water layer.
185 '''
187 return num.sum(self._thickness[3:]) + self._thickness[LICE]
189 def averages(self):
190 '''
191 Get crustal averages for vp, vs, density and total crustal thickness.
193 Takes into account ice layer.
194 Does not take into account water layer.
195 '''
197 vthi = self.crustal_thickness()
198 vvp = num.sum(self._thickness[3:] / self._vp[3:-1]) + \
199 self._thickness[LICE] / self._vp[LICE]
200 vvs = num.sum(self._thickness[3:] / self._vs[3:-1]) + \
201 self._thickness[LICE] / self._vs[LICE]
202 vrho = num.sum(self._thickness[3:] * self._rho[3:-1]) + \
203 self._thickness[LICE] * self._rho[LICE]
205 vvp = vthi / vvp
206 vvs = vthi / vvs
207 vrho = vrho / vthi
209 return vvp, vvs, vrho, vthi
212def _sa2arr(sa):
213 return num.array([float(x) for x in sa], dtype=float)
216def _wrap(x, mi, ma):
217 if mi <= x and x <= ma:
218 return x
220 return x - math.floor((x-mi)/(ma-mi)) * (ma-mi)
223def _clip(x, mi, ma):
224 return min(max(mi, x), ma)
227class Crust2(object):
228 '''
229 Access CRUST2.0 model.
231 :param directory: Directory with the data files which contain the
232 CRUST2.0 model data. If this is set to ``None``, builtin CRUST2.0
233 files are used.
234 '''
236 fn_keys = 'CNtype2_key.txt'
237 fn_elevation = 'CNelevatio2.txt'
238 fn_map = 'CNtype2.txt'
240 nlo = 180
241 nla = 90
243 _instance = None
245 def __init__(self, directory=None):
247 self.profile_keys = []
248 self._directory = directory
249 self._typemap = None
250 self._load_crustal_model()
252 def get_profile(self, *args, **kwargs):
253 '''
254 Get crustal profile at a specific location or raw profile for given
255 key.
257 Get profile for location ``(lat, lon)``, or raw profile for given
258 string key.
260 :rtype: instance of :py:class:`Crust2Profile`
261 '''
263 lat = kwargs.pop('lat', None)
264 lon = kwargs.pop('lon', None)
266 if len(args) == 2:
267 lat, lon = args
269 if lat is not None and lon is not None:
270 return self._typemap[self._indices(float(lat), float(lon))]
271 else:
272 return self._raw_profiles[args[0]]
274 def _indices(self, lat, lon):
275 lat = _clip(lat, -90., 90.)
276 lon = _wrap(lon, -180., 180.)
277 dlo = 360./Crust2.nlo
278 dla = 180./Crust2.nla
279 cola = 90.-lat
280 ilat = _clip(int(cola/dla), 0, Crust2.nla-1)
281 ilon = int((lon+180.)/dlo) % Crust2.nlo
282 return ilat, ilon
284 def _load_crustal_model(self):
286 if self._directory is not None:
287 path_keys = os.path.join(self._directory, Crust2.fn_keys)
288 f = open(path_keys, 'r')
289 else:
290 from .crust2x2_data import decode, type2_key, type2, elevation
292 f = StringIO(decode(type2_key))
294 # skip header
295 for i in range(5):
296 f.readline()
298 profiles = {}
299 while True:
300 line = f.readline()
301 if not line:
302 break
303 ident, name = line.split(None, 1)
304 line = f.readline()
305 vp = _sa2arr(line.split()) * 1000.
306 line = f.readline()
307 vs = _sa2arr(line.split()) * 1000.
308 line = f.readline()
309 rho = _sa2arr(line.split()) * 1000.
310 line = f.readline()
311 toks = line.split()
312 thickness = _sa2arr(toks[:-2]) * 1000.
314 assert ident not in profiles
316 profiles[ident] = Crust2Profile(
317 ident.strip(), name.strip(), vp, vs, rho, thickness, 0.0)
319 f.close()
321 self._raw_profiles = profiles
322 self.profile_keys = sorted(profiles.keys())
324 if self._directory is not None:
325 path_map = os.path.join(self._directory, Crust2.fn_map)
326 f = open(path_map, 'r')
327 else:
328 f = StringIO(decode(type2))
330 f.readline() # header
332 amap = {}
333 for ila, line in enumerate(f):
334 keys = line.split()[1:]
335 for ilo, key in enumerate(keys):
336 amap[ila, ilo] = copy.deepcopy(profiles[key])
338 f.close()
340 if self._directory is not None:
341 path_elevation = os.path.join(self._directory, Crust2.fn_elevation)
342 f = open(path_elevation, 'r')
344 else:
345 f = StringIO(decode(elevation))
347 f.readline()
348 for ila, line in enumerate(f):
349 for ilo, s in enumerate(line.split()[1:]):
350 p = amap[ila, ilo]
351 p.set_elevation(float(s))
352 if p.elevation() < 0.:
353 p.set_layer_thickness(LWATER, -p.elevation())
355 f.close()
357 self._typemap = amap
359 @staticmethod
360 def instance():
361 '''
362 Get the global default Crust2 instance.
363 '''
365 if Crust2._instance is None:
366 Crust2._instance = Crust2()
368 return Crust2._instance
371def get_profile_keys():
372 '''
373 Get list of all profile keys.
374 '''
376 crust2 = Crust2.instance()
377 return list(crust2.profile_keys)
380def get_profile(*args, **kwargs):
381 '''
382 Get Crust2x2 profile for given location or profile key.
384 Get profile for (lat,lon) or raw profile for given string key.
385 '''
387 crust2 = Crust2.instance()
388 return crust2.get_profile(*args, **kwargs)
391def plot_crustal_thickness(crust2=None, filename='crustal_thickness.pdf'):
392 '''
393 Create a quick and dirty plot of the crustal thicknesses defined in
394 CRUST2.0.
395 '''
397 if crust2 is None:
398 crust2 = Crust2.instance()
400 def func(lat, lon):
401 return crust2.get_profile(lat, lon).crustal_thickness(),
403 plot(func, filename, zscaled_unit='km', zscaled_unit_factor=0.001)
406def plot_vp_belowcrust(crust2=None, filename='vp_below_crust.pdf'):
407 '''
408 Create a quick and dirty plot of vp below the crust, as defined in
409 CRUST2.0.
410 '''
412 if crust2 is None:
413 crust2 = Crust2.instance()
415 def func(lat, lon):
416 return crust2.get_profile(lat, lon).get_layer(LBELOWCRUST)[1]
418 plot(func, filename, zscaled_unit='km/s', zscaled_unit_factor=0.001)
421def plot(func, filename, **kwargs):
422 nlats, nlons = 91, 181
423 lats = num.linspace(-90., 90., nlats)
424 lons = num.linspace(-180., 180., nlons)
426 vecfunc = num.vectorize(func, [float])
427 latss, lonss = num.meshgrid(lats, lons)
428 thickness = vecfunc(latss, lonss)
430 from pyrocko.plot import gmtpy
431 cm = gmtpy.cm
432 marg = (1.5*cm, 2.5*cm, 1.5*cm, 1.5*cm)
433 p = gmtpy.Simple(
434 width=20*cm, height=10*cm, margins=marg,
435 with_palette=True, **kwargs)
437 p.density_plot(gmtpy.tabledata(lons, lats, thickness.T))
438 p.save(filename)