diff --git a/tools/python/tests/test_udbase_vis.py b/tools/python/tests/test_udbase_vis.py index d15fd825a..9f8022795 100644 --- a/tools/python/tests/test_udbase_vis.py +++ b/tools/python/tests/test_udbase_vis.py @@ -91,6 +91,14 @@ def plot_veg(self, **kwargs): self.assertEqual(result, "plot_veg_result") self.assertEqual(sim.vis.calls, [("plot_veg", {"show": False})]) + # regression: the aliases must forward backend/color/opacity kwargs + # instead of raising TypeError on cases exercising them (issue #337) + sim.vis.calls.clear() + sim.plot_trees(show=False, backend="plotly", opacity=0.5) + self.assertEqual(sim.vis.calls, [ + ("plot_veg", {"show": False, "backend": "plotly", "opacity": 0.5}), + ]) + def test_udbase_plot_2dmap_forwards_to_vis_facade(self): sim = UDBase.__new__(UDBase) sim.vis = RecordingVis() @@ -557,6 +565,7 @@ def _vis_with_stub_sim(tmp_path): sim = SimpleNamespace( _lfgeom=True, geom=geom, expnr="001", path=Path(tmp_path), xt=grid, yt=grid, zt=grid, xm=grid, ym=grid, zm=grid, + dx=0.4, dy=0.4, dzt=_np.full(6, 0.4), ffacets="facets", ffactypes="factypes", veg={"points": _np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]])}, Sc=_np.zeros((6, 6, 6), dtype=bool), @@ -588,6 +597,7 @@ def test_pyvista_backend_builds_all_overlay_plots(self): # These default to the instance backend ("pyvista"); no per-call backend. for label, plotter in [ ("plot_veg", vis.plot_veg(show=False)), + ("plot_veg_outline", vis.plot_veg_outline(show=False)), ("plot_scalar_source", vis.plot_scalar_source(show=False)), ("plot_solid", vis.plot_solid("c", show=False)), ("plot_fluid_boundary", vis.plot_fluid_boundary("c", show=False)), @@ -622,5 +632,166 @@ def test_pyvista_backend_missing_names_install_command(self): self.assertIn("pyvista", msg.lower()) self.assertIn("pip install", msg) + +class TestNiceTicks(unittest.TestCase): + """Tick locator used by the hand-drawn PyVista axes (issue #337).""" + + def test_round_steps_for_typical_domains(self): + from udvis.scene import _nice_ticks + + values, decimals = _nice_ticks(0.0, 340.0) + self.assertEqual(values, [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0]) + self.assertEqual(decimals, 0) + + values, decimals = _nice_ticks(0.0, 240.0) + self.assertEqual(values, [0.0, 50.0, 100.0, 150.0, 200.0]) + self.assertEqual(decimals, 0) + + def test_z_range_gets_intermediate_ticks(self): + from udvis.scene import _nice_ticks + + values, decimals = _nice_ticks(0.0, 104.0) + self.assertEqual(values, [0.0, 20.0, 40.0, 60.0, 80.0, 100.0]) + self.assertEqual(decimals, 0) + + def test_at_most_max_ticks_values(self): + from udvis.scene import _nice_ticks + + for vmax in (1.0, 7.3, 55.0, 104.0, 240.0, 333.0, 1024.0): + values, _ = _nice_ticks(0.0, vmax, max_ticks=8) + self.assertLessEqual(len(values), 8) + self.assertGreaterEqual(len(values), 4) + self.assertTrue(all(0.0 <= v <= vmax for v in values)) + + def test_fractional_steps_report_decimals(self): + from udvis.scene import _nice_ticks + + values, decimals = _nice_ticks(0.0, 1.0) + self.assertEqual(values, [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]) + self.assertEqual(decimals, 1) + + def test_offset_range_aligns_to_step_multiples(self): + from udvis.scene import _nice_ticks + + values, _ = _nice_ticks(-104.0, 104.0) + self.assertIn(0.0, values) + self.assertEqual(values[0], -100.0) + self.assertEqual(values[-1], 100.0) + + def test_degenerate_range_returns_single_tick(self): + from udvis.scene import _nice_ticks + + values, decimals = _nice_ticks(5.0, 5.0) + self.assertEqual(values, [5.0]) + self.assertEqual(decimals, 0) + + +class TestPlotlySceneAxes(unittest.TestCase): + """Plotly 3-D axes must draw visible lines, clear of the data (issue #337).""" + + @unittest.skipUnless( + importlib.util.find_spec("plotly"), "plotly backend is optional and not installed" + ) + def test_axis_lines_shown_on_bounds_with_label_standoff(self): + from udvis.scene import MeshPrimitive, Scene, render_scene + + verts = np.array([[0.0, 0.0, 0.0], [10.0, 0.0, 0.0], [10.0, 20.0, 5.0], + [0.0, 20.0, 5.0]]) + faces = np.array([[0, 1, 2], [0, 2, 3]]) + scene = Scene(meshes=[MeshPrimitive(verts, faces)], + bounds=(verts.min(axis=0), verts.max(axis=0))) + fig = render_scene(scene, backend="plotly", show=False) + + for name, hi in (("xaxis", 10.0), ("yaxis", 20.0), ("zaxis", 5.0)): + axis = getattr(fig.layout.scene, name) + self.assertTrue(axis.showline, name) + # axes anchored exactly on the domain bounds (no padding); the + # long outside ticks provide the label standoff instead + self.assertEqual(axis.range[0], 0.0, name) + self.assertEqual(axis.range[1], hi, name) + self.assertEqual(axis.ticks, "outside", name) + self.assertGreaterEqual(axis.ticklen, 8, name) + + +class TestVegVoxelMesh(unittest.TestCase): + """plot_veg draws each vegetation point as a box filling its grid cell.""" + + def test_single_cell_spans_exact_cell_bounds(self): + from udvis.udbase_vis import _veg_voxel_mesh + + xm = np.array([0.0, 2.0, 4.0]) + ym = np.array([0.0, 3.0]) + zm = np.array([0.0, 1.0, 2.5]) + dzt = np.array([1.0, 1.5, 2.0]) + verts, faces = _veg_voxel_mesh([[1, 0, 2]], xm, ym, zm, dx=2.0, dy=3.0, dzt=dzt) + self.assertEqual(verts.shape, (8, 3)) + self.assertEqual(faces.shape, (12, 3)) + np.testing.assert_allclose(verts.min(axis=0), [2.0, 0.0, 2.5]) + np.testing.assert_allclose(verts.max(axis=0), [4.0, 3.0, 4.5]) # z1 = zm[2] + dzt[2] + + def test_multiple_cells_index_into_own_vertex_block(self): + from udvis.udbase_vis import _veg_voxel_mesh + + grid = np.arange(5, dtype=float) + verts, faces = _veg_voxel_mesh( + [[0, 0, 0], [3, 2, 1]], grid, grid, grid, dx=1.0, dy=1.0, + dzt=np.ones(5)) + self.assertEqual(verts.shape, (16, 3)) + self.assertEqual(faces.shape, (24, 3)) + # second box's triangles reference only its own 8 vertices + self.assertTrue((faces[12:] >= 8).all() and (faces[12:] < 16).all()) + + def test_boxes_are_closed_watertight_triangulations(self): + import trimesh + + from udvis.udbase_vis import _veg_voxel_mesh + + grid = np.arange(3, dtype=float) + verts, faces = _veg_voxel_mesh([[1, 1, 1]], grid, grid, grid, + dx=1.0, dy=1.0, dzt=np.ones(3)) + box = trimesh.Trimesh(vertices=verts, faces=faces, process=False) + self.assertTrue(box.is_watertight) + self.assertGreater(box.volume, 0) # outward winding + + def test_empty_points_give_empty_mesh(self): + from udvis.udbase_vis import _veg_voxel_mesh + + verts, faces = _veg_voxel_mesh( + np.empty((0, 3)), np.array([0.0]), np.array([0.0]), np.array([0.0]), + dx=1.0, dy=1.0, dzt=np.array([1.0])) + self.assertEqual(len(verts), 0) + self.assertEqual(len(faces), 0) + + +class TestZTickValues(unittest.TestCase): + """Vertical-axis ticks skip the base label that collides with the y axis.""" + + def test_base_tick_dropped_at_exact_zero(self): + from udvis.scene import _z_tick_values + + values, _ = _z_tick_values(0.0, 104.0) + self.assertEqual(values, [20.0, 40.0, 60.0, 80.0, 100.0]) + + def test_base_tick_dropped_despite_float_noise_below_zero(self): + # STL meshes commonly have zmin at -1e-9 rather than exactly 0; the + # base tick must still be excluded (regression for issue #337). + from udvis.scene import _z_tick_values + + values, _ = _z_tick_values(-1.28e-9, 104.42) + self.assertEqual(values, [20.0, 40.0, 60.0, 80.0, 100.0]) + + def test_first_tick_kept_for_elevated_ranges(self): + from udvis.scene import _z_tick_values + + values, _ = _z_tick_values(50.0, 104.0) + self.assertEqual(values[0], 60.0) + + def test_degenerate_range_yields_no_ticks(self): + from udvis.scene import _z_tick_values + + values, _ = _z_tick_values(0.0, 0.0) + self.assertEqual(values, []) + + if __name__ == "__main__": unittest.main() diff --git a/tools/python/udbase.py b/tools/python/udbase.py index 34323deed..6bdb1e411 100644 --- a/tools/python/udbase.py +++ b/tools/python/udbase.py @@ -1306,14 +1306,15 @@ def coarsegrain_field(var: np.ndarray, Lflt: np.ndarray, return udstats.coarsegrain_field(var, Lflt, xm, ym) def plot_veg(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, - backend: Optional[str] = None): - """Plot vegetation points on top of the geometry using the visualization facade. + backend: Optional[str] = None, **kwargs): + """Plot vegetation on top of the geometry using the visualization facade. ``show`` defaults to ``False`` (returns the figure/plotter without displaying it), consistent with the other overlay plots - (``plot_scalar_source``, ``plot_trees``). + (``plot_scalar_source``, ``plot_trees``). Remaining keyword arguments + (e.g. ``color``, ``opacity``) are forwarded to :meth:`UDVis.plot_veg`. """ - return self.vis.plot_veg(veg=veg, show=show, backend=backend) + return self.vis.plot_veg(veg=veg, show=show, backend=backend, **kwargs) def plot_scalar_source( self, @@ -1334,14 +1335,28 @@ def plot_scalar_source( backend=backend, ) - def plot_trees(self, show: bool = False): + def plot_veg_outline(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, + backend: Optional[str] = None, **kwargs): + """Plot vegetation voxels over outline-style geometry via the facade. + + Same vegetation rendering as :meth:`plot_veg` but on the outline base + of ``show_outline`` instead of the full wireframe. ``show`` defaults + to ``False``; remaining keyword arguments (``color``, ``opacity``, + ``line_width``, ``angle_threshold``) are forwarded to + :meth:`UDVis.plot_veg_outline`. + """ + return self.vis.plot_veg_outline(veg=veg, show=show, backend=backend, **kwargs) + + def plot_trees(self, show: bool = False, **kwargs): """Backward-compatible alias for :meth:`plot_veg`. ``UDVis`` exposes only ``plot_veg`` (vegetation is the current name for what legacy cases called "trees"), so forward there rather than to a - nonexistent ``UDVis.plot_trees``. ``show`` defaults to ``False``. + nonexistent ``UDVis.plot_trees``. ``show`` defaults to ``False``; + remaining keyword arguments (``backend``, ``color``, ``opacity``, ...) + are forwarded unchanged. """ - return self.vis.plot_veg(show=show) + return self.vis.plot_veg(show=show, **kwargs) def plot_fac(self, var: np.ndarray, building_ids: Optional[np.ndarray] = None, show: bool = True, backend: Optional[str] = None): diff --git a/tools/python/udvis/scene.py b/tools/python/udvis/scene.py index f2078c21d..d2fa7bdf7 100644 --- a/tools/python/udvis/scene.py +++ b/tools/python/udvis/scene.py @@ -93,10 +93,16 @@ class LineSet: segments: np.ndarray # (S, 2) int color: str = "black" width: float = 2.0 + opacity: float = 1.0 name: Optional[str] = None # Plotly cosmetic: lift edges that sit on the ground plane (z≈0) slightly so # they are not hidden by the ground mesh. lift_ground: bool = False + # Optional Plotly-specific overrides: WebGL lines rasterize much cruder + # than VTK's, so widths/opacities tuned for PyVista often need different + # values in Plotly. None means "use width/opacity". + width_plotly: Optional[float] = None + opacity_plotly: Optional[float] = None @dataclass @@ -179,6 +185,12 @@ def compute_bounds(self) -> Tuple[np.ndarray, np.ndarray]: # optional lightweight backend (``backend="plotly"``). DEFAULT_BACKEND = "pyvista" +# Supersampling ratio for notebook (trame) still frames. Text is rendered at +# a fixed pixel size by VTK, so every font size in the PyVista renderer is +# scaled by this same factor to keep its apparent size once the browser +# scales the frame back down to the widget size. +PYVISTA_STILL_RATIO = 1.5 + def normalize_backend(backend: Optional[str]) -> str: """Return the lower-cased backend name, validating it is supported. @@ -224,6 +236,29 @@ def render_scene(scene: Scene, backend: str = DEFAULT_BACKEND, show: bool = True _PLOTLY_COLORSCALE = {"viridis": "Viridis", "greys": "Greys", "greys_r": "Greys_r"} +def _bake_face_shading(vertices, faces, rgba, light_dir=(-0.45, -0.45, 0.78), + ambient=0.55, diffuse=0.45): + """Multiply a directional Lambert term into per-face colours. + + plotly.js does not apply its lighting model to Mesh3d ``facecolor`` + arrays, so meshes that should look shaded need the shading baked into + the colours themselves. Two-sided (|n.l|) like VTK's default, so STL + winding direction cannot black out faces. + """ + v = np.asarray(vertices, dtype=float) + f = np.asarray(faces, dtype=int) + n = np.cross(v[f[:, 1]] - v[f[:, 0]], v[f[:, 2]] - v[f[:, 0]]) + norms = np.linalg.norm(n, axis=1) + norms[norms == 0] = 1.0 + n /= norms[:, None] + ld = np.asarray(light_dir, dtype=float) + ld /= np.linalg.norm(ld) + shade = ambient + diffuse * np.abs(n @ ld) + out = np.array(rgba, dtype=float, copy=True) + out[:, :3] *= shade[:, None] + return np.clip(out, 0.0, 1.0) + + def _render_plotly(scene: Scene, show: bool = True): try: import plotly.graph_objects as go @@ -244,21 +279,25 @@ def _render_plotly(scene: Scene, show: bool = True): continue verts = np.asarray(mesh.vertices, dtype=float) rgba = mesh.resolved_face_colors() + if mesh.lighting != "flat": + # "flat" meshes carry pre-shaded colours (e.g. the two-tone + # outline views); everything else gets shading baked in, since + # plotly.js ignores its lighting model for facecolor meshes. + rgba = _bake_face_shading(verts, faces, rgba) opacity = float(mesh.opacity) trace = go.Mesh3d( x=verts[:, 0], y=verts[:, 1], z=verts[:, 2], i=faces[:, 0], j=faces[:, 1], k=faces[:, 2], facecolor=[f"rgb({int(round(255 * c[0]))},{int(round(255 * c[1]))},{int(round(255 * c[2]))})" for c in rgba], opacity=opacity, - flatshading=True, + # Ambient-only, set explicitly: some plotly versions apply + # partial lighting to facecolor meshes, most apply none. All + # shading intent is already baked into the colours here. + flatshading=False, + lighting=dict(ambient=1.0, diffuse=0.0, specular=0.0, roughness=1.0, fresnel=0.0), + lightposition=dict(x=0, y=0, z=1), name=mesh.name, ) - if mesh.lighting == "flat": - trace.update( - flatshading=False, - lighting=dict(ambient=1.0, diffuse=0.0, specular=0.0, roughness=1.0, fresnel=0.0), - lightposition=dict(x=0, y=0, z=1), - ) fig.add_trace(trace) for ln in scene.lines: @@ -274,9 +313,16 @@ def _render_plotly(scene: Scene, show: bool = True): ex.extend([p0[0], p1[0], None]) ey.extend([p0[1], p1[1], None]) ez.extend([z0, z1, None]) + # Plotly's WebGL lines rasterize very differently from VTK's (no + # depth offset, cruder joins), so LineSets carry optional + # plotly-specific width/opacity to keep the original bold look here + # while PyVista uses thinner, translucent wireframes. + p_width = ln.width if ln.width_plotly is None else ln.width_plotly + p_opacity = ln.opacity if ln.opacity_plotly is None else ln.opacity_plotly fig.add_trace(go.Scatter3d( x=ex, y=ey, z=ez, mode="lines", - line=dict(color=ln.color, width=ln.width), + line=dict(color=ln.color, width=p_width), + opacity=float(p_opacity), name=ln.name, showlegend=False, hoverinfo="skip", )) @@ -353,6 +399,17 @@ def _render_plotly(scene: Scene, show: bool = True): mins, maxs = scene.compute_bounds() az, el, dist = np.deg2rad(225.0), np.deg2rad(20.0), 1.75 lx, ly, lz = scene.axis_labels + x_range = [float(mins[0]), float(maxs[0])] + y_range = [float(mins[1]), float(maxs[1])] + z_range = [float(mins[2]), float(maxs[2])] + # Explicit axis lines and outside ticks: plotly 3-D axes draw no line at + # all by default, leaving the z axis in particular invisible. The axes sit + # exactly on the domain bounds; the long outside ticks are what keep the + # tick labels clear of the plotted field (3-D axes have no label + # standoff property). + _axis_common = dict(showgrid=False, showbackground=False, showline=True, + linecolor="black", linewidth=2, ticks="outside", + tickcolor="black", ticklen=10) # visible=False hides the axis line, ticks, labels and title together. _axis_extra = {} if scene.show_axes else dict(visible=False) _legend = dict(title=scene.legend_title, itemsizing="constant") if scene.legend else {} @@ -363,9 +420,9 @@ def _render_plotly(scene: Scene, show: bool = True): margin=dict(l=0, r=0, b=0, t=40, pad=0), scene=dict( aspectmode="data", - xaxis=dict(title=lx, range=[float(mins[0]), float(maxs[0])], showgrid=False, showbackground=False, **_axis_extra), - yaxis=dict(title=ly, range=[float(mins[1]), float(maxs[1])], showgrid=False, showbackground=False, **_axis_extra), - zaxis=dict(title=lz, range=[float(mins[2]), float(maxs[2])], showgrid=False, showbackground=False, **_axis_extra), + xaxis=dict(title=lx, range=x_range, **_axis_common, **_axis_extra), + yaxis=dict(title=ly, range=y_range, **_axis_common, **_axis_extra), + zaxis=dict(title=lz, range=z_range, **_axis_common, **_axis_extra), camera=dict( projection=dict(type="orthographic"), eye=dict( @@ -435,6 +492,41 @@ def _render_pyvista(scene: Scene, show: bool = True): plotter = pv.Plotter(image_scale=4) plotter.set_background("white") plotter.enable_parallel_projection() + try: + # 8-sample MSAA: sub-pixel edge sampling keeps outlines crisp where + # the shader-based FXAA visibly softens them. (A widget hang once + # attributed to MSAA turned out to be VS Code webview corruption.) + # FXAA remains the fallback for GL stacks without multisample + # support. + plotter.enable_anti_aliasing("msaa", multi_samples=16) + except Exception: + try: + plotter.enable_anti_aliasing("fxaa") + except Exception: + pass + try: + # Screen-space ambient occlusion: contact shading at building bases + # and in street canyons that gives the flat-shaded scene depth. The + # radius scales with the scene so the effect survives any domain + # size. + _mins, _maxs = scene.compute_bounds() + _span = float(np.max(np.asarray(_maxs, dtype=float) - np.asarray(_mins, dtype=float))) + if _span > 0: + plotter.enable_ssao(radius=0.02 * _span, bias=0.005 * _span, + kernel_size=128, blur=True) + except Exception: + pass + try: + # Stream notebook stills above the widget's CSS size and at higher + # JPEG quality: at the pyvista defaults (ratio 1, quality 85) any + # OS/browser scaling above 100% shows an undersampled, soft image at + # rest. plotter._theme is a per-plotter copy, so nothing leaks + # globally. Font sizes below are scaled by the same ratio. + plotter._theme.trame.still_ratio = PYVISTA_STILL_RATIO + plotter._theme.trame.jpeg_quality = 95 + except AttributeError: + pass + ts = PYVISTA_STILL_RATIO for mesh in scene.meshes: faces = np.asarray(mesh.faces, dtype=np.int64) @@ -442,7 +534,8 @@ def _render_pyvista(scene: Scene, show: bool = True): continue verts = np.asarray(mesh.vertices, dtype=float) poly = pv.PolyData(verts, _faces_to_pyvista(faces)) - common = dict(show_edges=mesh.show_edges, name=mesh.name or None) + common = dict(show_edges=mesh.show_edges, name=mesh.name or None, + opacity=float(mesh.opacity)) if mesh.show_edges: common.update(edge_color=mesh.edge_color, line_width=mesh.edge_width) if mesh.scalars is not None: @@ -478,7 +571,8 @@ def _render_pyvista(scene: Scene, show: bool = True): line_poly = pv.PolyData() line_poly.points = verts line_poly.lines = cells.ravel() - plotter.add_mesh(line_poly, color=_pyvista_color(ln.color), line_width=ln.width, name=ln.name or None) + plotter.add_mesh(line_poly, color=_pyvista_color(ln.color), line_width=ln.width, + opacity=float(ln.opacity), name=ln.name or None) for ps in scene.points: pts = np.asarray(ps.points, dtype=float) @@ -489,7 +583,7 @@ def _render_pyvista(scene: Scene, show: bool = True): render_points_as_spheres=True, name=ps.name or None) if ps.labels: plotter.add_point_labels( - pts, list(ps.labels), font_size=12, text_color="black", + pts, list(ps.labels), font_size=round(12 * ts), text_color="black", font_family="times", shape=None, fill_shape=False, show_points=False, always_visible=True, name=(ps.name or "labels") + "-labels", @@ -505,7 +599,8 @@ def _render_pyvista(scene: Scene, show: bool = True): plotter.add_mesh(glyphs, color=_pyvista_color(g.color), name=g.name or None) if scene.title: - plotter.add_text(scene.title, position="upper_edge", font_size=10, color="black") + plotter.add_text(scene.title, position="upper_edge", font_size=round(10 * ts), + color="black") if scene.legend: plotter.add_legend( @@ -517,13 +612,14 @@ def _render_pyvista(scene: Scene, show: bool = True): plotter.add_scalar_bar( title=scene.colorbar.title, color="black", vertical=True, position_x=0.92, position_y=0.15, width=0.04, height=0.7, - title_font_size=16, label_font_size=14, font_family="times", + title_font_size=round(16 * ts), label_font_size=round(14 * ts), + font_family="times", fmt=scene.colorbar.fmt, ) if scene.show_axes: mins, maxs = scene.compute_bounds() - _draw_pyvista_axes(plotter, mins, maxs, scene.axis_labels) + _draw_pyvista_axes(plotter, mins, maxs, scene.axis_labels, text_scale=ts) plotter.view_isometric() plotter.camera.azimuth = 180 @@ -535,7 +631,50 @@ def _render_pyvista(scene: Scene, show: bool = True): return plotter -def _draw_pyvista_axes(plotter, mins, maxs, labels=("x (m)", "y (m)", "z (m)")) -> None: +def _nice_ticks(vmin, vmax, max_ticks: int = 8): + """Round tick values covering [vmin, vmax] using a 1/2/5 step sequence. + + Returns ``(values, decimals)`` where ``values`` are the tick positions + (at most ``max_ticks`` of them, aligned to multiples of the step) and + ``decimals`` is the number of fractional digits needed to print the step + exactly (0 for steps >= 1). + """ + span = float(vmax) - float(vmin) + if not np.isfinite(span) or span <= 0: + return [float(vmin)], 0 + raw_step = span / max(max_ticks - 1, 1) + magnitude = 10.0 ** np.floor(np.log10(raw_step)) + step = 10.0 * magnitude + for mult in (1.0, 2.0, 5.0): + if mult * magnitude >= raw_step - 1e-12 * magnitude: + step = mult * magnitude + break + first = np.ceil(float(vmin) / step - 1e-9) * step + count = int(np.floor((float(vmax) - first) / step + 1e-9)) + 1 + decimals = 0 if step >= 1.0 else int(np.ceil(-np.log10(step) - 1e-9)) + # Snap away float noise (0 + 3*0.2 -> 0.6000000000000001) so labels and + # equality checks see exact values. + snap = max(decimals + 2, 9) + values = [round(float(first + i * step), snap) for i in range(max(count, 0))] + return values, decimals + + +def _z_tick_values(zmin, zmax, max_ticks: int = 8): + """Nice ticks for the vertical axis with the base tick dropped. + + The z axis shares its base corner with the y axis, so a tick label at + (or near) the base collides with the y tick labels there. Mesh float + noise can put zmin fractionally below zero, so the base is excluded by + a span-relative threshold rather than an exact comparison. + """ + values, decimals = _nice_ticks(zmin, zmax, max_ticks) + span = float(zmax) - float(zmin) + values = [v for v in values if (v - float(zmin)) > 0.05 * span] + return values, decimals + + +def _draw_pyvista_axes(plotter, mins, maxs, labels=("x (m)", "y (m)", "z (m)"), + text_scale=1.0) -> None: """Draw x/y/z axes (lines, ticks, labels), working around a VTK StaticTriad z-axis rendering bug by drawing plain line meshes plus point labels.""" import pyvista as pv @@ -544,34 +683,97 @@ def _draw_pyvista_axes(plotter, mins, maxs, labels=("x (m)", "y (m)", "z (m)")) xmax_f, ymax_f, zmax_f = float(maxs[0]), float(maxs[1]), float(maxs[2]) span_xy = float(np.asarray(maxs - mins)[:2].max()) - def _draw_axis(p0, p1, ticks, label, tick_dir, font_size=14): + tick_len = 0.02 * span_xy + label_offset = 6.0 * tick_len + title_offset = 14.0 * tick_len + + def _draw_axis(p0, p1, ticks, label, tick_dir, font_size=16, decimals=0, + skip_first_label=False, label_dir=None, label_dist=None): + font_size = round(font_size * text_scale) plotter.add_mesh(pv.Line(p0, p1), color="black", line_width=3) - tick_len = 0.02 * span_xy - label_offset = 4.0 * tick_len pts, lbls = [], [] - for val, pos in ticks: + for i, (val, pos) in enumerate(ticks): t_end = [pos[j] + tick_dir[j] * tick_len for j in range(3)] plotter.add_mesh(pv.Line(pos, t_end), color="black", line_width=2) - pts.append([pos[j] + tick_dir[j] * label_offset for j in range(3)]) - lbls.append(f"{val:.0f}") + if skip_first_label and i == 0: + continue + if label_dir is None: + # default: label continues along the tick direction + pts.append([pos[j] + tick_dir[j] * label_offset for j in range(3)]) + else: + # label hangs off the tick tip in label_dir (e.g. straight + # down on screen), centred under the tick mark + pts.append([t_end[j] + label_dir[j] * label_dist for j in range(3)]) + lbls.append(f"{val:.{decimals}f}") if pts: + # Centre the text on its anchor: the default left/bottom + # justification grows the text box rightward, which visibly + # offsets labels from their tick lines (most on the x axis, + # whose ticks point down-right on screen). plotter.add_point_labels( np.array(pts), lbls, point_size=0, render_points_as_spheres=False, font_size=font_size, text_color="black", font_family="times", shape=None, fill_shape=False, show_points=False, margin=3, always_visible=True, + justification_horizontal="center", justification_vertical="center", ) mid = [0.5 * (p0[j] + p1[j]) for j in range(3)] - title_pos = [mid[j] + tick_dir[j] * tick_len * 8 for j in range(3)] + title_pos = [mid[j] + tick_dir[j] * title_offset for j in range(3)] plotter.add_point_labels( np.array([title_pos]), [label], point_size=0, render_points_as_spheres=False, show_points=False, font_size=font_size + 2, text_color="black", font_family="times", shape=None, fill_shape=False, bold=True, margin=3, always_visible=True, + justification_horizontal="center", justification_vertical="center", ) lx, ly, lz = labels - x_ticks = [(v, (v, ymin_f, zmin_f)) for v in np.linspace(xmin_f, xmax_f, 9)] - _draw_axis((xmin_f, ymin_f, zmin_f), (xmax_f, ymin_f, zmin_f), x_ticks, lx, (0, -1, 0)) - y_ticks = [(v, (xmin_f, v, zmin_f)) for v in np.linspace(ymin_f, ymax_f, 9)] - _draw_axis((xmin_f, ymin_f, zmin_f), (xmin_f, ymax_f, zmin_f), y_ticks, ly, (-1, 0, 0)) - z_ticks = [(v, (xmin_f, ymax_f, v)) for v in [zmin_f, zmax_f]] - _draw_axis((xmin_f, ymax_f, zmin_f), (xmin_f, ymax_f, zmax_f), z_ticks, lz, (-1, 0, 0)) + x_vals, x_dec = _nice_ticks(xmin_f, xmax_f) + y_vals, y_dec = _nice_ticks(ymin_f, ymax_f) + # When both tick sets start on the shared front corner with the same + # value, the two labels print on top of each other there. Replace them + # with one label hung straight below the corner: the world direction + # (-1, -1, 0) projects screen-vertical under the default camera (the + # sideways components of the two axis offsets cancel). + share_origin = ( + len(x_vals) > 0 and len(y_vals) > 0 + and (x_vals[0] - xmin_f) <= 0.05 * (xmax_f - xmin_f) + and (y_vals[0] - ymin_f) <= 0.05 * (ymax_f - ymin_f) + and x_vals[0] == y_vals[0] + ) + # World direction (-1, -1, 0) projects straight down on screen under the + # default camera; x/y tick labels hang below their tick tips along it so + # each number reads as centred under its tick mark. (1, -1, 0) projects + # pure screen-right and provides the small x-label nudge. + screen_down = (-1.0 / np.sqrt(2.0), -1.0 / np.sqrt(2.0), 0.0) + screen_right = (1.0 / np.sqrt(2.0), -1.0 / np.sqrt(2.0), 0.0) + label_drop = 3.5 * tick_len + x_label_vec = tuple(screen_down[j] * label_drop + screen_right[j] * 0.75 * tick_len + for j in range(3)) + x_ticks = [(v, (v, ymin_f, zmin_f)) for v in x_vals] + _draw_axis((xmin_f, ymin_f, zmin_f), (xmax_f, ymin_f, zmin_f), x_ticks, lx, (0, -1, 0), + decimals=x_dec, skip_first_label=share_origin, + label_dir=x_label_vec, label_dist=1.0) + y_ticks = [(v, (xmin_f, v, zmin_f)) for v in y_vals] + _draw_axis((xmin_f, ymin_f, zmin_f), (xmin_f, ymax_f, zmin_f), y_ticks, ly, (-1, 0, 0), + decimals=y_dec, skip_first_label=share_origin, + label_dir=screen_down, label_dist=label_drop) + if share_origin: + offd = label_offset / np.sqrt(2.0) + corner_pos = [xmin_f - offd, ymin_f - offd, zmin_f] + corner_dec = max(x_dec, y_dec) + plotter.add_point_labels( + np.array([corner_pos]), [f"{x_vals[0]:.{corner_dec}f}"], + point_size=0, render_points_as_spheres=False, + font_size=round(16 * text_scale), text_color="black", font_family="times", + shape=None, fill_shape=False, show_points=False, margin=3, always_visible=True, + justification_horizontal="center", justification_vertical="center", + ) + z_vals, z_dec = _z_tick_values(zmin_f, zmax_f) + z_ticks = [(v, (xmin_f, ymax_f, v)) for v in z_vals] + # The default camera (isometric, azimuth 180, elevation -10, parallel) + # projects pure -x downward on screen by 0.334 world-z units per offset + # unit (vtkCoordinate measurement). Compensate the z tick direction so + # ticks, labels and title extend horizontally on screen: diagonal tick + # stubs on the (screen-vertical) z axis read as kinks in the axis line. + z_tick_dir = (-1.0, 0.0, 0.334) + _draw_axis((xmin_f, ymax_f, zmin_f), (xmin_f, ymax_f, zmax_f), z_ticks, lz, z_tick_dir, + decimals=z_dec) diff --git a/tools/python/udvis/udbase_vis.py b/tools/python/udvis/udbase_vis.py index 07079e4a9..d60319a40 100644 --- a/tools/python/udvis/udbase_vis.py +++ b/tools/python/udvis/udbase_vis.py @@ -1,4 +1,4 @@ -"""Visualization facade for uDALES postprocessing. +"""Visualization facade for uDALES postprocessing. Provides the :class:`UDVis` class attached to each :class:`udbase.UDBase` instance as ``sim.vis``, offering geometry, field, and statistics @@ -14,6 +14,7 @@ import matplotlib.pyplot as plt from matplotlib.collections import PatchCollection +from matplotlib.colors import to_rgb from matplotlib.patches import Polygon as mplPolygon import numpy as np @@ -33,6 +34,69 @@ logger = logging.getLogger(__name__) +VEG_RGB = (34 / 255, 139 / 255, 34 / 255) # forest green + +# Faces whose centre sits within this height above z=0 still count as ground: +# STL float noise puts nominally flat ground at ~1e-9..1e-5 m, and a strict +# z > 0 test scatters ground triangles into the buildings group, painting +# most of the ground in the building colour (issue #337). No real building +# is a millimetre tall, so 1 mm cleanly separates noise from geometry. +GROUND_Z_TOL = 1e-3 # metres + + +def _is_building_faces(face_centers: np.ndarray) -> np.ndarray: + """Boolean mask of faces whose centre lies above the ground tolerance.""" + return np.asarray(face_centers)[:, 2] > GROUND_Z_TOL + + +def _as_rgb(color): + """Return an (r, g, b) tuple in 0..1 from a tuple or a matplotlib colour.""" + if isinstance(color, str): + return tuple(to_rgb(color)) + return tuple(float(c) for c in color[:3]) + + +# Triangulation of one box: 12 outward-wound triangles over the 8 corners +# ordered (x0,y0,z0),(x1,y0,z0),(x1,y1,z0),(x0,y1,z0), then the same at z1. +_BOX_TRIANGLES = np.array([ + [0, 2, 1], [0, 3, 2], # bottom + [4, 5, 6], [4, 6, 7], # top + [0, 1, 5], [0, 5, 4], # south (y0) + [1, 2, 6], [1, 6, 5], # east (x1) + [2, 3, 7], [2, 7, 6], # north (y1) + [3, 0, 4], [3, 4, 7], # west (x0) +], dtype=np.int64) + + +def _veg_voxel_mesh(points, xm, ym, zm, dx, dy, dzt): + """Triangulated boxes filling the grid cell of each vegetation point. + + The grid edge arrays hold the lower edge of every cell (length = number + of cells), so cell (i, j, k) spans [xm[i], xm[i] + dx] x [ym[j], ym[j] + dy] + x [zm[k], zm[k] + dzt[k]] (dzt carries stretched-grid spacings). + + Returns ``(vertices, faces)`` with shapes (8n, 3) and (12n, 3). + """ + pts = np.asarray(points, dtype=int) + if pts.size == 0: + return np.empty((0, 3), dtype=float), np.empty((0, 3), dtype=np.int64) + x0 = np.asarray(xm, dtype=float)[pts[:, 0]] + y0 = np.asarray(ym, dtype=float)[pts[:, 1]] + z0 = np.asarray(zm, dtype=float)[pts[:, 2]] + x1 = x0 + float(dx) + y1 = y0 + float(dy) + z1 = z0 + np.asarray(dzt, dtype=float)[pts[:, 2]] + n = len(pts) + corners = np.empty((n, 8, 3), dtype=float) + for idx, (cx, cy, cz) in enumerate( + [(x0, y0, z0), (x1, y0, z0), (x1, y1, z0), (x0, y1, z0), + (x0, y0, z1), (x1, y0, z1), (x1, y1, z1), (x0, y1, z1)]): + corners[:, idx, 0] = cx + corners[:, idx, 1] = cy + corners[:, idx, 2] = cz + faces = _BOX_TRIANGLES[None, :, :] + 8 * np.arange(n, dtype=np.int64)[:, None, None] + return corners.reshape(-1, 3), faces.reshape(-1, 3) + class UDVis: """ @@ -166,17 +230,21 @@ def show_geometry( faces = np.asarray(stl.faces, dtype=int) face_centers = np.asarray(stl.triangles_center, dtype=float) face_normals = np.asarray(stl.face_normals, dtype=float) - is_building = face_centers[:, 2] > 0 + is_building = _is_building_faces(face_centers) - lighting = "flat" if (color_buildings and show_ground) else None + # No "flat" lighting preset here: these meshes carry uniform solid + # colours, so ambient-only lighting collapses the Plotly view into an + # unshaded silhouette (issue #337). Default lighting shades them per + # face, matching the PyVista backend. The outline plots keep "flat" + # because their shading is pre-baked into two-tone face colours. meshes = [] if color_buildings: if show_ground and np.any(~is_building): meshes.append(MeshPrimitive(vertices, faces[~is_building], - solid_color=GROUND_RGB, name="ground", lighting=lighting)) + solid_color=GROUND_RGB, name="ground")) if np.any(is_building): meshes.append(MeshPrimitive(vertices, faces[is_building], - solid_color=BUILDING_RGB, name="buildings", lighting=lighting)) + solid_color=BUILDING_RGB, name="buildings")) else: selected_faces = faces if show_ground else faces[is_building] meshes.append(MeshPrimitive(vertices, selected_faces, solid_color=GROUND_RGB, name="geometry")) @@ -191,7 +259,13 @@ def show_geometry( edge_faces = faces if show_ground else faces[is_building] segments = np.asarray(self._collect_mesh_edges(edge_faces), dtype=int) if len(segments): - scene.lines.append(LineSet(vertices, segments, color="black", width=2, lift_ground=True)) + # PyVista: thin and semi-transparent, so dense wireframes read + # as shading instead of a black mass (issue #337). Plotly keeps + # its original bold opaque lines - its WebGL renderer turns + # thin translucent lines into stipple. + scene.lines.append(LineSet(vertices, segments, color="black", width=1.0, + opacity=0.4, lift_ground=True, + width_plotly=2.0, opacity_plotly=1.0)) if plot_quiver: scene.glyphs.append(GlyphSet(points=face_centers, vectors=face_normals, @@ -231,6 +305,20 @@ def show_geometry_outline( ------- plotly.graph_objects.Figure or pyvista.Plotter or None """ + scene = self._outline_scene(angle_threshold=angle_threshold, + show_ground=show_ground, + color_buildings=color_buildings) + if scene is None: + return None + return render_scene(scene, backend=self._resolve_backend(backend), show=show) + + def _outline_scene(self, angle_threshold: float = 45.0, show_ground: bool = True, + color_buildings: bool = False, line_width: float = 2.0): + """Geometry Scene with outline edges, shared by the outline plots. + + Returns ``None`` (after a warning) when the mesh yields no outline + edges; raises when no geometry is loaded. + """ if self.geom is None or getattr(self.geom, "stl", None) is None: raise ValueError("No geometry loaded. Cannot visualize.") @@ -244,7 +332,7 @@ def show_geometry_outline( faces = np.asarray(stl.faces, dtype=int) face_centers = np.asarray(stl.triangles_center, dtype=float) face_normals = np.asarray(stl.face_normals, dtype=float) - is_building = face_centers[:, 2] > 0 + is_building = _is_building_faces(face_centers) meshes = [] if color_buildings: @@ -264,13 +352,13 @@ def show_geometry_outline( name="geometry", lighting="flat")) segments = np.asarray(outline_edges, dtype=int) - scene = Scene( + return Scene( meshes=meshes, - lines=[LineSet(vertices, segments, color="black", width=2, lift_ground=True)], + lines=[LineSet(vertices, segments, color="black", width=line_width, + lift_ground=True)], title=f"Geometry Outline ({len(outline_edges)} edges)", bounds=(vertices.min(axis=0), vertices.max(axis=0)), ) - return render_scene(scene, backend=self._resolve_backend(backend), show=show) def _base_overlay_scene(self, geom, title: Optional[str] = None) -> Scene: """Scene with the geometry as a grey base mesh plus its outline edges. @@ -289,12 +377,21 @@ def _base_overlay_scene(self, geom, title: Optional[str] = None) -> Scene: ) segments = np.asarray(self._collect_mesh_edges(faces), dtype=int) if len(segments): - scene.lines.append(LineSet(vertices, segments, color="black", width=1.5, lift_ground=True)) + # Thin, semi-transparent full wireframe in PyVista; original bold + # opaque lines in Plotly - see show_geometry for the rationale + # (issue #337). + scene.lines.append(LineSet(vertices, segments, color="black", width=1.0, + opacity=0.4, lift_ground=True, + width_plotly=1.5, opacity_plotly=1.0)) return scene def plot_veg(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, - backend: Optional[str] = None): - """Plot vegetation points on top of the geometry. + backend: Optional[str] = None, color=VEG_RGB, opacity: float = 1.0): + """Plot vegetation on top of the geometry. + + Every vegetation point is drawn as an opaque box filling its grid + cell, so canopy volumes read at their true size instead of as + near-invisible dots (issue #337). Parameters ---------- @@ -305,6 +402,10 @@ def plot_veg(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, unlike sibling methods), return the figure/plotter. backend : {"plotly", "pyvista"}, optional Rendering backend; defaults to this UDVis instance's backend. + color : matplotlib colour or (r, g, b) tuple, optional + Vegetation colour; defaults to forest green. + opacity : float, default=1.0 + Vegetation opacity in 0..1. Returns ------- @@ -312,6 +413,28 @@ def plot_veg(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, """ if not self.sim._lfgeom or self.sim.geom is None: return self._missing_plot_data("Geometry data not found for plot_veg") + points = self._resolve_veg_points(veg, "plot_veg") + if points is None: + return None + + verts, faces = _veg_voxel_mesh( + points, self.sim.xm, self.sim.ym, self.sim.zm, + self.sim.dx, self.sim.dy, self.sim.dzt) + + scene = self._base_overlay_scene( + self.sim.geom, title=f"Geometry with Vegetation ({len(points)} points)") + scene.meshes.append(MeshPrimitive( + verts, faces, solid_color=_as_rgb(color), opacity=float(opacity), + name="vegetation")) + return render_scene(scene, backend=self._resolve_backend(backend), show=show) + + def _resolve_veg_points(self, veg: Optional[Dict[str, Any]], caller: str): + """Vegetation points from the argument or the case, subsampled. + + Returns an (n, 3) int array, or ``None`` after a single warning when + the case has no usable vegetation data (shared by the vegetation + plots so they degrade identically on cases without vegetation). + """ if veg is None: if not hasattr(self.sim, "veg") or self.sim.veg is None: try: @@ -329,17 +452,60 @@ def plot_veg(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, if len(points) > max_points: rng = np.random.default_rng(0) points = points[rng.choice(len(points), size=max_points, replace=False)] - logger.info("plot_veg: showing %d of %d points", max_points, len(veg['points'])) + logger.info("%s: showing %d of %d points", caller, max_points, len(veg['points'])) + return points - xs = self.sim.xt[points[:, 0].astype(int)] - ys = self.sim.yt[points[:, 1].astype(int)] - zs = self.sim.zt[points[:, 2].astype(int)] + def plot_veg_outline(self, veg: Optional[Dict[str, Any]] = None, show: bool = False, + backend: Optional[str] = None, color=VEG_RGB, + opacity: float = 1.0, line_width: float = 2.0, + angle_threshold: float = 45.0): + """Plot vegetation voxels over outline-style geometry. - scene = self._base_overlay_scene( - self.sim.geom, title=f"Geometry with Vegetation ({len(points)} points)") - scene.points.append(PointSet( - np.column_stack([xs, ys, zs]), - color="rgb(34,139,34)", size=2, opacity=0.2, name="vegetation")) + Same vegetation rendering as :meth:`plot_veg` (opaque boxes filling + their grid cells), but on the two-tone outline base of + :meth:`show_geometry_outline` instead of the full wireframe - the + most readable vegetation view on fine meshes (issue #337). + + Parameters + ---------- + veg : dict, optional + Vegetation data; loaded from the case when omitted. + show : bool, default=False + If True, display immediately and return None. If False (default, + like the other overlay plots), return the figure/plotter. + backend : {"plotly", "pyvista"}, optional + Rendering backend; defaults to this UDVis instance's backend. + color : matplotlib colour or (r, g, b) tuple, optional + Vegetation colour; defaults to forest green. + opacity : float, default=1.0 + Vegetation opacity in 0..1. + line_width : float, default=2.0 + Outline edge width. + angle_threshold : float, default=45.0 + Angle threshold used to detect outline edges. + + Returns + ------- + plotly.graph_objects.Figure or pyvista.Plotter or None + """ + if self.sim is None or not self.sim._lfgeom or self.sim.geom is None: + return self._missing_plot_data("Geometry data not found for plot_veg_outline") + points = self._resolve_veg_points(veg, "plot_veg_outline") + if points is None: + return None + + scene = self._outline_scene(angle_threshold=angle_threshold, + line_width=line_width) + if scene is None: + return None + scene.title = f"Geometry Outline with Vegetation ({len(points)} points)" + + verts, faces = _veg_voxel_mesh( + points, self.sim.xm, self.sim.ym, self.sim.zm, + self.sim.dx, self.sim.dy, self.sim.dzt) + scene.meshes.append(MeshPrimitive( + verts, faces, solid_color=_as_rgb(color), opacity=float(opacity), + name="vegetation")) return render_scene(scene, backend=self._resolve_backend(backend), show=show) @staticmethod