Browse Source

faults: added fault database

pull/4/head
mmetz 10 months ago
parent
commit
d22faf275c
  1. 1
      .gitattributes
  2. 4
      setup.py
  3. 13
      src/dataset/fault_loader/__init__.py
  4. 160
      src/dataset/fault_loader/active_faults.py
  5. 278
      src/dataset/fault_loader/diss.py
  6. 188
      src/dataset/fault_loader/edsf.py
  7. 73
      src/dataset/faults.py
  8. 24
      src/util.py

1
.gitattributes

@ -0,0 +1 @@
* text=auto

4
setup.py

@ -44,6 +44,7 @@ packages = [
'%s' % packname,
'%s.apps' % packname,
'%s.dataset' % packname,
'%s.dataset.fault_loader' % packname,
'%s.gm' % packname,
'%s.io' % packname,
'%s.ls' % packname,
@ -91,6 +92,7 @@ setup(
'Topic :: Scientific/Engineering'],
package_data={
packname: ['data/tsumaps_all_params.*',
'data/examples/*'
'data/examples/*',
'data/faults/*'
] + get_readme_paths()}
)

13
src/dataset/fault_loader/__init__.py

@ -0,0 +1,13 @@
#/usr/bin/python3
try:
from pyrocko.dataset import active_faults
except ImportError:
from . import active_faults
class SeismogenicSource(active_faults.ActiveFault):
def km_to_m(self):
pass
def to_fm(self):
return None

160
src/dataset/fault_loader/active_faults.py

@ -0,0 +1,160 @@
# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import absolute_import
import re
import logging
import numpy as num
from os import path as op
from collections import OrderedDict
import json
from pyrocko import util as putil
from ewrica import util
def parse_3tup(s):
m = re.match(r'^\(?([^,]+),([^,]*),([^,]*)\)$', s)
if m:
return [float(m.group(1)) if m.group(1) else None for i in range(3)]
else:
return [None, None, None]
logger = logging.getLogger('ewrica.dataset.fault_loader.activefaults')
def active_faults_fn():
name = 'gem_active_faults.geojson'
return op.join(util.HiddenDataDir.faults_path(), '%s' % name)
class Fault(object):
__fields__ = OrderedDict()
__slots__ = list(__fields__.keys())
def get_property(self, fault_obj, attr):
try:
values = float(fault_obj['properties'][attr][1:3])
except KeyError:
if attr == 'lower_seis_depth' or attr == 'upper_seis_depth':
values = 0
else:
values = -999
return values
def __init__(self, f):
nodes = f['geometry']['coordinates']
props = f['properties']
lons = [p[0] for p in nodes]
lats = [p[1] for p in nodes]
self.lon = lons
self.lat = lats
self.slip_type = props.get('slip_type', 'Unknown')
for attr, attr_type in self.__fields__.items():
if attr in props:
if props[attr] is None or props[attr] is '':
continue
if isinstance(props[attr], attr_type):
v = props[attr]
elif attr_type is float:
try:
v = parse_3tup(props[attr])[0]
except TypeError:
v = float(props[attr])
elif attr_type is int:
v = int(props[attr])
else:
v = None
setattr(self, attr, v)
def get_surface_line(self):
arr = num.empty((len(self.lat), 2))
for i in range(len(self.lat)):
arr[i, 0] = self.lat[i]
arr[i, 1] = self.lon[i]
return arr
def __str__(self):
d = {attr: getattr(self, attr) for attr in self.__fields__.keys()}
return '\n'.join(['%s: %s' % (attr, val) for attr, val in d.items()])
class ActiveFault(Fault):
__fields__ = OrderedDict([
('lat', list),
('lon', list),
('average_dip', float),
('average_rake', float),
('lower_seis_depth', float),
('upper_seis_depth', float),
('slip_type', str),
])
class ActiveFaults(object):
'''
Get GEM active faults database
Styron, Richard, and Marco Pagani.
The GEM Global Active Faults Database. Earthquake Spectra, vol. 36,
no. 1_suppl, Oct. 2020, pp. 160180, doi:10.1177/8755293020944182.
code copied from pyrocko python package:
Heimann, Sebastian; Kriegerowski, Marius; Isken, Marius; Cesca, Simone;
Daout, Simon; Grigoli, Francesco; Juretzek, Carina; Megies, Tobias;
Nooshiri, Nima; Steinberg, Andreas; Sudhaus, Henriette;
Vasyura-Bathke, Hannes; Willey, Timothy; Dahm, Torsten (2017):
Pyrocko - An open-source seismology toolbox and library.
V. 0.3. GFZ Data Services. http://doi.org/10.5880/GFZ.2.1.2017.001
'''
URL_GEM_ACTIVE_FAULTS = 'https://raw.githubusercontent.com/cossatot/gem-global-active-faults/master/geojson/gem_active_faults.geojson' # noqa
def __init__(self):
self.fname_active_faults = active_faults_fn()
if not op.exists(self.fname_active_faults):
self.download()
self.active_faults = []
self._load_faults(self.fname_active_faults, ActiveFault)
def _load_faults(self, fname, cls):
with open(fname, 'r') as f:
gj = json.load(f)
faults = gj['features']
for f in faults:
fault = cls(f)
self.active_faults.append(fault)
logger.debug('loaded %d fault', self.nactive_faults)
def download(self):
logger.info('Downloading GEM active faults database...')
putil.download_file(self.URL_GEM_ACTIVE_FAULTS,
self.fname_active_faults)
@property
def nactive_faults(self):
return len(self.active_faults)
def nactive_faults_nodes(self):
return int(sum(len(f.lat) for f in self.active_faults))
def get_coords(self):
return num.array([f['coordinates'] for f in self.active_faults])
if __name__ == '__main__':
logging.basicConfig(level=logging.DEBUG)
activefaults = ActiveFaults()

278
src/dataset/fault_loader/diss.py

@ -0,0 +1,278 @@
#/usr/bin/python3
import os
import os.path as op
from collections import OrderedDict
import shapefile # PyShp package
import numpy as num
from pyrocko import orthodrome as od, moment_tensor as pmt
from pyrocko.model import Event
from ewrica.util import HiddenDataDir, is_dir
from . import SeismogenicSource
km = 1000.
class DISSError(Exception):
pass
source_to_fn = dict(
individual='ISS{version}',
toptrace_composite='CSSTOP{version}',
planes_composite='CSSPLN{version}',
debated='DSS{version}',
subduction_contour='SUBDCNT{version}',
subduction_zone='SUBDZON{version}',
layer='DISS{version}')
class Contour(SeismogenicSource):
__fields__ = OrderedDict([
('IDSource', str),
('Depth_km', float),
('lat', list),
('lon', list)
])
def km_to_m(self):
self.Depth_km *= km
class SubductionZone(SeismogenicSource):
__fields__ = OrderedDict([
('IDSource', str),
('SourceName', str),
('lat', list),
('lon', list),
('MinDepth', float),
('MaxDepth', float),
('DipDir', str),
('ConvAzMin', float),
('ConvAzMax', float),
('ConvRMin', float),
('ConvRMax', float),
('MaxMag', float),
('Created', str),
('Updated', str),
('MinDepthQ', int),
('MaxDepthQ', int),
('DipDirQ', int),
('ConvAzQ', int),
('ConvRQ', int),
('MaxMagQ', int)
])
def km_to_m(self):
for attr in (self.MinDepth, self.MaxDepth):
attr *= km
class CompositeSeismogenicSource(SeismogenicSource):
__fields__ = OrderedDict([
('IDSource', str),
('SourceName', str),
('lat', list),
('lon', list),
('MinDepth', float),
('MaxDepth', float),
('DipMin', float),
('DipMax', float),
('RakeMin', float),
('RakeMax', float),
('SRmin', float),
('SRmax', float),
('MaxMag', float),
('Created', str),
('Updated', str),
('MinDepthQ', int),
('MaxDepthQ', int),
('StrikeQ', int),
('DipQ', int),
('RakeQ', int),
('SlipRateQ', int),
('MaxMagQ', int)
])
def km_to_m(self):
for attr in (self.MinDepth, self.MaxDepth):
attr *= km
def to_event(self):
dip = num.mean([self.DipMin, self.DipMax])
rake = num.mean([self.RakeMin, self.RakeMax])
mag = self.MaxMag
if len(self.lat) % 2:
idx = int((len(self.lat) - 1) / 2.)
lat, lon = self.lat[idx], self.lon[idx]
else:
idx1, idx2 = int(len(self.lat) / 2. - 1), int(len(self.lat) / 2.)
n, e = od.latlon_to_ne_numpy(self.lat[idx1], self.lon[idx1],
self.lat[idx2], self.lon[idx2])
lat, lon = od.ne_to_latlon(self.lat[idx1], self.lon[idx1],
n / 2., e / 2.)
lat, lon = lat[0], lon[0]
azi, _ = od.azibazi_numpy(
num.array(self.lat)[:-1], num.array(self.lon)[:-1],
num.array(self.lat)[1:], num.array(self.lon)[1:])
strike = num.mean(azi)
mt = pmt.MomentTensor(strike=strike, dip=dip, rake=rake, magnitude=mag)
return Event(
lat=lat, lon=lon, depth=num.mean([self.MinDepth, self.MaxDepth]),
name=self.IDSource, moment_tensor=mt)
class IndividualSeismogenicSource(SeismogenicSource):
__fields__ = OrderedDict([
('IDSource', str),
('SourceName', str),
('Length', float),
('Width', float),
('lat', list),
('lon', list),
('MinDepth', float),
('MaxDepth', float),
('Strike', float),
('Dip', float),
('Rake', float),
('AvgDispl', float),
('SRmin', float),
('SRmax', float),
('RecIntMin', float),
('RecIntMax', float),
('LatestEq', str),
('MaxMag', float),
('LatUL', float),
('LonUL', float),
('LatUR', float),
('LonUR', float),
('LatLR', float),
('LonLR', float),
('LatLL', float),
('LonLL', float),
('Created', str),
('Updated', str),
('LengthQ', int),
('WidthQ', int),
('MinDepthQ', int),
('MaxDepthQ', int),
('StrikeQ', int),
('DipQ', int),
('RakeQ', int),
('AvgDisplQ', int),
('SlipRateQ', int),
('RecIntQ', int),
('MaxMagQ', int),
('LocationQ', int)
])
def km_to_m(self):
for attr in (self.Length, self.Width, self.MinDepth, self.MaxDepth):
attr *= km
def to_event(self):
n, e = od.latlon_to_ne(self.LatUL, self.LonUL,
self.LatLR, self.LonLR)
lat, lon = od.ne_to_latlon(self.LatUL, self.LonUL,
n / 2., e / 2.)
mt = pmt.MomentTensor(strike=self.Strike, dip=self.Dip,
rake=self.Rake, magnitude=self.MaxMag)
return Event(
lat=lat, lon=lon, depth=num.mean([self.MinDepth, self.MaxDepth]),
name=self.IDSource, moment_tensor=mt)
def load_faults_from_shp(version='3.2.1', include='all'):
template = op.join(
HiddenDataDir.faults_path(),
'DISS_{version}_shp',
'SHP')
fault_path = template.format(version=version)
if not is_dir(fault_path):
raise DISSError(
'Could not find DISS shape directory "{fault_path}". If not '
'downloaded, get the DISS database from {url} in SHP format and '
'store the individual files in {path}.'.format(
fault_path=fault_path,
url='http://diss.rm.ingv.it/diss/',
path=template.format(version='<version_number>')))
if type(include) == str:
include = [include]
if 'all' in include:
files = [os.path.splitext(f)[0] for f in os.listdir(fault_path)
if f.endswith('.shp')]
else:
files = [
op.join(
fault_path,
source_to_fn[incl].format(version=version.replace('.', '')))
for incl in include]
faults = []
for f in files:
f = os.path.basename(f)
kwargs = dict([
(ext, open(os.path.join(fault_path, '.'.join([f, ext])), 'rb'))
for ext in ['shp', 'prj', 'dbf', 'shx']])
shape = shapefile.Reader(encoding='ISO-8859-1', **kwargs)
geoj = shape.shapeRecord(0).__geo_interface__
if 'Subduction' in geoj['properties']['LinkToInfo']:
if 'Polygon' == geoj['geometry']['type']:
source_type = SubductionZone
elif 'LineString' == geoj['geometry']['type']:
source_type = Contour
elif 'SASources' in geoj['properties']['LinkToInfo']:
source_type = CompositeSeismogenicSource
elif 'GGSources' in geoj['properties']['LinkToInfo']:
source_type = IndividualSeismogenicSource
else:
source_type = SeismogenicSource
for sh in shape.shapeRecords():
geoj = sh.__geo_interface__
if isinstance(geoj['geometry']['coordinates'][0][0], tuple):
pts = geoj['geometry']['coordinates']
geoj['geometry']['coordinates'] = [
(p[0], p[1])
for p in pts[0]]
fault = source_type(geoj)
fault.km_to_m()
faults.append(fault)
return faults
__all__ = [
'SeismogenicSource',
'Contour',
'IndividualSeismogenicSource',
'CompositeSeismogenicSource',
'SubductionZone',
'load_faults_from_shp']

188
src/dataset/fault_loader/edsf.py

@ -0,0 +1,188 @@
#/usr/bin/python3
import logging
import os
import os.path as op
from collections import OrderedDict
import shutil
import shapefile # PyShp package
import numpy as num
import datetime
from pyrocko import orthodrome as od, moment_tensor as pmt, util as putil
from pyrocko.model import Event
from ewrica.util import HiddenDataDir, is_dir
from . import SeismogenicSource
logger = logging.getLogger('ewrica.dataset.fault_loader.edsf')
km = 1000.
class EDSFError(Exception):
pass
source_to_fn = dict(
crustal='CFS',
subduction='SBD')
folder_template = 'DB-{database}-ESRI-shapefile'
class IndividualSeismogenicSourceTop(SeismogenicSource):
__fields__ = OrderedDict([
('IDSOURCE', str),
('SOURCENAME', str),
('MINDEPTH', float),
('MAXDEPTH', float),
('STRIKEMIN', int),
('STRIKEMAX', int),
('DIPMIN', int),
('DIPMAX', int),
('RAKEMIN', int),
('RAKEMAX', int),
('SRMIN', float),
('SRMAX', float),
('MAXMAG', float),
('REGION', str),
('CREATED', datetime.date),
('UPDATED', datetime.date),
('LINK', str)])
def km_to_m(self):
try:
for attr in (self.MINDEPTH, self.MAXDEPTH):
attr *= km
except AttributeError:
return
def to_event(self):
if len(self.lat) % 2:
idx = int((len(self.lat) - 1) / 2.)
lat, lon = self.lat[idx], self.lon[idx]
else:
idx1, idx2 = int(len(self.lat) / 2. - 1), int(len(self.lat) / 2.)
n, e = od.latlon_to_ne_numpy(self.lat[idx1], self.lon[idx1],
self.lat[idx2], self.lon[idx2])
lat, lon = od.ne_to_latlon(self.lat[idx1], self.lon[idx1],
n / 2., e / 2.)
lat, lon = lat[0], lon[0]
mt = pmt.MomentTensor(
strike=num.mean([self.STRIKEMIN, self.STRIKEMAX]),
dip=num.mean([self.DIPMIN, self.DIPMAX]),
rake=num.mean([self.RAKEMIN, self.RAKEMAX]),
magnitude=self.MAXMAG)
return Event(
lat=lat, lon=lon, depth=num.mean([self.MINDEPTH, self.MAXDEPTH]),
name=self.IDSOURCE, moment_tensor=mt)
class IndividualSeismogenicSourcePlane(IndividualSeismogenicSourceTop):
pass
def load_faults_from_shp(include='all'):
fault_path = op.join(HiddenDataDir.faults_path(), 'EDSF')
if not is_dir(fault_path):
download(fault_path)
if len(os.listdir(fault_path)) == 0:
download(fault_path)
if type(include) == str:
include = [include]
files = []
if 'all' in include:
for k, v in source_to_fn.items():
folder = op.join(fault_path, folder_template.format(database=v))
files += [
op.join(folder, os.path.splitext(f)[0]) for f in os.listdir(
folder) if f.endswith('.shp')]
else:
for k, v in source_to_fn.items():
if k not in include:
continue
folder = op.join(fault_path, folder_template.format(database=v))
files += [
op.join(folder, os.path.splitext(f)[0]) for f in os.listdir(
folder) if f.endswith('.shp')]
faults = []
for f in files:
kwargs = dict([
(ext, open('.'.join([f, ext]), 'rb'))
for ext in ['shp', 'prj', 'dbf', 'shx']])
shape = shapefile.Reader(encoding='ISO-8859-1', **kwargs)
geoj = shape.shapeRecord(0).__geo_interface__
if 'Polygon' == geoj['geometry']['type']:
source_type = IndividualSeismogenicSourcePlane
elif 'LineString' == geoj['geometry']['type']:
source_type = IndividualSeismogenicSourceTop
else:
source_type = SeismogenicSource
for sh in shape.shapeRecords():
geoj = sh.__geo_interface__
if isinstance(geoj['geometry']['coordinates'][0][0], tuple):
pts = geoj['geometry']['coordinates']
geoj['geometry']['coordinates'] = [
(p[0], p[1])
for p in pts[0]]
fault = source_type(geoj)
fault.km_to_m()
faults.append(fault)
return faults
def download(fault_path):
import zipfile
logger.info('Downloading EDSF database...')
putil.ensuredirs(fault_path)
url_template = 'http://diss.rm.ingv.it/share-edsf/sharedata/downloads/{}.zip'.format(folder_template) # noqa
fnzip_template = op.join(fault_path, '{}.zip'.format(folder_template))
for db in ('CFS', 'SBD'):
url = url_template.format(database=db)
fnzip = fnzip_template.format(database=db)
folder = folder_template.format(database=db)
temp_path = op.join(fault_path, 'temp')
putil.download_file(url, fnzip)
with zipfile.ZipFile(fnzip, 'r') as f:
f.extractall(temp_path)
shutil.move(op.join(temp_path, folder), op.join(fault_path, folder))
os.remove(fnzip)
__all__ = [
'IndividualSeismogenicSourcePlane',
'IndividualSeismogenicSourceTop',
'load_faults_from_shp']

73
src/dataset/faults.py

@ -0,0 +1,73 @@
#/usr/bin/python3
import os
import numpy as num
from pyrocko import orthodrome as pod
from .fault_loader import diss, edsf, active_faults
def load_diss(*args, **kwargs):
return diss.load_faults_from_shp(*args, **kwargs)
def load_edsf(*args, **kwargs):
return edsf.load_faults_from_shp(*args, **kwargs)
def load_gem(**kwargs):
return active_faults.ActiveFaults().active_faults
def extract_faults(
faults,
wesn=(),
lat=None,
lon=None,
radius=None,
strike=None,
dip=None,
rake=None):
if not wesn:
if lat is not None and lon is not None and radius is not None:
wesn = pod.radius_to_region(lat, lon, radius)
if wesn:
faults = [
f for f in faults
if pod.points_in_region(num.vstack((f.lat, f.lon)).T, wesn).any()]
if strike:
faults = [
f for f in faults
if hasattr(f, 'STRIKEMIN') and hasattr(f, 'STRIKEMAX')
if f.STRIKEMIN % 360. <= strike and f.STRIKEMAX % 360. >= strike]
if dip:
faults = [
f for f in faults
if f.DIPMIN % 90. <= dip and f.DIPMAX % 90. >= dip]
if rake:
faults = [
f for f in faults
if f.RAKEMIN % 360. <= rake and f.RAKEMAX % 360. >= rake]
return faults
def load_faults_from_database(database='GEM', **kwargs):
database = database.upper()
if database == 'DISS':
return load_diss(**kwargs)
elif database == 'EDSF':
return load_edsf(**kwargs)
elif database == 'GEM':
return load_gem(**kwargs)
__all__ = [
'load_diss',
'load_edsf',
'load_gem']

24
src/util.py

@ -7,6 +7,8 @@ import logging.handlers
import os
import os.path as op
from pyrocko import util as putil
def data_file(fn):
return os.path.join(os.path.split(__file__)[0], 'data', fn)
@ -16,6 +18,28 @@ def is_dir(dst):
return op.exists(dst)
ewrica_dir_tmpl = op.expanduser(op.join('~', '.ewrica'))
class HiddenDataDir(object):
@staticmethod
def create():
putil.ensuredirs(ewrica_dir_tmpl)
@staticmethod
def base_path():
HiddenDataDir.create()
return ewrica_dir_tmpl
@staticmethod
def faults_path():
path = op.join(ewrica_dir_tmpl, 'faults')
putil.ensuredirs(op.join(ewrica_dir_tmpl, 'faults'))
return path
def setup_file_logging(logger, fn):
logger.handlers = []
file_handler = logging.handlers.RotatingFileHandler(

Loading…
Cancel
Save