1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import
7import os.path as op
8import math
9from collections import defaultdict
11import numpy as num
13from pyrocko import util, config
14from pyrocko import orthodrome as od
15from pyrocko.guts_array import Array
16from pyrocko.guts import Object, String
17from .util import get_download_callback
20PI = math.pi
23class Plate(Object):
25 name = String.T(
26 help='Name of the tectonic plate.')
27 points = Array.T(
28 dtype=float, shape=(None, 2),
29 help='Points on the plate.')
31 def max_interpoint_distance(self):
32 p = od.latlon_to_xyz(self.points)
33 return math.sqrt(num.max(num.sum(
34 (p[num.newaxis, :, :] - p[:, num.newaxis, :])**2, axis=2)))
36 def contains_point(self, point):
37 return od.contains_point(self.points, point)
39 def contains_points(self, points):
40 return od.contains_points(self.points, points)
43class Boundary(Object):
45 plate_name1 = String.T()
46 plate_name2 = String.T()
47 kind = String.T()
48 points = Array.T(dtype=float, shape=(None, 2))
49 cpoints = Array.T(dtype=float, shape=(None, 2))
50 itypes = Array.T(dtype=int, shape=(None))
52 def split_types(self, groups=None):
53 xyz = od.latlon_to_xyz(self.points)
54 xyzmid = (xyz[1:] + xyz[:-1, :]) * 0.5
55 cxyz = od.latlon_to_xyz(self.cpoints)
56 d = od.distances3d(xyzmid[num.newaxis, :, :], cxyz[:, num.newaxis, :])
57 idmin = num.argmin(d, axis=0)
58 itypes = self.itypes[idmin]
60 if groups is None:
61 groupmap = num.arange(len(self._index_to_type))
62 else:
63 d = {}
64 for igroup, group in enumerate(groups):
65 for name in group:
66 d[name] = igroup
68 groupmap = num.array(
69 [d[name] for name in self._index_to_type],
70 dtype=int)
72 iswitch = num.concatenate(
73 ([0],
74 num.where(groupmap[itypes[1:]] != groupmap[itypes[:-1]])[0]+1,
75 [itypes.size]))
77 results = []
78 for ii in range(iswitch.size-1):
79 if groups is not None:
80 tt = [self._index_to_type[ityp] for ityp in num.unique(
81 itypes[iswitch[ii]:iswitch[ii+1]])]
82 else:
83 tt = self._index_to_type[itypes[iswitch[ii]]]
85 results.append((tt, self.points[iswitch[ii]:iswitch[ii+1]+1]))
87 return results
90class Dataset(object):
92 def __init__(self, name, data_dir, citation):
93 self.name = name
94 self.citation = citation
95 self.data_dir = data_dir
97 def fpath(self, filename):
98 return op.join(self.data_dir, filename)
100 def download_file(
101 self, url, fpath, username=None, password=None,
102 status_callback=None):
104 util.download_file(
105 url, fpath, username, password, status_callback=status_callback)
108class PlatesDataset(Dataset):
109 pass
112class PeterBird2003(PlatesDataset):
113 '''
114 An updated digital model of plate boundaries.
115 '''
116 __citation = '''
117 Bird, Peter. "An updated digital model of plate boundaries." Geochemistry,
118 Geophysics, Geosystems 4.3 (2003).
119 '''
121 def __init__(
122 self,
123 name='PeterBird2003',
124 data_dir=None,
125 raw_data_url=('https://mirror.pyrocko.org/peterbird.name/oldFTP/PB2002/%s')): # noqa
127 if data_dir is None:
128 data_dir = op.join(config.config().tectonics_dir, name)
130 PlatesDataset.__init__(
131 self,
132 name,
133 data_dir=data_dir,
134 citation=self.__citation)
136 self.raw_data_url = raw_data_url
138 self.filenames = [
139 '2001GC000252_readme.txt',
140 'PB2002_boundaries.dig.txt',
141 'PB2002_orogens.dig.txt',
142 'PB2002_plates.dig.txt',
143 'PB2002_poles.dat.txt',
144 'PB2002_steps.dat.txt']
146 self._full_names = None
148 def full_name(self, name):
149 if not self._full_names:
150 fn = util.data_file(op.join('tectonics', 'bird2003_plates.txt'))
152 with open(fn, 'r') as f:
153 self._full_names = dict(
154 line.strip().split(None, 1) for line in f)
156 return self._full_names[name]
158 def download_if_needed(self):
159 for fn in self.filenames:
160 fpath = self.fpath(fn)
161 if not op.exists(fpath):
162 self.download_file(
163 self.raw_data_url % fn, fpath,
164 status_callback=get_download_callback(
165 'Downloading Bird 2003 plate database...'))
167 def get_boundaries(self):
168 self.download_if_needed()
169 fpath = self.fpath('PB2002_steps.dat.txt')
171 d = defaultdict(list)
172 ntyp = 0
173 type_to_index = {}
174 index_to_type = []
175 with open(fpath, 'rb') as f:
176 data = []
177 for line in f:
178 lsplit = line.split()
179 name = lsplit[1].lstrip(b':').decode('ascii')
180 plate_name1 = str(name[0:2])
181 plate_name2 = str(name[3:5])
182 kind = name[2]
184 alon, alat, blon, blat = list(map(float, lsplit[2:6]))
185 mlat = (alat + blat) * 0.5
186 dlon = ((blon - alon) + 180.) % 360. - 180.
187 mlon = alon + dlon * 0.5
188 typ = str(lsplit[14].strip(b':*').decode('ascii'))
190 if typ not in type_to_index:
191 ntyp += 1
192 type_to_index[typ] = ntyp - 1
193 index_to_type.append(typ)
195 ityp = type_to_index[typ]
196 d[plate_name1, kind, plate_name2].append((mlat, mlon, ityp))
198 d2 = {}
199 for k in d:
200 d2[k] = (
201 num.array([x[:2] for x in d[k]], dtype=float),
202 num.array([x[2] for x in d[k]], dtype=int))
204 fpath = self.fpath('PB2002_boundaries.dig.txt')
205 boundaries = []
206 plate_name1 = ''
207 plate_name2 = ''
208 kind = '-'
209 with open(fpath, 'rb') as f:
210 data = []
211 for line in f:
212 if line.startswith(b'***'):
213 cpoints, itypes = d2[plate_name1, kind, plate_name2]
214 boundaries.append(Boundary(
215 plate_name1=plate_name1,
216 plate_name2=plate_name2,
217 kind=kind,
218 points=num.array(data, dtype=float),
219 cpoints=cpoints,
220 itypes=itypes))
222 boundaries[-1]._type_to_index = type_to_index
223 boundaries[-1]._index_to_type = index_to_type
225 data = []
226 elif line.startswith(b' '):
227 data.append(list(map(float, line.split(b',')))[::-1])
228 else:
229 s = line.strip().decode('ascii')
230 plate_name1 = str(s[0:2])
231 plate_name2 = str(s[3:5])
232 kind = s[2]
234 return boundaries
236 def get_plates(self):
237 self.download_if_needed()
238 fpath = self.fpath('PB2002_plates.dig.txt')
239 plates = []
240 name = ''
241 with open(fpath, 'rb') as f:
242 data = []
243 for line in f:
244 if line.startswith(b'***'):
245 plates.append(Plate(
246 name=name,
247 points=num.array(data, dtype=float)))
249 data = []
250 elif line.startswith(b' '):
251 data.append(list(map(float, line.split(b',')))[::-1])
252 else:
253 name = str(line.strip().decode('ascii'))
255 return plates
258class StrainRateDataset(Dataset):
259 pass
262class GSRM1(StrainRateDataset):
263 '''
264 Global Strain Rate Map. An integrated global model of present-day
265 plate motions and plate boundary deformation'''
266 __citation = '''Kreemer, C., W.E. Holt, and A.J. Haines, "An integrated
267 global model of present-day plate motions and plate boundary deformation",
268 Geophys. J. Int., 154, 8-34, 2003.
269 '''
271 def __init__(
272 self,
273 name='GSRM1.2',
274 data_dir=None,
275 raw_data_url=('https://mirror.pyrocko.org/gsrm.unavco.org/model'
276 '/files/1.2/%s')):
278 if data_dir is None:
279 data_dir = op.join(config.config().tectonics_dir, name)
281 StrainRateDataset.__init__(
282 self,
283 name,
284 data_dir=data_dir,
285 citation=self.__citation)
287 self.raw_data_url = raw_data_url
288 self._full_names = None
289 self._names = None
291 def full_names(self):
292 if not self._full_names:
293 fn = util.data_file(op.join('tectonics', 'gsrm1_plates.txt'))
295 with open(fn, 'r') as f:
296 self._full_names = dict(
297 line.strip().split(None, 1) for line in f)
299 return self._full_names
301 def full_name(self, name):
302 name = self.plate_alt_names().get(name, name)
303 return self.full_names()[name]
305 def plate_names(self):
306 if self._names is None:
307 self._names = sorted(self.full_names().keys())
309 return self._names
311 def plate_alt_names(self):
312 # 'African Plate' is named 'Nubian Plate' in GSRM1
313 return {'AF': 'NU'}
315 def download_if_needed(self, fn):
316 fpath = self.fpath(fn)
317 if not op.exists(fpath):
318 self.download_file(
319 self.raw_data_url % fn, fpath,
320 status_callback=get_download_callback(
321 'Downloading global strain rate map...'))
323 def get_velocities(self, reference_name=None, region=None):
324 reference_name = self.plate_alt_names().get(
325 reference_name, reference_name)
327 if reference_name is None:
328 reference_name = 'NNR'
330 fn = 'velocity_%s.dat' % reference_name
332 self.download_if_needed(fn)
333 fpath = self.fpath(fn)
334 data = []
335 with open(fpath, 'rb') as f:
336 for line in f:
337 if line.strip().startswith(b'#'):
338 continue
339 t = line.split()
340 data.append(list(map(float, t)))
342 arr = num.array(data, dtype=float)
344 if region is not None:
345 points = arr[:, 1::-1]
346 mask = od.points_in_region(points, region)
347 arr = arr[mask, :]
349 lons, lats, veast, vnorth, veast_err, vnorth_err, corr = arr.T
350 return lats, lons, vnorth, veast, vnorth_err, veast_err, corr
353__all__ = '''
354GSRM1
355PeterBird2003
356Plate'''.split()