Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/tectonics.py: 99%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Strain rate (`GSRM1 <https://gsrm.unavco.org/>`_) and plate boundaries 

8(`PeterBird2003 

9<http://peterbird.name/publications/2003_PB2002/2003_PB2002.htm>`_). 

10''' 

11 

12import os.path as op 

13import math 

14from collections import defaultdict 

15 

16import numpy as num 

17 

18from pyrocko import util, config 

19from pyrocko import orthodrome as od 

20from pyrocko.guts_array import Array 

21from pyrocko.guts import Object, String 

22from .util import get_download_callback 

23 

24 

25PI = math.pi 

26 

27 

28class Plate(Object): 

29 

30 name = String.T( 

31 help='Name of the tectonic plate.') 

32 points = Array.T( 

33 dtype=float, shape=(None, 2), 

34 help='Points on the plate.') 

35 

36 def max_interpoint_distance(self): 

37 p = od.latlon_to_xyz(self.points) 

38 return math.sqrt(num.max(num.sum( 

39 (p[num.newaxis, :, :] - p[:, num.newaxis, :])**2, axis=2))) 

40 

41 def contains_point(self, point): 

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

43 

44 def contains_points(self, points): 

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

46 

47 

48class Boundary(Object): 

49 

50 plate_name1 = String.T() 

51 plate_name2 = String.T() 

52 kind = String.T() 

53 points = Array.T(dtype=float, shape=(None, 2)) 

54 cpoints = Array.T(dtype=float, shape=(None, 2)) 

55 itypes = Array.T(dtype=int, shape=(None)) 

56 

57 def split_types(self, groups=None): 

58 xyz = od.latlon_to_xyz(self.points) 

59 xyzmid = (xyz[1:] + xyz[:-1, :]) * 0.5 

60 cxyz = od.latlon_to_xyz(self.cpoints) 

61 d = od.distances3d(xyzmid[num.newaxis, :, :], cxyz[:, num.newaxis, :]) 

62 idmin = num.argmin(d, axis=0) 

63 itypes = self.itypes[idmin] 

64 

65 if groups is None: 

66 groupmap = num.arange(len(self._index_to_type)) 

67 else: 

68 d = {} 

69 for igroup, group in enumerate(groups): 

70 for name in group: 

71 d[name] = igroup 

72 

73 groupmap = num.array( 

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

75 dtype=int) 

76 

77 iswitch = num.concatenate( 

78 ([0], 

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

80 [itypes.size])) 

81 

82 results = [] 

83 for ii in range(iswitch.size-1): 

84 if groups is not None: 

85 tt = [self._index_to_type[ityp] for ityp in num.unique( 

86 itypes[iswitch[ii]:iswitch[ii+1]])] 

87 else: 

88 tt = self._index_to_type[itypes[iswitch[ii]]] 

89 

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

91 

92 return results 

93 

94 

95class TectonicsDataset(object): 

96 ''' 

97 Base class for datasets in :py:mod:`pyrocko.dataset.tectonics`. 

98 ''' 

99 

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

101 self.name = name 

102 self.citation = citation 

103 self.data_dir = data_dir 

104 

105 def fpath(self, filename): 

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

107 

108 def download_file( 

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

110 status_callback=None): 

111 

112 util.download_file( 

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

114 

115 

116class PlatesDataset(TectonicsDataset): 

117 ''' 

118 Base class for datasets describing plate boundaries. 

119 ''' 

120 

121 

122class PeterBird2003(PlatesDataset): 

123 ''' 

124 An updated digital model of plate boundaries. 

125 ''' 

126 

127 __citation = \ 

128 'Bird, Peter. "An updated digital model of plate boundaries." ' \ 

129 'Geochemistry, Geophysics, Geosystems 4.3 (2003).' 

130 

131 def __init__( 

132 self, 

133 name='PeterBird2003', 

134 data_dir=None, 

135 raw_data_url=('https://mirror.pyrocko.org/peterbird.name/oldFTP/PB2002/%s')): # noqa 

136 

137 if data_dir is None: 

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

139 

140 PlatesDataset.__init__( 

141 self, 

142 name, 

143 data_dir=data_dir, 

144 citation=self.__citation) 

145 

146 self.raw_data_url = raw_data_url 

147 

148 self.filenames = [ 

149 '2001GC000252_readme.txt', 

150 'PB2002_boundaries.dig.txt', 

151 'PB2002_orogens.dig.txt', 

152 'PB2002_plates.dig.txt', 

153 'PB2002_poles.dat.txt', 

154 'PB2002_steps.dat.txt'] 

155 

156 self._full_names = None 

157 

158 def full_name(self, name): 

159 if not self._full_names: 

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

161 

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

163 self._full_names = dict( 

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

165 

166 return self._full_names[name] 

167 

168 def download_if_needed(self): 

169 for fn in self.filenames: 

170 fpath = self.fpath(fn) 

171 if not op.exists(fpath): 

172 self.download_file( 

173 self.raw_data_url % fn, fpath, 

174 status_callback=get_download_callback( 

175 'Downloading Bird 2003 plate database...')) 

176 

177 def get_boundaries(self): 

178 self.download_if_needed() 

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

180 

181 d = defaultdict(list) 

182 ntyp = 0 

183 type_to_index = {} 

184 index_to_type = [] 

185 with open(fpath, 'rb') as f: 

186 data = [] 

187 for line in f: 

188 lsplit = line.split() 

189 name = lsplit[1].lstrip(b':').decode('ascii') 

190 plate_name1 = str(name[0:2]) 

191 plate_name2 = str(name[3:5]) 

192 kind = name[2] 

193 

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

195 mlat = (alat + blat) * 0.5 

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

197 mlon = alon + dlon * 0.5 

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

199 

200 if typ not in type_to_index: 

201 ntyp += 1 

202 type_to_index[typ] = ntyp - 1 

203 index_to_type.append(typ) 

204 

205 ityp = type_to_index[typ] 

206 d[plate_name1, kind, plate_name2].append((mlat, mlon, ityp)) 

207 

208 d2 = {} 

209 for k in d: 

210 d2[k] = ( 

211 num.array([x[:2] for x in d[k]], dtype=float), 

212 num.array([x[2] for x in d[k]], dtype=int)) 

213 

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

215 boundaries = [] 

216 plate_name1 = '' 

217 plate_name2 = '' 

218 kind = '-' 

219 with open(fpath, 'rb') as f: 

220 data = [] 

221 for line in f: 

222 if line.startswith(b'***'): 

223 cpoints, itypes = d2[plate_name1, kind, plate_name2] 

224 boundaries.append(Boundary( 

225 plate_name1=plate_name1, 

226 plate_name2=plate_name2, 

227 kind=kind, 

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

229 cpoints=cpoints, 

230 itypes=itypes)) 

231 

232 boundaries[-1]._type_to_index = type_to_index 

233 boundaries[-1]._index_to_type = index_to_type 

234 

235 data = [] 

236 elif line.startswith(b' '): 

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

238 else: 

239 s = line.strip().decode('ascii') 

240 plate_name1 = str(s[0:2]) 

241 plate_name2 = str(s[3:5]) 

242 kind = s[2] 

243 

244 return boundaries 

245 

246 def get_plates(self): 

247 self.download_if_needed() 

248 fpath = self.fpath('PB2002_plates.dig.txt') 

249 plates = [] 

250 name = '' 

251 with open(fpath, 'rb') as f: 

252 data = [] 

253 for line in f: 

254 if line.startswith(b'***'): 

255 plates.append(Plate( 

256 name=name, 

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

258 

259 data = [] 

260 elif line.startswith(b' '): 

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

262 else: 

263 name = str(line.strip().decode('ascii')) 

264 

265 return plates 

266 

267 

268class StrainRateDataset(TectonicsDataset): 

269 ''' 

270 Base class for strain rate datasets. 

271 ''' 

272 

273 

274class GSRM1(StrainRateDataset): 

275 ''' 

276 Global Strain Rate Map. An integrated global model of present-day 

277 plate motions and plate boundary deformation. 

278 ''' 

279 

280 __citation = \ 

281 'Kreemer, C., W.E. Holt, and A.J. Haines, "An integrated global ' \ 

282 'model of present-day plate motions and plate boundary ' \ 

283 'deformation", Geophys. J. Int., 154, 8-34, 2003.' 

284 

285 def __init__( 

286 self, 

287 name='GSRM1.2', 

288 data_dir=None, 

289 raw_data_url=('https://mirror.pyrocko.org/gsrm.unavco.org/model' 

290 '/files/1.2/%s')): 

291 

292 if data_dir is None: 

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

294 

295 StrainRateDataset.__init__( 

296 self, 

297 name, 

298 data_dir=data_dir, 

299 citation=self.__citation) 

300 

301 self.raw_data_url = raw_data_url 

302 self._full_names = None 

303 self._names = None 

304 

305 def full_names(self): 

306 if not self._full_names: 

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

308 

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

310 self._full_names = dict( 

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

312 

313 return self._full_names 

314 

315 def full_name(self, name): 

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

317 return self.full_names()[name] 

318 

319 def plate_names(self): 

320 if self._names is None: 

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

322 

323 return self._names 

324 

325 def plate_alt_names(self): 

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

327 return {'AF': 'NU'} 

328 

329 def download_if_needed(self, fn): 

330 fpath = self.fpath(fn) 

331 if not op.exists(fpath): 

332 self.download_file( 

333 self.raw_data_url % fn, fpath, 

334 status_callback=get_download_callback( 

335 'Downloading global strain rate map...')) 

336 

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

338 reference_name = self.plate_alt_names().get( 

339 reference_name, reference_name) 

340 

341 if reference_name is None: 

342 reference_name = 'NNR' 

343 

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

345 

346 self.download_if_needed(fn) 

347 fpath = self.fpath(fn) 

348 data = [] 

349 with open(fpath, 'rb') as f: 

350 for line in f: 

351 if line.strip().startswith(b'#'): 

352 continue 

353 t = line.split() 

354 data.append(list(map(float, t))) 

355 

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

357 

358 if region is not None: 

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

360 mask = od.points_in_region(points, region) 

361 arr = arr[mask, :] 

362 

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

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

365 

366 

367for cls in PeterBird2003, GSRM1: 

368 cls.__doc__ += '\n\n .. note::\n Please cite: %s\n' \ 

369 % getattr(cls, '_%s__citation' % cls.__name__) 

370 

371 

372__all__ = ''' 

373GSRM1 

374PeterBird2003 

375Plate'''.split()