1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import numpy as num
8import scipy.signal
10from pyrocko.orthodrome import positive_region
12km = 1e3
15class OutOfBounds(Exception):
16 pass
19class Tile(object):
21 def __init__(self, xmin, ymin, dx, dy, data):
22 self.xmin = float(xmin)
23 self.ymin = float(ymin)
24 self.dx = float(dx)
25 self.dy = float(dy)
26 self.data = data
27 self._set_maxes()
29 def _set_maxes(self):
30 self.ny, self.nx = self.data.shape
31 self.xmax = self.xmin + (self.nx-1) * self.dx
32 self.ymax = self.ymin + (self.ny-1) * self.dy
34 def x(self):
35 return self.xmin + num.arange(self.nx) * self.dx
37 def y(self):
38 return self.ymin + num.arange(self.ny) * self.dy
40 def decimate(self, ndeci):
41 assert ndeci % 2 == 0
42 kernel = num.ones((ndeci+1, ndeci+1))
43 kernel /= num.sum(kernel)
44 data = scipy.signal.convolve2d(
45 self.data.astype(float), kernel, mode='valid')
47 self.data = data[::ndeci, ::ndeci].astype(self.data.dtype)
48 self.xmin += ndeci/2
49 self.ymin += ndeci/2
50 self.dx *= ndeci
51 self.dy *= ndeci
52 self._set_maxes()
54 def yextend_with_repeat(self, ymin, ymax):
55 assert ymax >= self.ymax
56 assert ymin <= self.ymin
58 nlo = int(round((self.ymin - ymin) / self.dy))
59 nhi = int(round((ymax - self.ymax) / self.dy))
61 nx, ny = self.nx, self.ny
62 data = num.zeros((ny+nlo+nhi, nx), dtype=self.data.dtype)
63 data[:nlo, :] = self.data[nlo, :]
64 data[nlo:nlo+ny, :] = self.data
65 data[nlo+ny:, :] = self.data[-1, :]
67 self.ymin = ymin
68 self.data = data
69 self._set_maxes()
71 def get(self, x, y):
72 ix = int(round((x - self.xmin) / self.dx))
73 iy = int(round((y - self.ymin) / self.dy))
74 if 0 <= ix < self.nx and 0 <= iy < self.ny:
75 return self.data[iy, ix]
76 else:
77 raise OutOfBounds()
79 @classmethod
80 def from_grd(cls, path):
81 '''Load from GMT .grd'''
82 from pyrocko.plot import gmtpy
83 lon, lat, z = gmtpy.loadgrd(path)
84 return cls(lon[0], lat[0], lon[1] - lon[0], lat[1] - lat[0], z)
87def multiple_of(x, dx, eps=1e-5):
88 return abs(int(round(x / dx))*dx - x) < dx * eps
91def combine(tiles, region=None):
92 if not tiles:
93 return None
95 dx = tiles[0].dx
96 dy = tiles[0].dy
97 dtype = tiles[0].data.dtype
99 assert all(t.dx == dx for t in tiles)
100 assert all(t.dy == dy for t in tiles)
101 assert all(t.data.dtype == dtype for t in tiles)
102 assert all(multiple_of(t.xmin, dx) for t in tiles)
103 assert all(multiple_of(t.ymin, dy) for t in tiles)
105 if region is None:
106 xmin = min(t.xmin for t in tiles)
107 xmax = max(t.xmax for t in tiles)
108 ymin = min(t.ymin for t in tiles)
109 ymax = max(t.ymax for t in tiles)
110 else:
111 xmin, xmax, ymin, ymax = positive_region(region)
113 if not multiple_of(xmin, dx):
114 xmin = math.floor(xmin / dx) * dx
115 if not multiple_of(xmax, dx):
116 xmax = math.ceil(xmax / dx) * dx
117 if not multiple_of(ymin, dy):
118 ymin = math.floor(ymin / dy) * dy
119 if not multiple_of(ymax, dy):
120 ymax = math.ceil(ymax / dy) * dy
122 nx = int(round((xmax - xmin) / dx)) + 1
123 ny = int(round((ymax - ymin) / dy)) + 1
125 data = num.zeros((ny, nx), dtype=dtype)
126 data[:, :] = 0
128 for t in tiles:
129 for txmin in (t.xmin, t.xmin + 360.):
130 ix = int(round((txmin - xmin) / dx))
131 iy = int(round((t.ymin - ymin) / dy))
132 ixlo = max(ix, 0)
133 ixhi = min(ix+t.nx, nx)
134 iylo = max(iy, 0)
135 iyhi = min(iy+t.ny, ny)
136 jxlo = ixlo-ix
137 jxhi = jxlo + max(0, ixhi - ixlo)
138 jylo = iylo-iy
139 jyhi = jylo + max(0, iyhi - iylo)
140 if iyhi > iylo and ixhi > ixlo:
141 data[iylo:iyhi, ixlo:ixhi] = t.data[jylo:jyhi, jxlo:jxhi]
143 if not num.any(num.isfinite(data)):
144 return None
146 return Tile(xmin, ymin, dx, dy, data)
149def tile_export(tle, path, format='stl', exaggeration=3., socket_scale=1.,
150 scale=100):
151 '''Export DEM to 3D printable formats'''
152 import trimesh
154 from pyrocko import geometry
155 from pyrocko.cake import earthradius
157 x = tle.x()
158 x = x - x.min()
160 y = tle.y()
161 y = y - y.min()
163 data = tle.data * exaggeration
165 vertices, faces = geometry.topo_to_mesh(
166 y, x, data, earthradius)
168 vertices[:, 0] -= vertices[:, 0].min()
170 nlat = y.size
171 nlon = x.size
172 socket_level = -vertices[:, 0].max() * socket_scale
173 nvertices = vertices.shape[0]
175 south_indices = num.arange(0, nlon)
176 north_indices = num.arange(nvertices - nlon, nvertices)
177 east_indices = num.arange(0, nlat) * nlon
178 west_indices = num.arange(0, nlat) * nlon + nlon - 1
180 north_indices = north_indices[::-1]
181 east_indices = east_indices[::-1]
183 socket_bottom_faces = []
184 for border_indices in [north_indices, east_indices,
185 south_indices, west_indices]:
186 nvertices = vertices.shape[0]
187 nfaces = border_indices.size - 1
189 socket_vertices = vertices[border_indices.astype(int)].copy()
190 socket_vertices[:, 0] = socket_level
192 socket_faces = num.empty((nfaces, 4), dtype=int)
193 for iface in range(nfaces):
194 socket_faces[iface, 0] = iface + nvertices
195 socket_faces[iface, 1] = iface + nvertices + 1
196 socket_faces[iface, 2] = border_indices[iface+1]
197 socket_faces[iface, 3] = border_indices[iface]
199 vertices = num.vstack((vertices, socket_vertices))
200 faces = num.vstack((faces, socket_faces))
202 socket_bottom_faces.append(vertices.shape[0] - nfaces)
204 faces = num.vstack((faces, socket_bottom_faces[::-1]))
206 mesh = trimesh.Trimesh(
207 vertices, faces)
208 mesh.fix_normals()
209 mesh.fill_holes()
211 mesh.apply_scale(scale)
213 with open(path, 'wb') as f:
214 mesh.export(f, file_type=format)
216 return mesh