A friendly earthquake detector.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

206 lines
5.9 KiB

import math
import numpy as num
import logging
from pyrocko.guts import Object, Float
from pyrocko import orthodrome as od
from lassie import receiver
logger = logging.getLogger('lassie.grid')
guts_prefix = 'lassie'
class Grid(Object):
pass
class Carthesian3DGrid(Grid):
lat = Float.T(default=0.0)
lon = Float.T(default=0.0)
xmin = Float.T()
xmax = Float.T()
ymin = Float.T()
ymax = Float.T()
zmin = Float.T()
zmax = Float.T()
dx = Float.T()
dy = Float.T()
dz = Float.T()
def __init__(self, **kwargs):
for k in kwargs.keys():
kwargs[k] = float(kwargs[k])
Grid.__init__(self, **kwargs)
self.update()
def size(self):
nx, ny, nz = self._shape()
return nx * ny * nz
def max_delta(self):
return max(self.dx, self.dy, self.dz)
def update(self):
self._coords = None
def _shape(self):
nx = int(round(((self.xmax - self.xmin) / self.dx))) + 1
ny = int(round(((self.ymax - self.ymin) / self.dy))) + 1
nz = int(round(((self.zmax - self.zmin) / self.dz))) + 1
return nx, ny, nz
def _get_coords(self):
if self._coords is None:
nx, ny, nz = self._shape()
x = num.linspace(self.xmin, self.xmax, nx)
y = num.linspace(self.ymin, self.ymax, ny)
z = num.linspace(self.zmin, self.zmax, nz)
self._coords = x, y, z
return self._coords
def index_to_location(self, i):
nx, ny, nz = self._shape()
iz, iy, ix = num.unravel_index(i, (nz, ny, nx))
x, y, z = self._get_coords()
return self.lat, self.lon, x[ix], y[iy], z[iz]
def distances(self, receivers):
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
rx, ry = receiver.receivers_coords(
receivers, system=('ne', self.lat, self.lon))
rz = num.array([r.z for r in receivers], dtype=num.float)
na = num.newaxis
nr = len(receivers)
distances = num.sqrt(
(x[na, na, :, na] - rx[na, na, na, :])**2 +
(y[na, :, na, na] - ry[na, na, na, :])**2 +
(z[:, na, na, na] - rz[na, na, na, :])**2).reshape((nx*ny*nz, nr))
return distances
def surface_points(self, system='latlon'):
x, y, z = self._get_coords()
xs = num.tile(x, y.size)
ys = num.repeat(y, x.size)
if system == 'latlon':
return od.ne_to_latlon(self.lat, self.lon, xs, ys)
elif system[0] == 'ne':
lat0, lon0 = system[1:]
if lat0 == self.lat and lon0 == self.lon:
return xs, ys
else:
elats, elons = od.ne_to_latlon(self.lat, self.lon, xs, ys)
return od.latlon_to_ne_numpy(lat0, lon0, elats, elons)
def plot_points(self, axes, system='latlon'):
x, y = self.surface_points(system=system)
axes.plot(y, x, '.', color='black', ms=1.0)
def plot(
self, axes, a,
amin=None,
amax=None,
z_slice=None,
cmap=None,
projection='2d',
system='latlon'):
if system == 'latlon':
assert False, 'not implemented yet'
if not (system[1] == self.lat and system[2] == self.lon):
assert False, 'not implemented yet'
nx, ny, nz = self._shape()
x, y, z = self._get_coords()
a3d = a.reshape((nz, ny, nx))
try:
import vtk
use_vtk = True
except ImportError as e:
logger.warn(e)
use_vtk = False
if projection == '2d' or use_vtk:
if z_slice is not None:
iz = num.argmin(num.abs(z-z_slice))
a2d = a3d[iz, :, :]
else:
a2d = num.max(a3d, axis=0)
axes.pcolormesh(y, x, a2d.T, vmin=amin, vmax=amax, cmap=cmap)
if projection == '3d':
#if use_vtk:
if True:
from lassie import vtk_graph
actors = vtk_graph.grid_actors(a3d.T, x, y, z, normalize=True)
vtk_graph.render_actors(actors)
else:
iy = len(y)
ix = len(x)
iz = len(z)
x = num.repeat(x, iy*iz)
y = num.repeat(num.tile(y, ix), iz)
z = num.tile(z, ix*iy)
a3d /= num.max(a3d)
#a3d = num.exp(a3d)
iplot = num.where(a3d.T!=num.nan)#>num.median(a3d))
iflat = num.ravel(iplot)
axes.scatter(x[iflat], y[iflat], z[iflat], s=a3d.T[iplot]*20.)
#, c=a3d[iflat],
#vmin=num.min(a3d),
#vmax=num.max(a3d),
#cmap=cmap)
axes.invert_zaxis()
axes.set_xlabel('E [m]')
axes.set_ylabel('N [m]')
axes.set_ylabel('D [m]')
def geometrical_normalization(grid, receivers):
distances = grid.distances(receivers)
delta_grid = grid.max_delta()
delta_ring = delta_grid * 3.0
ngridpoints, nstations = distances.shape
norm_map = num.zeros(ngridpoints)
for istation in xrange(nstations):
dists_station = distances[:, istation]
dist_min = num.floor(num.min(dists_station) / delta_ring) * delta_ring
dist_max = num.ceil(num.max(dists_station) / delta_ring) * delta_ring
dist = dist_min
while dist < dist_max:
indices = num.where(num.logical_and(
dist <= dists_station,
dists_station < dist + delta_ring))[0]
nexpect = math.pi * ((dist + delta_ring)**2 - dist**2) /\
delta_grid**2
norm_map[indices] += indices.size / nexpect
dist += delta_ring
norm_map /= nstations
return norm_map
__all__ = [
'Grid',
'Carthesian3DGrid',
]