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

# http://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

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

from __future__ import absolute_import 

from builtins import range 

from builtins import map 

import os.path as op 

import math 

from collections import defaultdict 

import numpy as num 

from pyrocko import util, config 

from pyrocko import orthodrome as od 

from pyrocko.guts_array import Array 

from pyrocko.guts import Object, String 

from .util import get_download_callback 

PI = math.pi 

class Plate(Object): 

name = String.T( 

help='Name of the tectonic plate.') 

points = Array.T( 

dtype=num.float, shape=(None, 2), 

help='Points on the plate.') 

def max_interpoint_distance(self): 

p = od.latlon_to_xyz(self.points) 

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

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

def contains_point(self, point): 

return od.contains_point(self.points, point) 

def contains_points(self, points): 

return od.contains_points(self.points, points) 

class Boundary(Object): 

name1 = String.T() 

name2 = String.T() 

kind = String.T() 

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

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

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

def split_types(self, groups=None): 

xyz = od.latlon_to_xyz(self.points) 

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

cxyz = od.latlon_to_xyz(self.cpoints) 

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

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

itypes = self.itypes[idmin] 

if groups is None: 

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

else: 

d = {} 

for igroup, group in enumerate(groups): 

for name in group: 

d[name] = igroup 

groupmap = num.array( 

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

dtype=num.int) 

iswitch = num.concatenate( 

([0], 

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

[itypes.size])) 

results = [] 

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

if groups is not None: 

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

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

else: 

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

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

return results 

class Dataset(object): 

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

self.name = name 

self.citation = citation 

self.data_dir = data_dir 

def fpath(self, filename): 

return op.join(self.data_dir, filename) 

def download_file( 

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

status_callback=None): 

util.download_file( 

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

class PlatesDataset(Dataset): 

pass 

class PeterBird2003(PlatesDataset): 

'''An updated digital model of plate boundaries.''' 

__citation = ''' 

Bird, Peter. "An updated digital model of plate boundaries." Geochemistry, 

Geophysics, Geosystems 4.3 (2003). 

''' 

def __init__( 

self, 

name='PeterBird2003', 

data_dir=None, 

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

if data_dir is None: 

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

PlatesDataset.__init__( 

self, 

name, 

data_dir=data_dir, 

citation=self.__citation) 

self.raw_data_url = raw_data_url 

self.filenames = [ 

'2001GC000252_readme.txt', 

'PB2002_boundaries.dig.txt', 

'PB2002_orogens.dig.txt', 

'PB2002_plates.dig.txt', 

'PB2002_poles.dat.txt', 

'PB2002_steps.dat.txt'] 

self._full_names = None 

def full_name(self, name): 

if not self._full_names: 

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

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

self._full_names = dict( 

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

return self._full_names[name] 

def download_if_needed(self): 

for fn in self.filenames: 

fpath = self.fpath(fn) 

if not op.exists(fpath): 

self.download_file( 

self.raw_data_url % fn, fpath, 

status_callback=get_download_callback( 

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

def get_boundaries(self): 

self.download_if_needed() 

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

d = defaultdict(list) 

ntyp = 0 

type_to_index = {} 

index_to_type = [] 

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

data = [] 

for line in f: 

t = line.split() 

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

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

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

kind = s[2] 

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

mlat = (alat + blat) * 0.5 

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

mlon = alon + dlon * 0.5 

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

if typ not in type_to_index: 

ntyp += 1 

type_to_index[typ] = ntyp - 1 

index_to_type.append(typ) 

ityp = type_to_index[typ] 

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

d2 = {} 

for k in d: 

d2[k] = ( 

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

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

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

boundaries = [] 

name1 = '' 

name2 = '' 

kind = '-' 

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

data = [] 

for line in f: 

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

cpoints, itypes = d2[name1, kind, name2] 

boundaries.append(Boundary( 

name1=name1, 

name2=name2, 

kind=kind, 

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

cpoints=cpoints, 

itypes=itypes)) 

boundaries[-1]._type_to_index = type_to_index 

boundaries[-1]._index_to_type = index_to_type 

data = [] 

elif line.startswith(b' '): 

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

else: 

s = line.strip() 

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

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

kind = s[2] 

return boundaries 

def get_plates(self): 

self.download_if_needed() 

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

plates = [] 

name = '' 

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

data = [] 

for line in f: 

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

plates.append(Plate( 

name=name, 

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

data = [] 

elif line.startswith(b' '): 

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

else: 

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

return plates 

class StrainRateDataset(Dataset): 

pass 

class GSRM1(StrainRateDataset): 

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

plate motions and plate boundary deformation''' 

__citation = '''Kreemer, C., W.E. Holt, and A.J. Haines, "An integrated 

global model of present-day plate motions and plate boundary deformation", 

Geophys. J. Int., 154, 8-34, 2003.''' 

def __init__( 

self, 

name='GSRM1.2', 

data_dir=None, 

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

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

if data_dir is None: 

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

StrainRateDataset.__init__( 

self, 

name, 

data_dir=data_dir, 

citation=self.__citation) 

self.raw_data_url = raw_data_url 

self._full_names = None 

self._names = None 

def full_names(self): 

if not self._full_names: 

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

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

self._full_names = dict( 

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

return self._full_names 

def full_name(self, name): 

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

return self.full_names()[name] 

def plate_names(self): 

if self._names is None: 

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

return self._names 

def plate_alt_names(self): 

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

return {'AF': 'NU'} 

def download_if_needed(self, fn): 

fpath = self.fpath(fn) 

if not op.exists(fpath): 

self.download_file( 

self.raw_data_url % fn, fpath, 

status_callback=get_download_callback( 

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

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

reference_name = self.plate_alt_names().get( 

reference_name, reference_name) 

if reference_name is None: 

reference_name = 'NNR' 

fn = 'velocity_%s.dat' % reference_name 

self.download_if_needed(fn) 

fpath = self.fpath(fn) 

data = [] 

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

for line in f: 

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

continue 

t = line.split() 

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

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

if region is not None: 

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

mask = od.points_in_region(points, region) 

arr = arr[mask, :] 

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

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

__all__ = ''' 

GSRM1 

PeterBird2003 

Plate'''.split()