1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import numpy as num
9from pyrocko.guts import Bool, Float, StringChoice
10from pyrocko import cake, automap, plot
11from pyrocko.dataset import topo
12from pyrocko.gui.qt_compat import qw, qc
14from pyrocko.gui.vtk_util import TrimeshPipe, faces_to_cells, \
15 cpt_to_vtk_lookuptable
17from .base import Element, ElementState
18from .. import state as vstate
19from pyrocko import geometry
21from .. import common
23guts_prefix = 'sparrow'
26def ticks(vmin, vmax, vstep):
27 vmin = num.floor(vmin / vstep) * vstep
28 vmax = num.ceil(vmax / vstep) * vstep
29 n = int(round((vmax - vmin) / vstep))
30 return vmin + num.arange(n+1) * vstep
33def nice_value_circle(step):
34 step = plot.nice_value(step)
35 if step > 30.:
36 return 30.
38 return step
41class TopoMeshPipe(TrimeshPipe):
43 def __init__(
44 self, tile, cells_cache=None,
45 mask_ocean=False, mask_land=False, **kwargs):
47 vertices, faces = geometry.topo_to_mesh(
48 tile.y(), tile.x(), tile.data,
49 cake.earthradius)
51 self._exaggeration = 1.0
52 self._tile = tile
53 self._raw_vertices = vertices
55 centers = geometry.face_centers(vertices, faces)
57 altitudes = (geometry.vnorm(centers) - 1.0) * cake.earthradius
58 if mask_ocean:
59 mask = num.all(tile.data.flatten()[faces] == 0, axis=1)
60 altitudes[mask] = None
62 if mask_land:
63 mask = num.all(tile.data.flatten()[faces] >= 0, axis=1)
64 altitudes[mask] = None
66 if cells_cache is not None:
67 if id(faces) not in cells_cache:
68 cells_cache[id(faces)] = faces_to_cells(faces)
70 cells = cells_cache[id(faces)]
71 else:
72 cells = faces_to_cells(faces)
74 TrimeshPipe.__init__(
75 self, self._raw_vertices, cells=cells, values=altitudes, **kwargs)
77 def set_exaggeration(self, exaggeration):
78 if self._exaggeration != exaggeration:
79 factors = \
80 (cake.earthradius + self._tile.data.flatten()*exaggeration) \
81 / (cake.earthradius + self._tile.data.flatten())
82 self.set_vertices(self._raw_vertices * factors[:, num.newaxis])
83 self._exaggeration = exaggeration
86class TopoCPTChoice(StringChoice):
87 choices = ['light', 'uniform']
90class TopoState(ElementState):
91 visible = Bool.T(default=True)
92 exaggeration = Float.T(default=1.0)
93 opacity = Float.T(default=1.0)
94 smooth = Bool.T(default=False)
95 cpt = TopoCPTChoice.T(default='light')
96 shading = vstate.ShadingChoice.T(default='phong')
97 resolution_max_factor = Float.T(default=1.0)
98 resolution_min_factor = Float.T(default=1.0)
99 coverage_factor = Float.T(default=1.0)
101 def create(self):
102 element = TopoElement()
103 return element
106class TopoElement(Element):
108 def __init__(self):
109 Element.__init__(self)
110 self._parent = None
111 self.mesh = None
112 self._controls = None
114 # region = (-180., 180, -90, 90)
115 # self._tile = topo.get('ETOPO1_D8', region)
117 self._visible = False
118 self._active_meshes = {}
119 self._meshes = {}
120 self._cells = {}
121 self._cpt_name = None
122 self._lookuptables = {}
124 def get_name(self):
125 return 'Topography'
127 def bind_state(self, state):
128 Element.bind_state(self, state)
129 for var in ['visible', 'exaggeration', 'opacity', 'smooth', 'shading',
130 'cpt', 'resolution_min_factor', 'resolution_max_factor',
131 'coverage_factor']:
132 self.register_state_listener3(self.update, state, var)
134 self._state = state
136 def unbind_state(self):
137 self._listeners.clear()
139 def set_parent(self, parent):
140 self._parent = parent
142 for var in ['distance', 'lat', 'lon']:
143 self.register_state_listener3(
144 self.update, self._parent.state, var)
146 self._parent.add_panel(
147 self.get_name(),
148 self._get_controls(),
149 visible=True,
150 title_controls=[
151 self.get_title_control_remove(),
152 self.get_title_control_visible()])
154 self.update()
156 def unset_parent(self):
157 self.unbind_state()
158 if self._parent:
159 for k in self._active_meshes:
160 self._parent.remove_actor(self._active_meshes[k].actor)
162 self._active_meshes.clear()
163 self._meshes.clear()
164 self._cells.clear()
166 if self._controls:
167 self._parent.remove_panel(self._controls)
168 self._controls = None
170 self._parent.update_view()
171 self._parent = None
173 def select_dems(self, delta, region):
174 if not self._parent:
175 return [], []
177 _, size_y = self._parent.renwin.GetSize()
179 dmin = 2.0 * delta * 1.0 / float(size_y) \
180 / self._state.resolution_max_factor # [deg]
181 dmax = 2.0 * delta * 20.0 \
182 * self._state.resolution_min_factor / float(size_y)
184 result = [
185 topo.select_dem_names(
186 k, dmin, dmax, topo.positive_region(region), mode='highest')
188 for k in ['ocean', 'land']]
190 if not any(result) and delta > 20.:
191 return ['ETOPO1_D8'], ['ETOPO1_D8']
192 else:
193 return result
195 def update_cpt(self, cpt_name):
196 if cpt_name not in self._lookuptables:
197 if cpt_name == 'light':
198 topo_cpt_wet = 'light_sea'
199 topo_cpt_dry = 'light_land'
201 elif cpt_name == 'uniform':
202 topo_cpt_wet = 'light_sea_uniform'
203 topo_cpt_dry = 'light_land_uniform'
205 cpt_wet = automap.read_cpt(topo.cpt(topo_cpt_wet))
206 cpt_dry = automap.read_cpt(topo.cpt(topo_cpt_dry))
207 cpt_combi = automap.cpt_merge_wet_dry(cpt_wet, cpt_dry)
209 lut_combi = cpt_to_vtk_lookuptable(cpt_combi)
210 lut_combi.SetNanColor(0.0, 0.0, 0.0, 0.0)
212 self._lookuptables[cpt_name] = lut_combi
214 def update(self, *args):
216 pstate = self._parent.state
217 delta = pstate.distance * geometry.r2d * 0.5
218 visible = self._state.visible
220 self.update_cpt(self._state.cpt)
222 step = nice_value_circle(
223 max(1./8., min(2**round(math.log(delta) / math.log(2.)), 10.)))
225 lat_min, lat_max, lon_min, lon_max, lon_closed = common.cover_region(
226 pstate.lat, pstate.lon, delta*self._state.coverage_factor, step)
228 if lon_closed:
229 lon_max = 180.
231 region = lon_min, lon_max, lat_min, lat_max
233 dems_ocean, dems_land = self.select_dems(delta, region)
235 lat_majors = ticks(lat_min, lat_max-step, step)
236 lon_majors = ticks(lon_min, lon_max-step, step)
238 wanted = set()
239 if visible:
240 for ilat, lat in enumerate(lat_majors):
241 for ilon, lon in enumerate(lon_majors):
242 lon = ((lon + 180.) % 360.) - 180.
244 region = topo.positive_region(
245 (lon, lon+step, lat, lat+step))
247 for demname in dems_land[:1] + dems_ocean[:1]:
248 mask_ocean = demname.startswith('SRTM') \
249 or demname.startswith('Iceland')
251 mask_land = demname.startswith('ETOPO') \
252 and (dems_land and dems_land[0] != demname)
254 k = (step, demname, region, mask_ocean, mask_land)
255 if k not in self._meshes:
256 tile = topo.get(demname, region)
257 if not tile:
258 continue
260 self._meshes[k] = TopoMeshPipe(
261 tile,
262 cells_cache=self._cells,
263 mask_ocean=mask_ocean,
264 mask_land=mask_land,
265 smooth=self._state.smooth,
266 lut=self._lookuptables[self._state.cpt])
268 wanted.add(k)
270 # prevent adding both, SRTM and ETOPO because
271 # vtk produces artifacts when showing the two masked
272 # meshes (like this, mask_land is never used):
273 break
275 unwanted = set()
276 for k in self._active_meshes:
277 if k not in wanted or not visible:
278 unwanted.add(k)
280 for k in unwanted:
281 self._parent.remove_actor(self._active_meshes[k].actor)
282 del self._active_meshes[k]
284 for k in wanted:
285 if k not in self._active_meshes:
286 m = self._meshes[k]
287 self._active_meshes[k] = m
288 self._parent.add_actor(m.actor)
290 self._active_meshes[k].set_exaggeration(self._state.exaggeration)
291 self._active_meshes[k].set_opacity(self._state.opacity)
292 self._active_meshes[k].set_smooth(self._state.smooth)
293 self._active_meshes[k].set_shading(self._state.shading)
294 self._active_meshes[k].set_lookuptable(
295 self._lookuptables[self._state.cpt])
297 self._parent.update_view()
299 def _get_controls(self):
300 state = self._state
301 if not self._controls:
302 from ..state import state_bind_slider, state_bind_checkbox, \
303 state_bind_combobox
305 frame = qw.QFrame()
306 layout = qw.QGridLayout()
307 frame.setLayout(layout)
309 iy = 0
311 # exaggeration
313 layout.addWidget(qw.QLabel('Exaggeration'), iy, 0)
315 slider = qw.QSlider(qc.Qt.Horizontal)
316 slider.setSizePolicy(
317 qw.QSizePolicy(
318 qw.QSizePolicy.Expanding, qw.QSizePolicy.Fixed))
319 slider.setMinimum(0)
320 slider.setMaximum(2000)
321 layout.addWidget(slider, iy, 1)
323 state_bind_slider(self, state, 'exaggeration', slider, factor=0.01)
325 iy += 1
327 # opacity
329 layout.addWidget(qw.QLabel('Opacity'), iy, 0)
331 slider = qw.QSlider(qc.Qt.Horizontal)
332 slider.setSizePolicy(
333 qw.QSizePolicy(
334 qw.QSizePolicy.Expanding, qw.QSizePolicy.Fixed))
335 slider.setMinimum(0)
336 slider.setMaximum(1000)
337 layout.addWidget(slider, iy, 1)
339 state_bind_slider(self, state, 'opacity', slider, factor=0.001)
341 iy += 1
343 # high resolution
345 layout.addWidget(qw.QLabel('High-Res Factor'), iy, 0)
347 slider = qw.QSlider(qc.Qt.Horizontal)
348 slider.setSizePolicy(
349 qw.QSizePolicy(
350 qw.QSizePolicy.Expanding, qw.QSizePolicy.Fixed))
351 slider.setMinimum(500)
352 slider.setMaximum(4000)
353 layout.addWidget(slider, iy, 1)
355 state_bind_slider(
356 self, state, 'resolution_max_factor', slider, factor=0.001)
358 iy += 1
360 # low resolution
362 layout.addWidget(qw.QLabel('Low-Res Factor'), iy, 0)
364 slider = qw.QSlider(qc.Qt.Horizontal)
365 slider.setSizePolicy(
366 qw.QSizePolicy(
367 qw.QSizePolicy.Expanding, qw.QSizePolicy.Fixed))
368 slider.setMinimum(500)
369 slider.setMaximum(4000)
370 layout.addWidget(slider, iy, 1)
372 state_bind_slider(
373 self, state, 'resolution_min_factor', slider, factor=0.001)
375 iy += 1
377 # low resolution
379 layout.addWidget(qw.QLabel('Coverage Factor'), iy, 0)
381 slider = qw.QSlider(qc.Qt.Horizontal)
382 slider.setSizePolicy(
383 qw.QSizePolicy(
384 qw.QSizePolicy.Expanding, qw.QSizePolicy.Fixed))
385 slider.setMinimum(500)
386 slider.setMaximum(4000)
387 layout.addWidget(slider, iy, 1)
389 state_bind_slider(
390 self, state, 'coverage_factor', slider, factor=0.001)
392 iy += 1
394 cb = common.string_choices_to_combobox(TopoCPTChoice)
395 layout.addWidget(qw.QLabel('CPT'), iy, 0)
396 layout.addWidget(cb, iy, 1)
397 state_bind_combobox(self, state, 'cpt', cb)
399 iy += 1
401 cb = qw.QCheckBox('Smooth')
402 layout.addWidget(cb, iy, 0)
403 state_bind_checkbox(self, state, 'smooth', cb)
405 cb = common.string_choices_to_combobox(vstate.ShadingChoice)
406 layout.addWidget(cb, iy, 1)
407 state_bind_combobox(self, state, 'shading', cb)
409 iy += 1
411 layout.addWidget(qw.QFrame(), iy, 0, 1, 2)
413 self._controls = frame
415 return self._controls
418__all__ = [
419 'TopoElement',
420 'TopoState',
421]