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'''
69from __future__ import absolute_import, division
71import os
72import copy
73import math
74from io import StringIO
76import numpy as num
79LICE, LWATER, LSOFTSED, LHARDSED, LUPPERCRUST, LMIDDLECRUST, \
80 LLOWERCRUST, LBELOWCRUST = list(range(8))
83class Crust2Profile(object):
84 '''
85 Representation of a CRUST2.0 key profile.
86 '''
88 layer_names = (
89 'ice', 'water', 'soft sed.', 'hard sed.', 'upper crust',
90 'middle crust', 'lower crust', 'mantle')
92 def __init__(self, ident, name, vp, vs, rho, thickness, elevation):
93 self._ident = ident
94 self._name = name
95 self._vp = vp
96 self._vs = vs
97 self._rho = rho
98 self._thickness = thickness
99 self._crustal_thickness = None
100 self._elevation = elevation
102 def get_weeded(self, include_waterlayer=False):
103 '''
104 Get layers used in the profile.
106 :param include_waterlayer: include water layer if ``True``. Default is
107 ``False``
109 :returns: NumPy array with rows ``depth``, ``vp``, ``vs``, ``density``
110 '''
111 depth = 0.
112 layers = []
113 for ilayer, thickness, vp, vs, rho in zip(
114 range(8),
115 self._thickness,
116 self._vp[:-1],
117 self._vs[:-1],
118 self._rho[:-1]):
120 if thickness == 0.0:
121 continue
123 if not include_waterlayer and ilayer == LWATER:
124 continue
126 layers.append([depth, vp, vs, rho])
127 layers.append([depth+thickness, vp, vs, rho])
128 depth += thickness
130 layers.append([
131 depth,
132 self._vp[LBELOWCRUST],
133 self._vs[LBELOWCRUST],
134 self._rho[LBELOWCRUST]])
136 return num.array(layers).T
138 def get_layer(self, ilayer):
139 '''
140 Get parameters for a layer.
142 :param ilayer: id of layer
143 :returns: thickness, vp, vs, density
144 '''
146 if ilayer == LBELOWCRUST:
147 thickness = num.Inf
148 else:
149 thickness = self._thickness[ilayer]
151 return thickness, self._vp[ilayer], self._vs[ilayer], self._rho[ilayer]
153 def set_elevation(self, elevation):
154 self._elevation = elevation
156 def set_layer_thickness(self, ilayer, thickness):
157 self._thickness[ilayer] = thickness
159 def elevation(self):
160 return self._elevation
162 def __str__(self):
164 vvp, vvs, vrho, vthi = self.averages()
166 return '''type, name: %s, %s
167elevation: %15.5g
168crustal thickness: %15.5g
169average vp, vs, rho: %15.5g %15.5g %15.5g
170mantle ave. vp, vs, rho: %15.5g %15.5g %15.5g
172%s''' % (self._ident, self._name, self._elevation, vthi, vvp, vvs, vrho,
173 self._vp[LBELOWCRUST], self._vs[LBELOWCRUST], self._rho[LBELOWCRUST],
174 '\n'.join([
175 '%15.5g %15.5g %15.5g %15.5g %s' % x
176 for x in zip(
177 self._thickness, self._vp[:-1], self._vs[:-1], self._rho[:-1],
178 Crust2Profile.layer_names)]))
180 def crustal_thickness(self):
181 '''
182 Get total crustal thickness.
184 Takes into account ice layer.
185 Does not take into account water layer.
186 '''
188 return num.sum(self._thickness[3:]) + self._thickness[LICE]
190 def averages(self):
191 '''
192 Get crustal averages for vp, vs, density and total crustal thickness.
194 Takes into account ice layer.
195 Does not take into account water layer.
196 '''
198 vthi = self.crustal_thickness()
199 vvp = num.sum(self._thickness[3:] / self._vp[3:-1]) + \
200 self._thickness[LICE] / self._vp[LICE]
201 vvs = num.sum(self._thickness[3:] / self._vs[3:-1]) + \
202 self._thickness[LICE] / self._vs[LICE]
203 vrho = num.sum(self._thickness[3:] * self._rho[3:-1]) + \
204 self._thickness[LICE] * self._rho[LICE]
206 vvp = vthi / vvp
207 vvs = vthi / vvs
208 vrho = vrho / vthi
210 return vvp, vvs, vrho, vthi
213def _sa2arr(sa):
214 return num.array([float(x) for x in sa], dtype=float)
217def _wrap(x, mi, ma):
218 if mi <= x and x <= ma:
219 return x
221 return x - math.floor((x-mi)/(ma-mi)) * (ma-mi)
224def _clip(x, mi, ma):
225 return min(max(mi, x), ma)
228class Crust2(object):
229 '''
230 Access CRUST2.0 model.
232 :param directory: Directory with the data files which contain the
233 CRUST2.0 model data. If this is set to ``None``, builtin CRUST2.0
234 files are used.
235 '''
237 fn_keys = 'CNtype2_key.txt'
238 fn_elevation = 'CNelevatio2.txt'
239 fn_map = 'CNtype2.txt'
241 nlo = 180
242 nla = 90
244 _instance = None
246 def __init__(self, directory=None):
248 self.profile_keys = []
249 self._directory = directory
250 self._typemap = None
251 self._load_crustal_model()
253 def get_profile(self, *args, **kwargs):
254 '''
255 Get crustal profile at a specific location or raw profile for given
256 key.
258 Get profile for location ``(lat, lon)``, or raw profile for given
259 string key.
261 :rtype: instance of :py:class:`Crust2Profile`
262 '''
264 lat = kwargs.pop('lat', None)
265 lon = kwargs.pop('lon', None)
267 if len(args) == 2:
268 lat, lon = args
270 if lat is not None and lon is not None:
271 return self._typemap[self._indices(float(lat), float(lon))]
272 else:
273 return self._raw_profiles[args[0]]
275 def _indices(self, lat, lon):
276 lat = _clip(lat, -90., 90.)
277 lon = _wrap(lon, -180., 180.)
278 dlo = 360./Crust2.nlo
279 dla = 180./Crust2.nla
280 cola = 90.-lat
281 ilat = _clip(int(cola/dla), 0, Crust2.nla-1)
282 ilon = int((lon+180.)/dlo) % Crust2.nlo
283 return ilat, ilon
285 def _load_crustal_model(self):
287 if self._directory is not None:
288 path_keys = os.path.join(self._directory, Crust2.fn_keys)
289 f = open(path_keys, 'r')
290 else:
291 from .crust2x2_data import decode, type2_key, type2, elevation
293 f = StringIO(decode(type2_key))
295 # skip header
296 for i in range(5):
297 f.readline()
299 profiles = {}
300 while True:
301 line = f.readline()
302 if not line:
303 break
304 ident, name = line.split(None, 1)
305 line = f.readline()
306 vp = _sa2arr(line.split()) * 1000.
307 line = f.readline()
308 vs = _sa2arr(line.split()) * 1000.
309 line = f.readline()
310 rho = _sa2arr(line.split()) * 1000.
311 line = f.readline()
312 toks = line.split()
313 thickness = _sa2arr(toks[:-2]) * 1000.
315 assert ident not in profiles
317 profiles[ident] = Crust2Profile(
318 ident.strip(), name.strip(), vp, vs, rho, thickness, 0.0)
320 f.close()
322 self._raw_profiles = profiles
323 self.profile_keys = sorted(profiles.keys())
325 if self._directory is not None:
326 path_map = os.path.join(self._directory, Crust2.fn_map)
327 f = open(path_map, 'r')
328 else:
329 f = StringIO(decode(type2))
331 f.readline() # header
333 amap = {}
334 for ila, line in enumerate(f):
335 keys = line.split()[1:]
336 for ilo, key in enumerate(keys):
337 amap[ila, ilo] = copy.deepcopy(profiles[key])
339 f.close()
341 if self._directory is not None:
342 path_elevation = os.path.join(self._directory, Crust2.fn_elevation)
343 f = open(path_elevation, 'r')
345 else:
346 f = StringIO(decode(elevation))
348 f.readline()
349 for ila, line in enumerate(f):
350 for ilo, s in enumerate(line.split()[1:]):
351 p = amap[ila, ilo]
352 p.set_elevation(float(s))
353 if p.elevation() < 0.:
354 p.set_layer_thickness(LWATER, -p.elevation())
356 f.close()
358 self._typemap = amap
360 @staticmethod
361 def instance():
362 '''
363 Get the global default Crust2 instance.
364 '''
366 if Crust2._instance is None:
367 Crust2._instance = Crust2()
369 return Crust2._instance
372def get_profile_keys():
373 '''
374 Get list of all profile keys.
375 '''
377 crust2 = Crust2.instance()
378 return list(crust2.profile_keys)
381def get_profile(*args, **kwargs):
382 '''
383 Get Crust2x2 profile for given location or profile key.
385 Get profile for (lat,lon) or raw profile for given string key.
386 '''
388 crust2 = Crust2.instance()
389 return crust2.get_profile(*args, **kwargs)
392def plot_crustal_thickness(crust2=None, filename='crustal_thickness.pdf'):
393 '''
394 Create a quick and dirty plot of the crustal thicknesses defined in
395 CRUST2.0.
396 '''
398 if crust2 is None:
399 crust2 = Crust2.instance()
401 def func(lat, lon):
402 return crust2.get_profile(lat, lon).crustal_thickness(),
404 plot(func, filename, zscaled_unit='km', zscaled_unit_factor=0.001)
407def plot_vp_belowcrust(crust2=None, filename='vp_below_crust.pdf'):
408 '''
409 Create a quick and dirty plot of vp below the crust, as defined in
410 CRUST2.0.
411 '''
413 if crust2 is None:
414 crust2 = Crust2.instance()
416 def func(lat, lon):
417 return crust2.get_profile(lat, lon).get_layer(LBELOWCRUST)[1]
419 plot(func, filename, zscaled_unit='km/s', zscaled_unit_factor=0.001)
422def plot(func, filename, **kwargs):
423 nlats, nlons = 91, 181
424 lats = num.linspace(-90., 90., nlats)
425 lons = num.linspace(-180., 180., nlons)
427 vecfunc = num.vectorize(func, [float])
428 latss, lonss = num.meshgrid(lats, lons)
429 thickness = vecfunc(latss, lonss)
431 from pyrocko.plot import gmtpy
432 cm = gmtpy.cm
433 marg = (1.5*cm, 2.5*cm, 1.5*cm, 1.5*cm)
434 p = gmtpy.Simple(
435 width=20*cm, height=10*cm, margins=marg,
436 with_palette=True, **kwargs)
438 p.density_plot(gmtpy.tabledata(lons, lats, thickness.T))
439 p.save(filename)