1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5from __future__ import absolute_import 

6 

7import os.path as op 

8import math 

9from collections import defaultdict 

10 

11import numpy as num 

12 

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 

18 

19 

20PI = math.pi 

21 

22 

23class Plate(Object): 

24 

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

30 

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

35 

36 def contains_point(self, point): 

37 return od.contains_point(self.points, point) 

38 

39 def contains_points(self, points): 

40 return od.contains_points(self.points, points) 

41 

42 

43class Boundary(Object): 

44 

45 name1 = String.T() 

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

51 

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] 

59 

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 

67 

68 groupmap = num.array( 

69 [d[name] for name in self._index_to_type], 

70 dtype=int) 

71 

72 iswitch = num.concatenate( 

73 ([0], 

74 num.where(groupmap[itypes[1:]] != groupmap[itypes[:-1]])[0]+1, 

75 [itypes.size])) 

76 

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

84 

85 results.append((tt, self.points[iswitch[ii]:iswitch[ii+1]+1])) 

86 

87 return results 

88 

89 

90class Dataset(object): 

91 

92 def __init__(self, name, data_dir, citation): 

93 self.name = name 

94 self.citation = citation 

95 self.data_dir = data_dir 

96 

97 def fpath(self, filename): 

98 return op.join(self.data_dir, filename) 

99 

100 def download_file( 

101 self, url, fpath, username=None, password=None, 

102 status_callback=None): 

103 

104 util.download_file( 

105 url, fpath, username, password, status_callback=status_callback) 

106 

107 

108class PlatesDataset(Dataset): 

109 pass 

110 

111 

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

120 

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 

126 

127 if data_dir is None: 

128 data_dir = op.join(config.config().tectonics_dir, name) 

129 

130 PlatesDataset.__init__( 

131 self, 

132 name, 

133 data_dir=data_dir, 

134 citation=self.__citation) 

135 

136 self.raw_data_url = raw_data_url 

137 

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

145 

146 self._full_names = None 

147 

148 def full_name(self, name): 

149 if not self._full_names: 

150 fn = util.data_file(op.join('tectonics', 'bird2003_plates.txt')) 

151 

152 with open(fn, 'r') as f: 

153 self._full_names = dict( 

154 line.strip().split(None, 1) for line in f) 

155 

156 return self._full_names[name] 

157 

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

166 

167 def get_boundaries(self): 

168 self.download_if_needed() 

169 fpath = self.fpath('PB2002_steps.dat.txt') 

170 

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 t = line.split() 

179 s = t[1].lstrip(b':') 

180 name1 = str(s[0:2].decode('ascii')) 

181 name2 = str(s[3:5].decode('ascii')) 

182 kind = s[2] 

183 

184 alon, alat, blon, blat = list(map(float, t[2:6])) 

185 mlat = (alat + blat) * 0.5 

186 dlon = ((blon - alon) + 180.) % 360. - 180. 

187 mlon = alon + dlon * 0.5 

188 typ = str(t[14].strip(b':*').decode('ascii')) 

189 

190 if typ not in type_to_index: 

191 ntyp += 1 

192 type_to_index[typ] = ntyp - 1 

193 index_to_type.append(typ) 

194 

195 ityp = type_to_index[typ] 

196 d[name1, kind, name2].append((mlat, mlon, ityp)) 

197 

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

203 

204 fpath = self.fpath('PB2002_boundaries.dig.txt') 

205 boundaries = [] 

206 name1 = '' 

207 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[name1, kind, name2] 

214 boundaries.append(Boundary( 

215 name1=name1, 

216 name2=name2, 

217 kind=kind, 

218 points=num.array(data, dtype=float), 

219 cpoints=cpoints, 

220 itypes=itypes)) 

221 

222 boundaries[-1]._type_to_index = type_to_index 

223 boundaries[-1]._index_to_type = index_to_type 

224 

225 data = [] 

226 elif line.startswith(b' '): 

227 data.append(list(map(float, line.split(b',')))[::-1]) 

228 else: 

229 s = line.strip() 

230 name1 = str(s[0:2].decode('ascii')) 

231 name2 = str(s[3:5].decode('ascii')) 

232 kind = s[2] 

233 

234 return boundaries 

235 

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

248 

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

254 

255 return plates 

256 

257 

258class StrainRateDataset(Dataset): 

259 pass 

260 

261 

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

270 

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

277 

278 if data_dir is None: 

279 data_dir = op.join(config.config().tectonics_dir, name) 

280 

281 StrainRateDataset.__init__( 

282 self, 

283 name, 

284 data_dir=data_dir, 

285 citation=self.__citation) 

286 

287 self.raw_data_url = raw_data_url 

288 self._full_names = None 

289 self._names = None 

290 

291 def full_names(self): 

292 if not self._full_names: 

293 fn = util.data_file(op.join('tectonics', 'gsrm1_plates.txt')) 

294 

295 with open(fn, 'r') as f: 

296 self._full_names = dict( 

297 line.strip().split(None, 1) for line in f) 

298 

299 return self._full_names 

300 

301 def full_name(self, name): 

302 name = self.plate_alt_names().get(name, name) 

303 return self.full_names()[name] 

304 

305 def plate_names(self): 

306 if self._names is None: 

307 self._names = sorted(self.full_names().keys()) 

308 

309 return self._names 

310 

311 def plate_alt_names(self): 

312 # 'African Plate' is named 'Nubian Plate' in GSRM1 

313 return {'AF': 'NU'} 

314 

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

322 

323 def get_velocities(self, reference_name=None, region=None): 

324 reference_name = self.plate_alt_names().get( 

325 reference_name, reference_name) 

326 

327 if reference_name is None: 

328 reference_name = 'NNR' 

329 

330 fn = 'velocity_%s.dat' % reference_name 

331 

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

341 

342 arr = num.array(data, dtype=float) 

343 

344 if region is not None: 

345 points = arr[:, 1::-1] 

346 mask = od.points_in_region(points, region) 

347 arr = arr[mask, :] 

348 

349 lons, lats, veast, vnorth, veast_err, vnorth_err, corr = arr.T 

350 return lats, lons, vnorth, veast, vnorth_err, veast_err, corr 

351 

352 

353__all__ = ''' 

354GSRM1 

355PeterBird2003 

356Plate'''.split()