Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/topo/tile.py: 58%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Topography tile data model. 

8''' 

9 

10import math 

11import numpy as num 

12import scipy.signal 

13 

14from pyrocko.orthodrome import positive_region 

15 

16km = 1e3 

17 

18 

19class OutOfBounds(Exception): 

20 pass 

21 

22 

23class Tile(object): 

24 ''' 

25 2D array of data values which are regularly gridded in lon/lat. 

26 ''' 

27 

28 def __init__(self, xmin, ymin, dx, dy, data): 

29 self.xmin = float(xmin) 

30 self.ymin = float(ymin) 

31 self.dx = float(dx) 

32 self.dy = float(dy) 

33 self.data = data 

34 self._set_maxes() 

35 

36 def _set_maxes(self): 

37 self.ny, self.nx = self.data.shape 

38 self.xmax = self.xmin + (self.nx-1) * self.dx 

39 self.ymax = self.ymin + (self.ny-1) * self.dy 

40 

41 def x(self): 

42 return self.xmin + num.arange(self.nx) * self.dx 

43 

44 def y(self): 

45 return self.ymin + num.arange(self.ny) * self.dy 

46 

47 def decimate(self, ndeci): 

48 assert ndeci % 2 == 0 

49 kernel = num.ones((ndeci+1, ndeci+1)) 

50 kernel /= num.sum(kernel) 

51 data = scipy.signal.convolve2d( 

52 self.data.astype(float), kernel, mode='valid') 

53 

54 self.data = data[::ndeci, ::ndeci].astype(self.data.dtype) 

55 self.xmin += ndeci/2 

56 self.ymin += ndeci/2 

57 self.dx *= ndeci 

58 self.dy *= ndeci 

59 self._set_maxes() 

60 

61 def yextend_with_repeat(self, ymin, ymax): 

62 assert ymax >= self.ymax 

63 assert ymin <= self.ymin 

64 

65 nlo = int(round((self.ymin - ymin) / self.dy)) 

66 nhi = int(round((ymax - self.ymax) / self.dy)) 

67 

68 nx, ny = self.nx, self.ny 

69 data = num.zeros((ny+nlo+nhi, nx), dtype=self.data.dtype) 

70 data[:nlo, :] = self.data[nlo, :] 

71 data[nlo:nlo+ny, :] = self.data 

72 data[nlo+ny:, :] = self.data[-1, :] 

73 

74 self.ymin = ymin 

75 self.data = data 

76 self._set_maxes() 

77 

78 def get(self, x, y): 

79 ix = int(round((x - self.xmin) / self.dx)) 

80 iy = int(round((y - self.ymin) / self.dy)) 

81 if 0 <= ix < self.nx and 0 <= iy < self.ny: 

82 return self.data[iy, ix] 

83 else: 

84 raise OutOfBounds() 

85 

86 @classmethod 

87 def from_grd(cls, path): 

88 '''Load from GMT .grd''' 

89 from pyrocko.plot import gmtpy 

90 lon, lat, z = gmtpy.loadgrd(path) 

91 return cls(lon[0], lat[0], lon[1] - lon[0], lat[1] - lat[0], z) 

92 

93 

94def multiple_of(x, dx, eps=1e-5): 

95 return abs(int(round(x / dx))*dx - x) < dx * eps 

96 

97 

98def combine(tiles, region=None): 

99 if not tiles: 

100 return None 

101 

102 dx = tiles[0].dx 

103 dy = tiles[0].dy 

104 dtype = tiles[0].data.dtype 

105 

106 assert all(t.dx == dx for t in tiles) 

107 assert all(t.dy == dy for t in tiles) 

108 assert all(t.data.dtype == dtype for t in tiles) 

109 assert all(multiple_of(t.xmin, dx) for t in tiles) 

110 assert all(multiple_of(t.ymin, dy) for t in tiles) 

111 

112 if region is None: 

113 xmin = min(t.xmin for t in tiles) 

114 xmax = max(t.xmax for t in tiles) 

115 ymin = min(t.ymin for t in tiles) 

116 ymax = max(t.ymax for t in tiles) 

117 else: 

118 xmin, xmax, ymin, ymax = positive_region(region) 

119 

120 if not multiple_of(xmin, dx): 

121 xmin = math.floor(xmin / dx) * dx 

122 if not multiple_of(xmax, dx): 

123 xmax = math.ceil(xmax / dx) * dx 

124 if not multiple_of(ymin, dy): 

125 ymin = math.floor(ymin / dy) * dy 

126 if not multiple_of(ymax, dy): 

127 ymax = math.ceil(ymax / dy) * dy 

128 

129 nx = int(round((xmax - xmin) / dx)) + 1 

130 ny = int(round((ymax - ymin) / dy)) + 1 

131 

132 data = num.zeros((ny, nx), dtype=dtype) 

133 data[:, :] = 0 

134 

135 for t in tiles: 

136 for txmin in (t.xmin, t.xmin + 360.): 

137 ix = int(round((txmin - xmin) / dx)) 

138 iy = int(round((t.ymin - ymin) / dy)) 

139 ixlo = max(ix, 0) 

140 ixhi = min(ix+t.nx, nx) 

141 iylo = max(iy, 0) 

142 iyhi = min(iy+t.ny, ny) 

143 jxlo = ixlo-ix 

144 jxhi = jxlo + max(0, ixhi - ixlo) 

145 jylo = iylo-iy 

146 jyhi = jylo + max(0, iyhi - iylo) 

147 if iyhi > iylo and ixhi > ixlo: 

148 data[iylo:iyhi, ixlo:ixhi] = t.data[jylo:jyhi, jxlo:jxhi] 

149 

150 if not num.any(num.isfinite(data)): 

151 return None 

152 

153 return Tile(xmin, ymin, dx, dy, data) 

154 

155 

156def tile_export(tle, path, format='stl', exaggeration=3., socket_scale=1., 

157 scale=100): 

158 '''Export DEM to 3D printable formats''' 

159 import trimesh 

160 

161 from pyrocko import geometry 

162 from pyrocko.cake import earthradius 

163 

164 x = tle.x() 

165 x = x - x.min() 

166 

167 y = tle.y() 

168 y = y - y.min() 

169 

170 data = tle.data * exaggeration 

171 

172 vertices, faces = geometry.topo_to_mesh( 

173 y, x, data, earthradius) 

174 

175 vertices[:, 0] -= vertices[:, 0].min() 

176 

177 nlat = y.size 

178 nlon = x.size 

179 socket_level = -vertices[:, 0].max() * socket_scale 

180 nvertices = vertices.shape[0] 

181 

182 south_indices = num.arange(0, nlon) 

183 north_indices = num.arange(nvertices - nlon, nvertices) 

184 east_indices = num.arange(0, nlat) * nlon 

185 west_indices = num.arange(0, nlat) * nlon + nlon - 1 

186 

187 north_indices = north_indices[::-1] 

188 east_indices = east_indices[::-1] 

189 

190 socket_bottom_faces = [] 

191 for border_indices in [north_indices, east_indices, 

192 south_indices, west_indices]: 

193 nvertices = vertices.shape[0] 

194 nfaces = border_indices.size - 1 

195 

196 socket_vertices = vertices[border_indices.astype(int)].copy() 

197 socket_vertices[:, 0] = socket_level 

198 

199 socket_faces = num.empty((nfaces, 4), dtype=int) 

200 for iface in range(nfaces): 

201 socket_faces[iface, 0] = iface + nvertices 

202 socket_faces[iface, 1] = iface + nvertices + 1 

203 socket_faces[iface, 2] = border_indices[iface+1] 

204 socket_faces[iface, 3] = border_indices[iface] 

205 

206 vertices = num.vstack((vertices, socket_vertices)) 

207 faces = num.vstack((faces, socket_faces)) 

208 

209 socket_bottom_faces.append(vertices.shape[0] - nfaces) 

210 

211 faces = num.vstack((faces, socket_bottom_faces[::-1])) 

212 

213 mesh = trimesh.Trimesh( 

214 vertices, faces) 

215 mesh.fix_normals() 

216 mesh.fill_holes() 

217 

218 mesh.apply_scale(scale) 

219 

220 with open(path, 'wb') as f: 

221 mesh.export(f, file_type=format) 

222 

223 return mesh