1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import os.path as op
7import math
8from collections import defaultdict
10import numpy as num
12from pyrocko import util, config
13from pyrocko import orthodrome as od
14from pyrocko.guts_array import Array
15from pyrocko.guts import Object, String
16from .util import get_download_callback
19PI = math.pi
22class Plate(Object):
24 name = String.T(
25 help='Name of the tectonic plate.')
26 points = Array.T(
27 dtype=float, shape=(None, 2),
28 help='Points on the plate.')
30 def max_interpoint_distance(self):
31 p = od.latlon_to_xyz(self.points)
32 return math.sqrt(num.max(num.sum(
33 (p[num.newaxis, :, :] - p[:, num.newaxis, :])**2, axis=2)))
35 def contains_point(self, point):
36 return od.contains_point(self.points, point)
38 def contains_points(self, points):
39 return od.contains_points(self.points, points)
42class Boundary(Object):
44 plate_name1 = String.T()
45 plate_name2 = String.T()
46 kind = String.T()
47 points = Array.T(dtype=float, shape=(None, 2))
48 cpoints = Array.T(dtype=float, shape=(None, 2))
49 itypes = Array.T(dtype=int, shape=(None))
51 def split_types(self, groups=None):
52 xyz = od.latlon_to_xyz(self.points)
53 xyzmid = (xyz[1:] + xyz[:-1, :]) * 0.5
54 cxyz = od.latlon_to_xyz(self.cpoints)
55 d = od.distances3d(xyzmid[num.newaxis, :, :], cxyz[:, num.newaxis, :])
56 idmin = num.argmin(d, axis=0)
57 itypes = self.itypes[idmin]
59 if groups is None:
60 groupmap = num.arange(len(self._index_to_type))
61 else:
62 d = {}
63 for igroup, group in enumerate(groups):
64 for name in group:
65 d[name] = igroup
67 groupmap = num.array(
68 [d[name] for name in self._index_to_type],
69 dtype=int)
71 iswitch = num.concatenate(
72 ([0],
73 num.where(groupmap[itypes[1:]] != groupmap[itypes[:-1]])[0]+1,
74 [itypes.size]))
76 results = []
77 for ii in range(iswitch.size-1):
78 if groups is not None:
79 tt = [self._index_to_type[ityp] for ityp in num.unique(
80 itypes[iswitch[ii]:iswitch[ii+1]])]
81 else:
82 tt = self._index_to_type[itypes[iswitch[ii]]]
84 results.append((tt, self.points[iswitch[ii]:iswitch[ii+1]+1]))
86 return results
89class Dataset(object):
91 def __init__(self, name, data_dir, citation):
92 self.name = name
93 self.citation = citation
94 self.data_dir = data_dir
96 def fpath(self, filename):
97 return op.join(self.data_dir, filename)
99 def download_file(
100 self, url, fpath, username=None, password=None,
101 status_callback=None):
103 util.download_file(
104 url, fpath, username, password, status_callback=status_callback)
107class PlatesDataset(Dataset):
108 pass
111class PeterBird2003(PlatesDataset):
112 '''
113 An updated digital model of plate boundaries.
114 '''
115 __citation = '''
116 Bird, Peter. "An updated digital model of plate boundaries." Geochemistry,
117 Geophysics, Geosystems 4.3 (2003).
118 '''
120 def __init__(
121 self,
122 name='PeterBird2003',
123 data_dir=None,
124 raw_data_url=('https://mirror.pyrocko.org/peterbird.name/oldFTP/PB2002/%s')): # noqa
126 if data_dir is None:
127 data_dir = op.join(config.config().tectonics_dir, name)
129 PlatesDataset.__init__(
130 self,
131 name,
132 data_dir=data_dir,
133 citation=self.__citation)
135 self.raw_data_url = raw_data_url
137 self.filenames = [
138 '2001GC000252_readme.txt',
139 'PB2002_boundaries.dig.txt',
140 'PB2002_orogens.dig.txt',
141 'PB2002_plates.dig.txt',
142 'PB2002_poles.dat.txt',
143 'PB2002_steps.dat.txt']
145 self._full_names = None
147 def full_name(self, name):
148 if not self._full_names:
149 fn = util.data_file(op.join('tectonics', 'bird2003_plates.txt'))
151 with open(fn, 'r') as f:
152 self._full_names = dict(
153 line.strip().split(None, 1) for line in f)
155 return self._full_names[name]
157 def download_if_needed(self):
158 for fn in self.filenames:
159 fpath = self.fpath(fn)
160 if not op.exists(fpath):
161 self.download_file(
162 self.raw_data_url % fn, fpath,
163 status_callback=get_download_callback(
164 'Downloading Bird 2003 plate database...'))
166 def get_boundaries(self):
167 self.download_if_needed()
168 fpath = self.fpath('PB2002_steps.dat.txt')
170 d = defaultdict(list)
171 ntyp = 0
172 type_to_index = {}
173 index_to_type = []
174 with open(fpath, 'rb') as f:
175 data = []
176 for line in f:
177 lsplit = line.split()
178 name = lsplit[1].lstrip(b':').decode('ascii')
179 plate_name1 = str(name[0:2])
180 plate_name2 = str(name[3:5])
181 kind = name[2]
183 alon, alat, blon, blat = list(map(float, lsplit[2:6]))
184 mlat = (alat + blat) * 0.5
185 dlon = ((blon - alon) + 180.) % 360. - 180.
186 mlon = alon + dlon * 0.5
187 typ = str(lsplit[14].strip(b':*').decode('ascii'))
189 if typ not in type_to_index:
190 ntyp += 1
191 type_to_index[typ] = ntyp - 1
192 index_to_type.append(typ)
194 ityp = type_to_index[typ]
195 d[plate_name1, kind, plate_name2].append((mlat, mlon, ityp))
197 d2 = {}
198 for k in d:
199 d2[k] = (
200 num.array([x[:2] for x in d[k]], dtype=float),
201 num.array([x[2] for x in d[k]], dtype=int))
203 fpath = self.fpath('PB2002_boundaries.dig.txt')
204 boundaries = []
205 plate_name1 = ''
206 plate_name2 = ''
207 kind = '-'
208 with open(fpath, 'rb') as f:
209 data = []
210 for line in f:
211 if line.startswith(b'***'):
212 cpoints, itypes = d2[plate_name1, kind, plate_name2]
213 boundaries.append(Boundary(
214 plate_name1=plate_name1,
215 plate_name2=plate_name2,
216 kind=kind,
217 points=num.array(data, dtype=float),
218 cpoints=cpoints,
219 itypes=itypes))
221 boundaries[-1]._type_to_index = type_to_index
222 boundaries[-1]._index_to_type = index_to_type
224 data = []
225 elif line.startswith(b' '):
226 data.append(list(map(float, line.split(b',')))[::-1])
227 else:
228 s = line.strip().decode('ascii')
229 plate_name1 = str(s[0:2])
230 plate_name2 = str(s[3:5])
231 kind = s[2]
233 return boundaries
235 def get_plates(self):
236 self.download_if_needed()
237 fpath = self.fpath('PB2002_plates.dig.txt')
238 plates = []
239 name = ''
240 with open(fpath, 'rb') as f:
241 data = []
242 for line in f:
243 if line.startswith(b'***'):
244 plates.append(Plate(
245 name=name,
246 points=num.array(data, dtype=float)))
248 data = []
249 elif line.startswith(b' '):
250 data.append(list(map(float, line.split(b',')))[::-1])
251 else:
252 name = str(line.strip().decode('ascii'))
254 return plates
257class StrainRateDataset(Dataset):
258 pass
261class GSRM1(StrainRateDataset):
262 '''
263 Global Strain Rate Map. An integrated global model of present-day
264 plate motions and plate boundary deformation'''
265 __citation = '''Kreemer, C., W.E. Holt, and A.J. Haines, "An integrated
266 global model of present-day plate motions and plate boundary deformation",
267 Geophys. J. Int., 154, 8-34, 2003.
268 '''
270 def __init__(
271 self,
272 name='GSRM1.2',
273 data_dir=None,
274 raw_data_url=('https://mirror.pyrocko.org/gsrm.unavco.org/model'
275 '/files/1.2/%s')):
277 if data_dir is None:
278 data_dir = op.join(config.config().tectonics_dir, name)
280 StrainRateDataset.__init__(
281 self,
282 name,
283 data_dir=data_dir,
284 citation=self.__citation)
286 self.raw_data_url = raw_data_url
287 self._full_names = None
288 self._names = None
290 def full_names(self):
291 if not self._full_names:
292 fn = util.data_file(op.join('tectonics', 'gsrm1_plates.txt'))
294 with open(fn, 'r') as f:
295 self._full_names = dict(
296 line.strip().split(None, 1) for line in f)
298 return self._full_names
300 def full_name(self, name):
301 name = self.plate_alt_names().get(name, name)
302 return self.full_names()[name]
304 def plate_names(self):
305 if self._names is None:
306 self._names = sorted(self.full_names().keys())
308 return self._names
310 def plate_alt_names(self):
311 # 'African Plate' is named 'Nubian Plate' in GSRM1
312 return {'AF': 'NU'}
314 def download_if_needed(self, fn):
315 fpath = self.fpath(fn)
316 if not op.exists(fpath):
317 self.download_file(
318 self.raw_data_url % fn, fpath,
319 status_callback=get_download_callback(
320 'Downloading global strain rate map...'))
322 def get_velocities(self, reference_name=None, region=None):
323 reference_name = self.plate_alt_names().get(
324 reference_name, reference_name)
326 if reference_name is None:
327 reference_name = 'NNR'
329 fn = 'velocity_%s.dat' % reference_name
331 self.download_if_needed(fn)
332 fpath = self.fpath(fn)
333 data = []
334 with open(fpath, 'rb') as f:
335 for line in f:
336 if line.strip().startswith(b'#'):
337 continue
338 t = line.split()
339 data.append(list(map(float, t)))
341 arr = num.array(data, dtype=float)
343 if region is not None:
344 points = arr[:, 1::-1]
345 mask = od.points_in_region(points, region)
346 arr = arr[mask, :]
348 lons, lats, veast, vnorth, veast_err, vnorth_err, corr = arr.T
349 return lats, lons, vnorth, veast, vnorth_err, veast_err, corr
352__all__ = '''
353GSRM1
354PeterBird2003
355Plate'''.split()