1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

# 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 

 

from builtins import range 

 

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)) 

 

 

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 

 

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 

 

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] 

 

def set_elevation(self, elevation): 

self._elevation = elevation 

 

def set_layer_thickness(self, ilayer, thickness): 

self._thickness[ilayer] = thickness 

 

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)])) 

 

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] 

 

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) 

 

 

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() 

 

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 

 

@staticmethod 

def instance(): 

'''Get the global default Crust2 instance.''' 

 

if Crust2._instance is None: 

Crust2._instance = Crust2() 

 

return Crust2._instance 

 

 

def get_profile_keys(): 

'''Get list of all profile keys.''' 

 

crust2 = Crust2.instance() 

return list(crust2.profile_keys) 

 

 

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) 

 

 

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) 

 

 

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) 

 

 

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)