# https://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
self.data.astype(num.float), kernel, mode='valid')
assert ymax >= self.ymax assert ymin <= self.ymin
nlo = int(round((self.ymin - ymin) / self.dy)) nhi = int(round((ymax - self.ymax) / self.dy))
nx, ny = self.nx, self.ny data = num.zeros((ny+nlo+nhi, nx), dtype=self.data.dtype) data[:nlo, :] = self.data[nlo, :] data[nlo:nlo+ny, :] = self.data data[nlo+ny:, :] = self.data[-1, :]
self.ymin = ymin self.data = data self._set_maxes()
ix = int(round((x - self.xmin) / self.dx)) iy = int(round((y - self.ymin) / self.dy)) if 0 <= ix < self.nx and 0 <= iy < self.ny: return self.data[iy, ix] else: raise OutOfBounds()
def from_grd(cls, path): '''Load from GMT .grd''' from pyrocko.plot import gmtpy lon, lat, z = gmtpy.loadgrd(path) return cls(lon[0], lat[0], lon[1] - lon[0], lat[1] - lat[0], z)
return None
else:
return None
scale=100): '''Export DEM to 3D printable formats''' import trimesh
from pyrocko import geometry from pyrocko.cake import earthradius
x = tle.x() x = x - x.min()
y = tle.y() y = y - y.min()
data = tle.data * exaggeration
vertices, faces = geometry.topo_to_mesh( y, x, data, earthradius)
vertices[:, 0] -= vertices[:, 0].min()
nlat = y.size nlon = x.size socket_level = -vertices[:, 0].max() * socket_scale nvertices = vertices.shape[0]
south_indices = num.arange(0, nlon) north_indices = num.arange(nvertices - nlon, nvertices) east_indices = num.arange(0, nlat) * nlon west_indices = num.arange(0, nlat) * nlon + nlon - 1
north_indices = north_indices[::-1] east_indices = east_indices[::-1]
socket_bottom_faces = [] for border_indices in [north_indices, east_indices, south_indices, west_indices]: nvertices = vertices.shape[0] nfaces = border_indices.size - 1
socket_vertices = vertices[border_indices.astype(num.int)].copy() socket_vertices[:, 0] = socket_level
socket_faces = num.empty((nfaces, 4), dtype=num.int) for iface in range(nfaces): socket_faces[iface, 0] = iface + nvertices socket_faces[iface, 1] = iface + nvertices + 1 socket_faces[iface, 2] = border_indices[iface+1] socket_faces[iface, 3] = border_indices[iface]
vertices = num.vstack((vertices, socket_vertices)) faces = num.vstack((faces, socket_faces))
socket_bottom_faces.append(vertices.shape[0] - nfaces)
faces = num.vstack((faces, socket_bottom_faces[::-1]))
mesh = trimesh.Trimesh( vertices, faces) mesh.fix_normals() mesh.fill_holes()
mesh.apply_scale(scale)
with open(path, 'wb') as f: mesh.export(f, file_type=format)
return mesh |