update

Author: geohang

PyHydroGeophysX/qt_apps/artifact_renderers.py

@@ -20,13 +20,15 @@ def select_renderer(
         fmt = Path(str(artifact.get("path") or "")).suffix.lower().lstrip(".")
     if "model_bundle" in kind or kind in {"pygimli_bundle", "mesh_model_bundle"}:
         return "mesh_bundle"
-    # A bare PyGIMLi mesh is the primary product of the 3-D mesh builder, so it
-    # needs its own viewer rather than falling through to "no renderer".
-    if fmt == "bms" or kind in {"mesh", "pygimli_mesh"}:
+    # A recognized file format is stronger evidence than a broad semantic kind.
+    # In particular, ``kind=volume`` with ``format=npy`` is a NumPy stack, not a
+    # VTK file, and a ``figure_*`` data artifact must not be sent to an image
+    # decoder merely because its kind contains the word "figure".
+    if fmt == "bms":
         return "mesh"
-    if fmt in _VTK_FORMATS or kind in {"velocity_model", "volume", "vtk"}:
+    if fmt in _VTK_FORMATS:
         return "vtk"
-    if fmt in _IMAGE_FORMATS or "figure" in kind or kind == "image":
+    if fmt in _IMAGE_FORMATS:
         return "image"
     if fmt in {"npy", "npz"}:
         shape = tuple(array_shape or artifact.get("shape") or ())
@@ -41,6 +43,14 @@ def select_renderer(
         return "table"
     if fmt == "json":
         return "json"
+    # Kind-only fallbacks cover virtual/in-record artifacts and files without a
+    # useful extension. A bare PyGIMLi mesh is the primary mesh-builder product.
+    if kind in {"mesh", "pygimli_mesh"}:
+        return "mesh"
+    if kind in {"velocity_model", "volume", "vtk"}:
+        return "vtk"
+    if "figure" in kind or kind == "image":
+        return "image"
     return "file"
 
 

PyHydroGeophysX/qt_apps/modules/em_input_format.md

@@ -106,6 +106,16 @@ A header-less file is read by column order: one column = position; two columns =
 
 ## Notes
 
+- A line inversion runs the per-sounding work on threads: reading the stations,
+  building their forward operators, and every forward and Jacobian evaluation
+  inside the coupled solve. The thread count comes from the machine; set
+  `parallel_workers` through `set_params` to pin it. Each sounding owns its
+  forward operator and the workers only read the shared model, so the models
+  come back bit for bit identical to a serial run, which the test suite checks
+  over repeated solves. The coupled solve itself measures 9.5 times faster on 20
+  threads; end to end the gain is smaller, because building the forward
+  operators is a fixed cost that threads do not remove (roughly 1.4 times on a
+  short run, 3.9 times on one that iterates 25 times).
 - **Auto-calibrate** (checkbox, on by default) estimates the data-scale
   calibration from the data before inverting. Leave it on for normalized airborne
   data (e.g. moment-normalized dB/dt); it returns ~1 for data already in the

PyHydroGeophysX/qt_apps/modules/ert_processing.py

@@ -157,14 +157,37 @@ def __init__(self, state: Any, log: LogFn, parent=None) -> None:
         self._plot.addItem(self._sel_scatter)
         self._plot.scene().sigMouseClicked.connect(self._on_click)
 
-        self._pseudo_widget = pg.PlotWidget()
+        # Use a GraphicsLayout rather than a bare PlotWidget so the colour scale
+        # is part of the pseudosection.  The plotted quantity is log10(rhoa), but
+        # the colour-bar axis formats those exponents back into physical ohm-m
+        # values; readers should not have to infer resistivity from colour alone.
+        self._pseudo_widget = pg.GraphicsLayoutWidget()
         self._pseudo_widget.setBackground("w")
-        self._pseudo_widget.showGrid(x=True, y=True, alpha=0.25)
-        self._pseudo_widget.setLabel("bottom", "x (m)")
-        self._pseudo_widget.setLabel("left", "pseudo-depth (m)")
-        self._pseudo_plot = self._pseudo_widget.getPlotItem()
-        self._pseudo_scatter = pg.ScatterPlotItem(size=11)
+        self._pseudo_plot = self._pseudo_widget.addPlot(row=0, col=0)
+        self._pseudo_plot.showGrid(x=True, y=True, alpha=0.25)
+        self._pseudo_plot.setLabel("bottom", "x (m)")
+        self._pseudo_plot.setLabel("left", "pseudo-depth (m, positive down)")
+        self._pseudo_plot.invertY(True)
+        self._pseudo_scatter = pg.ScatterPlotItem(size=11, hoverable=True)
         self._pseudo_plot.addItem(self._pseudo_scatter)
+        self._pseudo_scatter.sigHovered.connect(self._on_pseudo_hover)
+        self._pseudo_colorbar = pg.ColorBarItem(
+            values=(0.0, 1.0),
+            width=18,
+            colorMap=self._cmap,
+            label="Apparent resistivity (Ω·m)",
+            interactive=False,
+            colorMapMenu=False,
+        )
+        self._pseudo_colorbar.axis.setLogMode(True)
+        self._pseudo_widget.addItem(self._pseudo_colorbar, row=0, col=1)
+        self._pseudo_readout = self._pseudo_widget.addLabel(
+            "Hover a measurement to read apparent resistivity.",
+            row=1,
+            col=0,
+            colspan=2,
+            justify="left",
+        )
 
         self._model_view = MeshResultView()
         # The "Resistivity model" tab shows the single inversion OR any time step
@@ -2226,26 +2249,61 @@ def _draw_pseudosection(self) -> None:
         if not self._pseudo:
             self._pseudo_scatter.setData([])
             self._pseudo_plot.setTitle("")
+            self._pseudo_colorbar.setVisible(False)
+            self._pseudo_readout.setText("No apparent-resistivity measurements loaded.")
             return
         arr = np.asarray(self._pseudo, dtype=float)
         mid, depth, rhoa = arr[:, 0], arr[:, 1], arr[:, 2]
         valid = np.isfinite(rhoa) & (rhoa > 0)
         mid, depth, rhoa = mid[valid], depth[valid], rhoa[valid]
         if rhoa.size == 0:
             self._pseudo_scatter.setData([])
+            self._pseudo_plot.setTitle("No positive finite apparent resistivity to display")
+            self._pseudo_colorbar.setVisible(False)
+            self._pseudo_readout.setText(
+                "All apparent-resistivity values are missing, non-finite, or non-positive."
+            )
             return
         log_rhoa = np.log10(rhoa)
         lo, hi = np.percentile(log_rhoa, [3, 97])
-        rng = hi - lo if hi > lo else 1.0
+        if hi <= lo:
+            lo, hi = float(lo) - 0.5, float(hi) + 0.5
+        rng = hi - lo
         norm = np.clip((log_rhoa - lo) / rng, 0.0, 1.0)
         lut = self._cmap.map(norm, mode="byte")
         spots = [
-            {"pos": (float(mid[i]), -float(depth[i])),
-             "brush": pg.mkBrush(int(lut[i, 0]), int(lut[i, 1]), int(lut[i, 2])), "size": 11}
+            {"pos": (float(mid[i]), float(depth[i])),
+             "data": (float(mid[i]), float(depth[i]), float(rhoa[i])),
+             "brush": pg.mkBrush(int(lut[i, 0]), int(lut[i, 1]), int(lut[i, 2])),
+             "size": 11}
             for i in range(mid.size)
         ]
         self._pseudo_scatter.setData(spots)
-        self._pseudo_plot.setTitle(f"Apparent resistivity (Ω·m): {rhoa.min():.0f} – {rhoa.max():.0f}  (n={rhoa.size})")
+        self._pseudo_colorbar.setLevels((float(lo), float(hi)))
+        self._pseudo_colorbar.setVisible(True)
+        self._pseudo_readout.setText(
+            "Hover a measurement to read x, pseudo-depth, and apparent resistivity."
+        )
+        self._pseudo_plot.setTitle(
+            f"Apparent resistivity: {rhoa.min():.3g} – {rhoa.max():.3g} Ω·m  "
+            f"(n={rhoa.size})"
+        )
+
+    def _on_pseudo_hover(self, _item, points, _event) -> None:
+        """Report the physical value behind a pseudosection colour."""
+        if not points:
+            self._pseudo_readout.setText(
+                "Hover a measurement to read x, pseudo-depth, and apparent resistivity."
+            )
+            return
+        data = points[0].data()
+        if not data or len(data) != 3:
+            return
+        x, depth, rhoa = data
+        self._pseudo_readout.setText(
+            f"x = {float(x):.3g} m    ·    pseudo-depth = {float(depth):.3g} m"
+            f"    ·    ρa = {float(rhoa):.4g} Ω·m"
+        )
 
     # -- interaction ---------------------------------------------------------
     def _nearest(self, x: float, z: float) -> Optional[int]:

PyHydroGeophysX/qt_apps/modules/hydro_geophysics.py

@@ -959,7 +959,22 @@ def _update_display(self) -> None:
             return
         if arr2d is None:
             return
-        self._map.set_array(arr2d)
+        labels = {
+            "Water content": "Volumetric water content",
+            "Porosity": "Porosity",
+            "Top": "Top elevation / depth (m)",
+            "Bottom": "Bottom elevation / depth (m)",
+        }
+        title = var
+        if var in {"Water content", "Porosity", "Bottom"}:
+            title += f" — layer {layer}"
+        if var == "Water content" and self._wc is not None and self._wc.ndim == 4:
+            title += f", snapshot {self._snapshot.value()}"
+        self._map.set_array(
+            arr2d,
+            value_label=labels.get(var, var),
+            title=title,
+        )
         if self._point1 and self._point2:
             self._map.set_profile_points(self._point1, self._point2)
 
@@ -1005,7 +1020,13 @@ def _on_preview_ready(self, profile: dict) -> None:
         self._profile = profile
         try:
             wc = np.asarray(profile["water_content_profile"], dtype=float)
-            self._preview.set_array(wc)
+            self._preview.set_array(
+                wc,
+                x_label="Profile sample",
+                y_label="Layer",
+                value_label="Volumetric water content",
+                title="Water-content profile preview",
+            )
         except Exception as exc:  # noqa: BLE001
             self.log(f"Profile preview failed: {exc}", "warn")
 

PyHydroGeophysX/qt_apps/modules/model_viewer.py

@@ -127,6 +127,41 @@ def _artifact_label(artifact: Dict[str, Any], *, missing: bool = False) -> str:
     return f"{label} (missing)" if missing else label
 
 
+def _artifact_plot_options(artifact: Dict[str, Any], path: Path) -> tuple[str, str, bool]:
+    """Return ``(title, value_label, log_scale)`` for numeric previews.
+
+    Result artifacts already carry field metadata in several workflows.  Keep
+    that meaning when the file reaches the generic viewer instead of reducing
+    every grid to anonymous rows, columns, and an unlabeled linear colour bar.
+    """
+    metadata = artifact.get("metadata")
+    metadata = dict(metadata) if isinstance(metadata, dict) else {}
+    kind = str(artifact.get("kind") or "")
+    field = str(
+        metadata.get("label")
+        or metadata.get("field")
+        or artifact.get("label")
+        or (_pretty_kind(kind) if kind else "Value")
+    ).strip()
+    units = str(metadata.get("units") or "").strip()
+    tokens = " ".join((kind, path.stem, field)).lower()
+    is_resistivity = any(token in tokens for token in ("resistiv", "rhoa", "rho_a"))
+    if not units and is_resistivity:
+        units = "Ω·m"
+    value_label = field or "Value"
+    if units and units.lower() not in value_label.lower():
+        value_label = f"{value_label} ({units})"
+    explicit_log = metadata.get("log_scale")
+    if isinstance(explicit_log, str):
+        log_scale = explicit_log.strip().lower() in {"1", "true", "yes", "log", "log10"}
+    elif explicit_log is None:
+        log_scale = is_resistivity
+    else:
+        log_scale = bool(explicit_log)
+    title = str(metadata.get("title") or path.stem.replace("_", " ").strip()).strip()
+    return title, value_label, log_scale
+
+
 def _run_title(record: "RunRecord") -> str:
     """Show the user's own name for a run, or a short stable one.
 
@@ -789,11 +824,13 @@ def _render_selected_artifact(self, _index: int) -> None:
         renderer = select_renderer(artifact)
         try:
             if renderer in {"array", "array_stack", "curve"} and path.suffix.lower() in {".npy", ".npz"}:
-                self._render_numpy(path, renderer)
+                self._render_numpy(path, renderer, artifact)
             elif renderer == "curve":
                 self._render_curve_file(path)
             elif renderer == "image":
-                view = ZoomableImageView(); view.set_image_file(path)
+                view = ZoomableImageView()
+                if not view.set_image_file(path):
+                    raise ValueError("the image decoder could not read this file")
                 self._replace_visual(view)
             elif renderer == "vtk":
                 from PyHydroGeophysX.qt_apps.widgets.model3d_view import VTKVolumeView
@@ -824,7 +861,14 @@ def _replace_visual(
         self._visual_resources.extend(resources or [])
         self._visual_layout.addWidget(widget)
 
-    def _render_numpy(self, path: Path, renderer: str) -> None:
+    def _render_numpy(
+        self,
+        path: Path,
+        renderer: str,
+        artifact: Optional[Dict[str, Any]] = None,
+    ) -> None:
+        artifact = dict(artifact or {})
+        title, value_label, log_scale = _artifact_plot_options(artifact, path)
         loaded = np.load(
             path,
             allow_pickle=False,
@@ -839,20 +883,43 @@ def _render_numpy(self, path: Path, renderer: str) -> None:
                 ]
                 if not names:
                     raise ValueError("NPZ contains no numeric arrays.")
-                array = np.array(loaded[names[0]], copy=True)
+                metadata = artifact.get("metadata")
+                metadata = dict(metadata) if isinstance(metadata, dict) else {}
+                requested = str(metadata.get("array_key") or "")
+                chosen = requested if requested in names else names[0]
+                array = np.array(loaded[chosen], copy=True)
                 loaded.close()
             else:
                 array = loaded
-            renderer = select_renderer({"path": str(path)}, array.shape)
+            # Keep semantic dispatch from the artifact (notably a 2-D curve
+            # table).  Shape only upgrades/downgrades the generic array modes.
+            if renderer == "array" and array.ndim >= 3:
+                renderer = "array_stack"
+            elif renderer == "array_stack" and array.ndim < 3:
+                renderer = "curve" if array.ndim == 1 else "array"
             if renderer == "curve" or array.ndim == 1:
-                values = np.array(array, copy=True).ravel()
-                view = CurveViewer(); view.add_curve(
-                    np.arange(values.size), values, path.stem
-                )
+                values = np.array(array, copy=True)
+                view = CurveViewer()
+                if values.ndim == 2 and min(values.shape) >= 2:
+                    # A compact 2 x N / 3 x N array normally stores x and one
+                    # or two series by row; a tall N x K array stores columns.
+                    if values.shape[0] <= 4 and values.shape[1] > 2 * values.shape[0]:
+                        values = values.T
+                    x = values[:, 0]
+                    for column in range(1, values.shape[1]):
+                        view.add_curve(x, values[:, column], f"{value_label} {column}")
+                else:
+                    values = values.ravel()
+                    view.add_curve(np.arange(values.size), values, value_label)
                 self._replace_visual(view); return
             if array.ndim == 2:
                 values = np.array(array, copy=True)
-                view = ArrayViewer(); view.set_array(values)
+                view = ArrayViewer(); view.set_array(
+                    values,
+                    log=log_scale,
+                    value_label=value_label,
+                    title=title,
+                )
                 self._replace_visual(view); return
             if array.ndim >= 3:
                 host = QWidget(); layout = QVBoxLayout(host)
@@ -861,17 +928,29 @@ def _render_numpy(self, path: Path, renderer: str) -> None:
                 controls.addWidget(QLabel("Step / slice:")); controls.addWidget(step); controls.addStretch(1)
                 view = ArrayViewer()
                 sample = np.asarray(array).ravel()[::max(1, array.size // 250_000)]
+                if log_scale:
+                    with np.errstate(divide="ignore", invalid="ignore"):
+                        sample = np.log10(np.where(sample > 0, sample, np.nan))
                 finite = sample[np.isfinite(sample)]
-                limits = np.percentile(finite, [2, 98]) if finite.size else None
+                limits = None
+                if finite.size:
+                    lo, hi = np.percentile(finite, [2, 98])
+                    if hi <= lo:
+                        half_span = 0.5 if log_scale else max(abs(float(lo)) * 0.05, 0.5)
+                        lo, hi = float(lo) - half_span, float(hi) + half_span
+                    limits = (float(lo), float(hi))
                 def show_slice(value: int) -> None:
                     # Copy only the visible slice: the widget never owns a view
                     # of the file mapping, so Delete Run can close it first.
                     view.set_array(
                         np.array(array[value], copy=True),
                         autoscale=limits is None,
+                        log=log_scale,
+                        value_label=value_label,
+                        title=f"{title} — slice {value + 1}",
                     )
-                    if limits is not None and limits[1] > limits[0]:
-                        view.set_levels(float(limits[0]), float(limits[1]))
+                    if limits is not None:
+                        view.set_levels(*limits)
                 step.valueChanged.connect(show_slice); show_slice(0)
                 layout.addLayout(controls); layout.addWidget(view, stretch=1)
                 self._replace_visual(host, resources=[loaded])