1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import math
8import numpy as num
9import scipy.signal
11from pyrocko.orthodrome import positive_region
14class OutOfBounds(Exception):
15 pass
18class Tile(object):
20 def __init__(self, xmin, ymin, dx, dy, data):
21 self.xmin = float(xmin)
22 self.ymin = float(ymin)
23 self.dx = float(dx)
24 self.dy = float(dy)
25 self.data = data
26 self._set_maxes()
28 def _set_maxes(self):
29 self.ny, self.nx = self.data.shape
30 self.xmax = self.xmin + (self.nx-1) * self.dx
31 self.ymax = self.ymin + (self.ny-1) * self.dy
33 def x(self):
34 return self.xmin + num.arange(self.nx) * self.dx
36 def y(self):
37 return self.ymin + num.arange(self.ny) * self.dy
39 def decimate(self, ndeci):
40 assert ndeci % 2 == 0
41 kernel = num.ones((ndeci+1, ndeci+1))
42 kernel /= num.sum(kernel)
43 data = scipy.signal.convolve2d(
44 self.data.astype(float), kernel, mode='valid')
46 self.data = data[::ndeci, ::ndeci].astype(self.data.dtype)
47 self.xmin += ndeci/2
48 self.ymin += ndeci/2
49 self.dx *= ndeci
50 self.dy *= ndeci
51 self._set_maxes()
53 def yextend_with_repeat(self, ymin, ymax):
54 assert ymax >= self.ymax
55 assert ymin <= self.ymin
57 nlo = int(round((self.ymin - ymin) / self.dy))
58 nhi = int(round((ymax - self.ymax) / self.dy))
60 nx, ny = self.nx, self.ny
61 data = num.zeros((ny+nlo+nhi, nx), dtype=self.data.dtype)
62 data[:nlo, :] = self.data[nlo, :]
63 data[nlo:nlo+ny, :] = self.data
64 data[nlo+ny:, :] = self.data[-1, :]
66 self.ymin = ymin
67 self.data = data
68 self._set_maxes()
70 def get(self, x, y):
71 ix = int(round((x - self.xmin) / self.dx))
72 iy = int(round((y - self.ymin) / self.dy))
73 if 0 <= ix < self.nx and 0 <= iy < self.ny:
74 return self.data[iy, ix]
75 else:
76 raise OutOfBounds()
79def multiple_of(x, dx, eps=1e-5):
80 return abs(int(round(x / dx))*dx - x) < dx * eps
83def combine(tiles, region=None):
84 if not tiles:
85 return None
87 dx = tiles[0].dx
88 dy = tiles[0].dy
89 dtype = tiles[0].data.dtype
91 assert all(t.dx == dx for t in tiles)
92 assert all(t.dy == dy for t in tiles)
93 assert all(t.data.dtype == dtype for t in tiles)
94 assert all(multiple_of(t.xmin, dx) for t in tiles)
95 assert all(multiple_of(t.ymin, dy) for t in tiles)
97 if region is None:
98 xmin = min(t.xmin for t in tiles)
99 xmax = max(t.xmax for t in tiles)
100 ymin = min(t.ymin for t in tiles)
101 ymax = max(t.ymax for t in tiles)
102 else:
103 xmin, xmax, ymin, ymax = positive_region(region)
105 if not multiple_of(xmin, dx):
106 xmin = math.floor(xmin / dx) * dx
107 if not multiple_of(xmax, dx):
108 xmax = math.ceil(xmax / dx) * dx
109 if not multiple_of(ymin, dy):
110 ymin = math.floor(ymin / dy) * dy
111 if not multiple_of(ymax, dy):
112 ymax = math.ceil(ymax / dy) * dy
114 nx = int(round((xmax - xmin) / dx)) + 1
115 ny = int(round((ymax - ymin) / dy)) + 1
117 data = num.zeros((ny, nx), dtype=dtype)
118 data[:, :] = 0
120 for t in tiles:
121 for txmin in (t.xmin, t.xmin + 360.):
122 ix = int(round((txmin - xmin) / dx))
123 iy = int(round((t.ymin - ymin) / dy))
124 ixlo = max(ix, 0)
125 ixhi = min(ix+t.nx, nx)
126 iylo = max(iy, 0)
127 iyhi = min(iy+t.ny, ny)
128 jxlo = ixlo-ix
129 jxhi = jxlo + max(0, ixhi - ixlo)
130 jylo = iylo-iy
131 jyhi = jylo + max(0, iyhi - iylo)
132 if iyhi > iylo and ixhi > ixlo:
133 data[iylo:iyhi, ixlo:ixhi] = t.data[jylo:jyhi, jxlo:jxhi]
135 if not num.any(num.isfinite(data)):
136 return None
138 return Tile(xmin, ymin, dx, dy, data)