canvod.viz API Reference¶

class HemisphereVisualizer:
    """Unified hemisphere visualizer combining 2D and 3D capabilities.

    Provides consistent API for both publication-quality matplotlib plots
    and interactive plotly visualizations. Handles styling coordination
    between different rendering backends.

    Parameters
    ----------
    grid : HemiGrid
        Hemisphere grid to visualize

    Examples
    --------
    Create both 2D and 3D visualizations::

        from canvod.grids import create_hemigrid
        from canvod.viz import HemisphereVisualizer

        grid = create_hemigrid(grid_type='equal_area', angular_resolution=10.0)
        viz = HemisphereVisualizer(grid)

        # Publication-quality 2D plot
        fig_2d, ax_2d = viz.plot_2d(
            data=vod_data,
            title="VOD Distribution",
            save_path="publication.png"
        )

        # Interactive 3D plot
        fig_3d = viz.plot_3d(
            data=vod_data,
            title="Interactive VOD Explorer"
        )
        fig_3d.show()

    Switch styles easily::

        # Publication style
        pub_style = create_publication_style()
        viz.set_style(pub_style)
        fig, ax = viz.plot_2d(data=vod_data)

        # Interactive style
        int_style = create_interactive_style(dark_mode=True)
        viz.set_style(int_style)
        fig = viz.plot_3d(data=vod_data)

    """

    def __init__(self, grid: HemiGrid) -> None:
        """Initialize unified visualizer.

        Parameters
        ----------
        grid : HemiGrid
            Hemisphere grid to visualize

        """
        self.grid = grid

        # Initialize specialized visualizers
        self.viz_2d = HemisphereVisualizer2D(grid)
        self.viz_3d = HemisphereVisualizer3D(grid)

        # Default styling
        self.style = PlotStyle()

    def set_style(self, style: PlotStyle) -> None:
        """Set unified styling for both 2D and 3D plots.

        Parameters
        ----------
        style : PlotStyle
            Styling configuration

        Examples
        --------
        >>> pub_style = create_publication_style()
        >>> viz.set_style(pub_style)
        >>> fig, ax = viz.plot_2d(data=vod_data)

        """
        self.style = style

    def plot_2d(
        self,
        data: np.ndarray | None = None,
        title: str | None = None,
        ax: Axes | None = None,
        save_path: Path | str | None = None,
        style: PolarPlotStyle | None = None,
        **kwargs: Any,
    ) -> tuple[Figure, Axes]:
        """Create 2D publication-quality plot.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, optional
            Plot title
        ax : matplotlib.axes.Axes, optional
            Existing axes to plot on
        save_path : Path or str, optional
            Save figure to this path
        style : PolarPlotStyle, optional
            Override default 2D style
        **kwargs
            Additional styling parameters

        Returns
        -------
        fig : matplotlib.figure.Figure
            Figure object
        ax : matplotlib.axes.Axes
            Polar axes with plot

        Examples
        --------
        >>> fig, ax = viz.plot_2d(
        ...     data=vod_data,
        ...     title="VOD Distribution",
        ...     cmap='plasma',
        ...     save_path="output.png",
        ...     dpi=300
        ... )

        """
        if style is None:
            style = self.style.to_polar_style()

        if title:
            style.title = title

        return self.viz_2d.plot_grid_patches(
            data=data, style=style, ax=ax, save_path=save_path, **kwargs
        )

    def plot_3d(
        self,
        data: np.ndarray | None = None,
        title: str | None = None,
        style: PlotStyle | None = None,
        **kwargs: Any,
    ) -> go.Figure:
        """Create 3D interactive plot.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, optional
            Plot title
        style : PlotStyle, optional
            Override default 3D style
        **kwargs
            Additional plotly parameters

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive 3D figure

        Examples
        --------
        >>> fig = viz.plot_3d(
        ...     data=vod_data,
        ...     title="Interactive VOD",
        ...     opacity=0.9,
        ...     width=1000,
        ...     height=800
        ... )
        >>> fig.show()
        >>> fig.write_html("interactive.html")

        """
        if style is None:
            style = self.style

        return self.viz_3d.plot_hemisphere_surface(
            data=data, style=style, title=title, **kwargs
        )

    def plot_3d_mesh(
        self,
        data: np.ndarray | None = None,
        title: str | None = None,
        **kwargs: Any,
    ) -> go.Figure:
        """Create 3D mesh plot showing cell boundaries.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, optional
            Plot title
        **kwargs
            Additional plotly parameters

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive mesh figure

        Examples
        --------
        >>> fig = viz.plot_3d_mesh(
        ...     data=vod_data,
        ...     title="VOD Mesh View",
        ...     opacity=0.7
        ... )

        """
        return self.viz_3d.plot_cell_mesh(data=data, title=title, **kwargs)

    def create_comparison_plot(
        self,
        data: np.ndarray | None = None,
        title_2d: str = "2D Polar View",
        title_3d: str = "3D Hemisphere View",
        save_2d: Path | str | None = None,
        save_3d: Path | str | None = None,
    ) -> tuple[tuple[Figure, Axes], go.Figure]:
        """Create both 2D and 3D plots for comparison.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title_2d : str, default "2D Polar View"
            Title for 2D plot
        title_3d : str, default "3D Hemisphere View"
            Title for 3D plot
        save_2d : Path or str, optional
            Save 2D figure to this path
        save_3d : Path or str, optional
            Save 3D figure to this path (HTML)

        Returns
        -------
        plot_2d : tuple of Figure and Axes
            2D matplotlib plot
        plot_3d : plotly.graph_objects.Figure
            3D plotly plot

        Examples
        --------
        >>> (fig_2d, ax_2d), fig_3d = viz.create_comparison_plot(
        ...     data=vod_data,
        ...     save_2d="comparison_2d.png",
        ...     save_3d="comparison_3d.html"
        ... )
        >>> plt.show()  # Show 2D
        >>> fig_3d.show()  # Show 3D

        """
        # Create 2D plot
        fig_2d, ax_2d = self.plot_2d(data=data, title=title_2d, save_path=save_2d)

        # Create 3D plot
        fig_3d = self.plot_3d(data=data, title=title_3d)

        # Save 3D if requested
        if save_3d:
            save_3d = Path(save_3d)
            save_3d.parent.mkdir(parents=True, exist_ok=True)
            fig_3d.write_html(str(save_3d))

        return (fig_2d, ax_2d), fig_3d

    def create_publication_figure(
        self,
        data: np.ndarray | None = None,
        title: str = "Hemispherical Data Distribution",
        save_path: Path | str | None = None,
        dpi: int = 300,
        **kwargs: Any,
    ) -> tuple[Figure, Axes]:
        """Create publication-ready figure with optimal styling.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, default "Hemispherical Data Distribution"
            Plot title
        save_path : Path or str, optional
            Save figure to this path
        dpi : int, default 300
            Resolution in dots per inch
        **kwargs
            Additional styling parameters

        Returns
        -------
        fig : matplotlib.figure.Figure
            Publication-ready figure
        ax : matplotlib.axes.Axes
            Styled polar axes

        Examples
        --------
        >>> fig, ax = viz.create_publication_figure(
        ...     data=vod_data,
        ...     title="VOD Distribution Over Rosalia Site",
        ...     save_path="paper_figure_3.png",
        ...     dpi=600
        ... )

        """
        # Use publication style and convert to PolarPlotStyle
        pub_plot_style = create_publication_style()
        polar_style = pub_plot_style.to_polar_style()
        polar_style.title = title
        polar_style.dpi = dpi

        # Override with kwargs
        for key, value in kwargs.items():
            if hasattr(polar_style, key):
                setattr(polar_style, key, value)

        return self.plot_2d(data=data, style=polar_style, save_path=save_path)

    def create_interactive_explorer(
        self,
        data: np.ndarray | None = None,
        title: str = "Interactive Data Explorer",
        dark_mode: bool = True,
        save_html: Path | str | None = None,
    ) -> go.Figure:
        """Create interactive explorer with optimal settings.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, default "Interactive Data Explorer"
            Plot title
        dark_mode : bool, default True
            Use dark theme
        save_html : Path or str, optional
            Save HTML to this path

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive explorer figure

        Examples
        --------
        >>> fig = viz.create_interactive_explorer(
        ...     data=vod_data,
        ...     title="VOD Explorer",
        ...     dark_mode=True,
        ...     save_html="explorer.html"
        ... )
        >>> fig.show()

        """
        # Use interactive style
        int_style = create_interactive_style(dark_mode=dark_mode)

        fig = self.plot_3d(data=data, title=title, style=int_style)

        # Save if requested
        if save_html:
            save_html = Path(save_html)
            save_html.parent.mkdir(parents=True, exist_ok=True)
            fig.write_html(str(save_html))

        return fig

class HemisphereVisualizer2D:
    """2D hemisphere visualization using matplotlib.

    Creates publication-quality polar projection plots of hemispherical grids.
    Supports multiple grid types and rendering methods.

    Parameters
    ----------
    grid : HemiGrid
        Hemisphere grid to visualize

    Examples
    --------
    >>> from canvod.grids import create_hemigrid
    >>> from canvod.viz import HemisphereVisualizer2D
    >>>
    >>> grid = create_hemigrid(grid_type='equal_area', angular_resolution=10.0)
    >>> viz = HemisphereVisualizer2D(grid)
    >>> fig, ax = viz.plot_grid_patches(data=vod_data, title="VOD Distribution")
    >>> plt.savefig("vod_plot.png", dpi=300, bbox_inches='tight')

    """

    def __init__(self, grid: HemiGrid) -> None:
        """Initialize 2D hemisphere visualizer.

        Parameters
        ----------
        grid : HemiGrid
            Hemisphere grid to visualize

        """
        self.grid = grid
        self._patches_cache: list[Polygon] | None = None
        self._cell_indices_cache: np.ndarray | None = None

    def plot_grid_patches(
        self,
        data: np.ndarray | None = None,
        style: PolarPlotStyle | None = None,
        ax: Axes | None = None,
        save_path: Path | str | None = None,
        **style_kwargs: Any,
    ) -> tuple[Figure, Axes]:
        """Plot hemisphere grid as colored patches in polar projection.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell. If None, plots uniform grid.
        style : PolarPlotStyle, optional
            Styling configuration. If None, uses defaults.
        ax : matplotlib.axes.Axes, optional
            Existing polar axes to plot on. If None, creates new figure.
        save_path : Path or str, optional
            If provided, saves figure to this path
        **style_kwargs
            Override individual style parameters

        Returns
        -------
        fig : matplotlib.figure.Figure
            Figure object
        ax : matplotlib.axes.Axes
            Polar axes with plot

        Examples
        --------
        >>> fig, ax = viz.plot_grid_patches(
        ...     data=vod_data,
        ...     title="VOD Distribution",
        ...     cmap='plasma',
        ...     save_path="output.png"
        ... )

        """
        # Initialize style
        if style is None:
            style = PolarPlotStyle(**style_kwargs)
        else:
            # Override style with kwargs
            for key, value in style_kwargs.items():
                if hasattr(style, key):
                    setattr(style, key, value)

        # Create figure if needed
        if ax is None:
            fig, ax = plt.subplots(
                figsize=style.figsize, dpi=style.dpi, subplot_kw={"projection": "polar"}
            )
        else:
            fig = cast("Figure", ax.figure)
        ax_polar = cast("PolarAxes", ax)

        # Get patches for grid
        patches, cell_indices = self._extract_grid_patches()

        # Map data to patches
        patch_data = self._map_data_to_patches(data, cell_indices)

        # Determine color limits
        vmin = style.vmin if style.vmin is not None else np.nanmin(patch_data)
        vmax = style.vmax if style.vmax is not None else np.nanmax(patch_data)

        # Create patch collection
        pc = PatchCollection(
            patches,
            cmap=style.cmap,
            edgecolor=style.edgecolor,
            linewidth=style.linewidth,
            alpha=style.alpha,
        )
        pc.set_array(np.ma.masked_invalid(patch_data))
        pc.set_clim(vmin, vmax)

        # Add to axes
        ax.add_collection(pc)

        # Style polar axes
        self._apply_polar_styling(ax_polar, style)

        # Add colorbar
        cbar = fig.colorbar(
            pc,
            ax=ax,
            shrink=style.colorbar_shrink,
            pad=style.colorbar_pad,
        )
        cbar.set_label(style.colorbar_label, fontsize=style.colorbar_fontsize)

        # Set title
        if style.title:
            ax.set_title(style.title, y=1.08, fontsize=14)

        # Save if requested
        if save_path:
            save_path = Path(save_path)
            save_path.parent.mkdir(parents=True, exist_ok=True)
            fig.savefig(
                save_path,
                dpi=style.dpi,
                bbox_inches="tight",
                facecolor="white",
                edgecolor="none",
            )

        return fig, ax

    def _extract_grid_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract 2D polygon patches from hemispherical grid.

        Returns
        -------
        patches : list of Polygon
            Matplotlib polygon patches
        cell_indices : np.ndarray
            Corresponding cell indices in grid

        """
        # Use cache if available
        if self._patches_cache is not None and self._cell_indices_cache is not None:
            return self._patches_cache, self._cell_indices_cache

        grid_type = self.grid.grid_type.lower()

        _rectangular_types = {"equal_area", "equal_angle", "equirectangular"}
        if grid_type in _rectangular_types:
            patches, indices = self._extract_rectangular_patches()
        elif grid_type == "htm":
            patches, indices = self._extract_htm_patches()
        elif grid_type == "geodesic":
            patches, indices = self._extract_geodesic_patches()
        elif grid_type == "healpix":
            patches, indices = self._extract_healpix_patches()
        elif grid_type == "fibonacci":
            patches, indices = self._extract_fibonacci_patches()
        else:
            raise ValueError(f"Unsupported grid type: {grid_type}")

        # Cache results
        self._patches_cache = patches
        self._cell_indices_cache = indices

        return patches, indices

    def _extract_rectangular_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract patches from rectangular/equal-area grid."""
        patches = []
        cell_indices = []

        # Access the grid DataFrame from GridData
        grid_df = self.grid.grid

        for idx, row in enumerate(grid_df.iter_rows(named=True)):
            phi_min = row["phi_min"]
            phi_max = row["phi_max"]
            theta_min = row["theta_min"]
            theta_max = row["theta_max"]

            # Skip cells beyond hemisphere
            if theta_min > np.pi / 2:
                continue

            # Convert to polar coordinates (rho = sin(theta))
            rho_min = np.sin(theta_min)
            rho_max = np.sin(theta_max)

            # Create rectangular patch in polar coordinates
            vertices = np.array(
                [
                    [phi_min, rho_min],
                    [phi_max, rho_min],
                    [phi_max, rho_max],
                    [phi_min, rho_max],
                ]
            )

            patches.append(Polygon(vertices, closed=True))
            cell_indices.append(idx)

        return patches, np.array(cell_indices)

    def _extract_htm_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract triangular patches from HTM grid."""
        patches = []
        cell_indices = []

        # Access the grid DataFrame from GridData
        grid_df = self.grid.grid

        for idx, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                # HTM stores vertices as columns htm_vertex_0, htm_vertex_1,
                # htm_vertex_2.
                v0 = np.array(row["htm_vertex_0"], dtype=float)
                v1 = np.array(row["htm_vertex_1"], dtype=float)
                v2 = np.array(row["htm_vertex_2"], dtype=float)

                vertices_3d = np.array([v0, v1, v2])
                x, y, z = (
                    vertices_3d[:, 0],
                    vertices_3d[:, 1],
                    vertices_3d[:, 2],
                )

                # Convert to spherical coordinates
                r = np.sqrt(x**2 + y**2 + z**2)
                theta = np.arccos(np.clip(z / r, -1, 1))
                phi = np.arctan2(y, x)
                phi = np.mod(phi, 2 * np.pi)

                # Skip if beyond hemisphere
                if np.all(theta > np.pi / 2):
                    continue

                # Convert to polar coordinates (rho = sin(theta))
                rho = np.sin(theta)
                vertices_2d = np.column_stack([phi, rho])

                patches.append(Polygon(vertices_2d, closed=True))
                cell_indices.append(idx)

            except KeyError, TypeError:
                # Skip cells that don't have proper HTM vertex data
                continue

        return patches, np.array(cell_indices)

    def _extract_geodesic_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract triangular patches from geodesic grid.

        The ``geodesic_vertices`` column stores vertex **indices** into the
        shared ``grid.vertices`` coordinate array (shape ``(n_vertices, 3)``).
        """
        patches = []
        cell_indices = []

        grid_df = self.grid.grid
        shared_vertices = self.grid.vertices  # (n_vertices, 3) or None

        if shared_vertices is None or "geodesic_vertices" not in grid_df.columns:
            # No vertex data — fall back to bounding-box rectangles
            return self._extract_rectangular_patches()

        for idx, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                v_indices = np.array(row["geodesic_vertices"], dtype=int)
                if len(v_indices) < 3:
                    continue

                # Look up actual 3D coordinates from shared vertex array
                verts_3d = shared_vertices[v_indices]  # (3, 3)
                x, y, z = verts_3d[:, 0], verts_3d[:, 1], verts_3d[:, 2]

                r = np.sqrt(x**2 + y**2 + z**2)
                theta = np.arccos(np.clip(z / r, -1, 1))
                phi = np.arctan2(y, x)
                phi = np.mod(phi, 2 * np.pi)

                if np.all(theta > np.pi / 2):
                    continue

                rho = np.sin(theta)
                vertices_2d = np.column_stack([phi, rho])
                patches.append(Polygon(vertices_2d, closed=True))
                cell_indices.append(idx)

            except IndexError, KeyError, TypeError, ValueError:
                continue

        return patches, np.array(cell_indices)

    def _extract_healpix_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract patches from HEALPix grid via ``healpy.boundaries``."""
        try:
            import healpy as hp
        except ImportError as e:
            raise ImportError(
                "healpy is required for HEALPix 2D visualization. "
                "Install with: pip install healpy"
            ) from e

        patches = []
        cell_indices = []
        grid_df = self.grid.grid
        nside = int(grid_df["healpix_nside"][0])

        for idx, row in enumerate(grid_df.iter_rows(named=True)):
            ipix = int(row["healpix_ipix"])
            # boundaries returns shape (3, n_vertices) in Cartesian
            boundary = hp.boundaries(nside, ipix, step=4)
            x, y, z = boundary[0], boundary[1], boundary[2]

            # Skip pixels entirely below horizon
            if np.all(z < -0.01):
                continue

            r = np.sqrt(x**2 + y**2 + z**2)
            theta = np.arccos(np.clip(z / r, -1, 1))
            phi = np.arctan2(y, x)
            phi = np.mod(phi, 2 * np.pi)

            # Keep only vertices in upper hemisphere
            mask = theta <= np.pi / 2 + 0.01
            if not np.any(mask):
                continue

            rho = np.sin(theta)
            vertices_2d = np.column_stack([phi, rho])
            patches.append(Polygon(vertices_2d, closed=True))
            cell_indices.append(idx)

        return patches, np.array(cell_indices)

    def _extract_fibonacci_patches(self) -> tuple[list[Polygon], np.ndarray]:
        """Extract patches from Fibonacci grid using Voronoi regions."""
        patches = []
        cell_indices = []
        grid_df = self.grid.grid
        voronoi = self.grid.voronoi  # scipy.spatial.SphericalVoronoi or None

        if voronoi is None or "voronoi_region" not in grid_df.columns:
            # No Voronoi data — fall back to bounding-box rectangles
            return self._extract_rectangular_patches()

        for idx, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                region_indices = row["voronoi_region"]
                if region_indices is None or len(region_indices) < 3:
                    continue

                verts_3d = voronoi.vertices[region_indices]
                x, y, z = verts_3d[:, 0], verts_3d[:, 1], verts_3d[:, 2]

                # Skip cells entirely below horizon
                if np.all(z < -0.01):
                    continue

                r = np.sqrt(x**2 + y**2 + z**2)
                theta = np.arccos(np.clip(z / r, -1, 1))
                phi = np.arctan2(y, x)
                phi = np.mod(phi, 2 * np.pi)

                # Vertices are already in polygon winding order from
                # sort_vertices_of_regions() — use directly.
                rho = np.sin(theta)
                vertices_2d = np.column_stack([phi, rho])
                patches.append(Polygon(vertices_2d, closed=True))
                cell_indices.append(idx)

            except IndexError, KeyError, TypeError, ValueError:
                continue

        return patches, np.array(cell_indices)

    def _map_data_to_patches(
        self,
        data: np.ndarray | None,
        cell_indices: np.ndarray,
    ) -> np.ndarray:
        """Map data values to patches.

        Parameters
        ----------
        data : np.ndarray or None
            Data per grid cell
        cell_indices : np.ndarray
            Cell indices corresponding to patches

        Returns
        -------
        np.ndarray
            Data values for each patch

        """
        if data is None:
            return np.ones(len(cell_indices)) * 0.5

        return data[cell_indices]

    def _apply_polar_styling(
        self,
        ax: PolarAxes,
        style: PolarPlotStyle,
    ) -> None:
        """Apply styling to polar axes.

        Parameters
        ----------
        ax : matplotlib.axes.Axes
            Polar axes to style
        style : PolarPlotStyle
            Styling configuration

        """
        # Set rho limits (0 to 1 for hemisphere projection)
        ax.set_ylim(0, 1.0)

        # Configure polar axis orientation
        ax.set_theta_zero_location("N")  # North at top
        ax.set_theta_direction(-1)  # Clockwise (azimuth convention)

        # Add degree labels on radial axis
        if style.show_degree_labels:
            theta_labels = style.theta_labels
            rho_ticks = [np.sin(np.radians(t)) for t in theta_labels]
            ax.set_yticks(rho_ticks)
            ax.set_yticklabels([f"{t}°" for t in theta_labels])

        # Grid styling
        if style.show_grid:
            ax.grid(
                True,
                alpha=style.grid_alpha,
                linestyle=style.grid_linestyle,
                color="gray",
            )
        else:
            ax.grid(False)

class HemisphereVisualizer3D:
    """3D hemisphere visualization using plotly.

    Creates interactive 3D plots with rotation, zoom, and hover capabilities.
    Designed for exploratory data analysis and presentations.

    Parameters
    ----------
    grid : HemiGrid
        Hemisphere grid to visualize

    Examples
    --------
    >>> from canvod.grids import create_hemigrid
    >>> from canvod.viz import HemisphereVisualizer3D
    >>>
    >>> grid = create_hemigrid(grid_type='equal_area', angular_resolution=10.0)
    >>> viz = HemisphereVisualizer3D(grid)
    >>> fig = viz.plot_hemisphere_surface(data=vod_data, title="Interactive VOD")
    >>> fig.show()

    """

    def __init__(self, grid: HemiGrid) -> None:
        """Initialize 3D hemisphere visualizer.

        Parameters
        ----------
        grid : HemiGrid
            Hemisphere grid to visualize

        """
        self.grid = grid

    def plot_hemisphere_surface(
        self,
        data: np.ndarray | None = None,
        style: PlotStyle | None = None,
        title: str | None = None,
        colorscale: str = "Viridis",
        opacity: float = 0.8,
        show_wireframe: bool = True,
        show_colorbar: bool = True,
        width: int = 800,
        height: int = 600,
        **kwargs: Any,
    ) -> go.Figure:
        """Create 3D surface plot on hemisphere with actual cell patches.

        Renders grid cells as colored 3D patches (not just points).

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell. If None, shows grid structure.
        style : PlotStyle, optional
            Styling configuration. If None, uses defaults.
        title : str, optional
            Plot title
        colorscale : str, default 'Viridis'
            Plotly colorscale name
        opacity : float, default 0.8
            Surface opacity (0=transparent, 1=opaque)
        show_wireframe : bool, default True
            Show grid lines on surface
        show_colorbar : bool, default True
            Display colorbar
        width : int, default 800
            Figure width in pixels
        height : int, default 600
            Figure height in pixels
        **kwargs
            Additional plotly trace parameters

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive 3D figure with cell patches

        Examples
        --------
        >>> fig = viz.plot_hemisphere_surface(
        ...     data=vod_data,
        ...     title="VOD Distribution 3D",
        ...     colorscale='Plasma',
        ...     opacity=0.9
        ... )
        >>> fig.write_html("vod_3d.html")

        """
        # Initialize style
        if style is None:
            style = PlotStyle()
            colorscale = colorscale  # Use parameter
        else:
            colorscale = style.colorscale

        # Render grid based on type
        grid_type = self.grid.grid_type.lower()

        if grid_type in ["equal_area", "equal_angle", "equirectangular"]:
            trace = self._render_rectangular_cells(
                data,
                colorscale,
                opacity,
                show_colorbar,
            )
        elif grid_type == "htm":
            trace = self._render_htm_cells(data, colorscale, opacity, show_colorbar)
        elif grid_type == "geodesic":
            trace = self._render_geodesic_cells(
                data, colorscale, opacity, show_colorbar
            )
        elif grid_type == "healpix":
            trace = self._render_healpix_cells(data, colorscale, opacity, show_colorbar)
        elif grid_type == "fibonacci":
            trace = self._render_fibonacci_cells(
                data, colorscale, opacity, show_colorbar
            )
        else:
            # Fallback to scatter for unknown types
            trace = self._render_scatter_fallback(
                data,
                colorscale,
                opacity,
                show_colorbar,
            )

        fig = go.Figure(data=[trace])

        # Apply layout
        layout_config = style.to_plotly_layout() if style else {}
        layout_config.update(
            {
                "title": title or "Hemisphere 3D",
                "scene": dict(
                    aspectmode="data",
                    xaxis=dict(title="East", showbackground=False),
                    yaxis=dict(title="North", showbackground=False),
                    zaxis=dict(title="Up", showbackground=False),
                    bgcolor=layout_config.get("plot_bgcolor", "white"),
                ),
                "width": width,
                "height": height,
                "margin": dict(l=0, r=0, b=0, t=40),
            }
        )

        fig.update_layout(**layout_config)

        return fig

    def _render_rectangular_cells(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Mesh3d:
        """Render rectangular grid cells as 3D mesh patches."""
        grid_df = self.grid.grid

        all_x, all_y, all_z = [], [], []
        all_i, all_j, all_k = [], [], []
        all_colors = []

        vertex_count = 0

        for i, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                phi_min, phi_max = row["phi_min"], row["phi_max"]
                theta_min, theta_max = row["theta_min"], row["theta_max"]

                # Skip if beyond hemisphere
                if theta_min > np.pi / 2:
                    continue

                # Create 4 corners of rectangular cell
                phi_corners = [phi_min, phi_max, phi_max, phi_min]
                theta_corners = [theta_min, theta_min, theta_max, theta_max]

                patch_x, patch_y, patch_z = [], [], []
                for phi, theta in zip(phi_corners, theta_corners):
                    # Convert to 3D Cartesian
                    x = np.sin(theta) * np.sin(phi)
                    y = np.sin(theta) * np.cos(phi)
                    z = np.cos(theta)
                    patch_x.append(x)
                    patch_y.append(y)
                    patch_z.append(z)

                all_x.extend(patch_x)
                all_y.extend(patch_y)
                all_z.extend(patch_z)

                # Color value for this cell
                color_val = data[i] if data is not None else 0.5
                all_colors.extend([color_val] * 4)

                # Two triangles per rectangle
                all_i.extend([vertex_count, vertex_count])
                all_j.extend([vertex_count + 1, vertex_count + 2])
                all_k.extend([vertex_count + 2, vertex_count + 3])

                vertex_count += 4

            except KeyError, IndexError:
                continue

        return go.Mesh3d(
            x=all_x,
            y=all_y,
            z=all_z,
            i=all_i,
            j=all_j,
            k=all_k,
            intensity=all_colors,
            colorscale=colorscale,
            showscale=show_colorbar,
            colorbar=dict(title="Value") if show_colorbar else None,
            opacity=opacity,
            flatshading=True,
            name=f"{self.grid.grid_type.title()} Grid",
        )

    def _render_htm_cells(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Mesh3d:
        """Render HTM triangular cells as 3D mesh."""
        grid_df = self.grid.grid

        all_x, all_y, all_z = [], [], []
        all_i, all_j, all_k = [], [], []
        all_colors = []

        vertex_count = 0

        for i, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                v0 = np.array(row["htm_vertex_0"], dtype=float)
                v1 = np.array(row["htm_vertex_1"], dtype=float)
                v2 = np.array(row["htm_vertex_2"], dtype=float)

                # Skip if beyond hemisphere
                if np.all([v[2] < 0 for v in [v0, v1, v2]]):
                    continue

                all_x.extend([v0[1], v1[1], v2[1]])
                all_y.extend([v0[0], v1[0], v2[0]])
                all_z.extend([v0[2], v1[2], v2[2]])

                color_val = data[i] if data is not None else 0.5
                all_colors.extend([color_val] * 3)

                # One triangle per cell
                all_i.append(vertex_count)
                all_j.append(vertex_count + 1)
                all_k.append(vertex_count + 2)

                vertex_count += 3

            except KeyError, TypeError, ValueError:
                continue

        return go.Mesh3d(
            x=all_x,
            y=all_y,
            z=all_z,
            i=all_i,
            j=all_j,
            k=all_k,
            intensity=all_colors,
            colorscale=colorscale,
            showscale=show_colorbar,
            colorbar=dict(title="Value") if show_colorbar else None,
            opacity=opacity,
            flatshading=True,
            name=f"{self.grid.grid_type.title()} Grid",
        )

    def _render_geodesic_cells(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Mesh3d | go.Scatter3d:
        """Render geodesic triangular cells as 3D mesh.

        Reads ``geodesic_vertices`` (3 vertex indices per cell) and looks up
        3D Cartesian coordinates from the shared ``grid.vertices`` array.
        """
        grid_df = self.grid.grid
        shared_vertices = self.grid.vertices

        if shared_vertices is None or "geodesic_vertices" not in grid_df.columns:
            return self._render_scatter_fallback(
                data, colorscale, opacity, show_colorbar
            )

        all_x, all_y, all_z = [], [], []
        all_i, all_j, all_k = [], [], []
        all_colors = []
        vertex_count = 0

        for i, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                v_indices = np.array(row["geodesic_vertices"], dtype=int)
                if len(v_indices) < 3:
                    continue

                verts = shared_vertices[v_indices]  # (3, 3)

                if np.all(verts[:, 2] < 0):
                    continue

                all_x.extend(verts[:, 1].tolist())
                all_y.extend(verts[:, 0].tolist())
                all_z.extend(verts[:, 2].tolist())

                color_val = data[i] if data is not None else 0.5
                all_colors.extend([color_val] * 3)

                all_i.append(vertex_count)
                all_j.append(vertex_count + 1)
                all_k.append(vertex_count + 2)
                vertex_count += 3

            except KeyError, IndexError, TypeError, ValueError:
                continue

        return go.Mesh3d(
            x=all_x,
            y=all_y,
            z=all_z,
            i=all_i,
            j=all_j,
            k=all_k,
            intensity=all_colors,
            colorscale=colorscale,
            showscale=show_colorbar,
            colorbar=dict(title="Value") if show_colorbar else None,
            opacity=opacity,
            flatshading=True,
            name="Geodesic Grid",
        )

    def _render_healpix_cells(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Mesh3d | go.Scatter3d:
        """Render HEALPix curvilinear cells as 3D mesh.

        Uses ``healpy.boundaries()`` to obtain true pixel boundaries,
        then fan-triangulates each quadrilateral pixel.
        """
        try:
            import healpy as hp
        except ImportError:
            return self._render_scatter_fallback(
                data, colorscale, opacity, show_colorbar
            )

        grid_df = self.grid.grid
        if "healpix_ipix" not in grid_df.columns:
            return self._render_scatter_fallback(
                data, colorscale, opacity, show_colorbar
            )

        nside = int(grid_df["healpix_nside"][0])
        step = 4  # 4 sub-points per edge → 16 boundary vertices

        all_x, all_y, all_z = [], [], []
        all_i, all_j, all_k = [], [], []
        all_colors = []
        vertex_count = 0

        for i, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                ipix = int(row["healpix_ipix"])
                boundary = hp.boundaries(nside, ipix, step=step)
                x, y, z = boundary[1], boundary[0], boundary[2]

                if np.all(z < -0.01):
                    continue

                n_verts = len(x)
                all_x.extend(x.tolist())
                all_y.extend(y.tolist())
                all_z.extend(z.tolist())

                color_val = data[i] if data is not None else 0.5
                all_colors.extend([color_val] * n_verts)

                # Fan triangulation from first vertex
                for j in range(1, n_verts - 1):
                    all_i.append(vertex_count)
                    all_j.append(vertex_count + j)
                    all_k.append(vertex_count + j + 1)

                vertex_count += n_verts

            except KeyError, IndexError, TypeError, ValueError:
                continue

        return go.Mesh3d(
            x=all_x,
            y=all_y,
            z=all_z,
            i=all_i,
            j=all_j,
            k=all_k,
            intensity=all_colors,
            colorscale=colorscale,
            showscale=show_colorbar,
            colorbar=dict(title="Value") if show_colorbar else None,
            opacity=opacity,
            flatshading=True,
            name="HEALPix Grid",
        )

    def _render_fibonacci_cells(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Mesh3d | go.Scatter3d:
        """Render Fibonacci Voronoi cells as 3D mesh.

        Reads ``voronoi_region`` (variable-length vertex index list) and
        looks up 3D coordinates from ``grid.voronoi.vertices``.  The
        vertex indices are already in correct polygon winding order
        (``SphericalVoronoi.sort_vertices_of_regions()`` was called
        during grid construction), so no re-sorting is needed.
        Fan-triangulates each polygon for ``go.Mesh3d``.
        """
        grid_df = self.grid.grid
        voronoi = self.grid.voronoi

        if voronoi is None or "voronoi_region" not in grid_df.columns:
            return self._render_scatter_fallback(
                data, colorscale, opacity, show_colorbar
            )

        all_x, all_y, all_z = [], [], []
        all_i, all_j, all_k = [], [], []
        all_colors = []
        vertex_count = 0

        for i, row in enumerate(grid_df.iter_rows(named=True)):
            try:
                region_indices = row["voronoi_region"]
                if region_indices is None or len(region_indices) < 3:
                    continue

                verts = voronoi.vertices[region_indices]

                if np.all(verts[:, 2] < -0.01):
                    continue

                # Vertices are already in polygon winding order from
                # sort_vertices_of_regions() — use directly.
                n_verts = len(verts)

                all_x.extend(verts[:, 1].tolist())
                all_y.extend(verts[:, 0].tolist())
                all_z.extend(verts[:, 2].tolist())

                color_val = data[i] if data is not None else 0.5
                all_colors.extend([color_val] * n_verts)

                # Fan triangulation from first vertex
                for j in range(1, n_verts - 1):
                    all_i.append(vertex_count)
                    all_j.append(vertex_count + j)
                    all_k.append(vertex_count + j + 1)

                vertex_count += n_verts

            except KeyError, IndexError, TypeError, ValueError:
                continue

        return go.Mesh3d(
            x=all_x,
            y=all_y,
            z=all_z,
            i=all_i,
            j=all_j,
            k=all_k,
            intensity=all_colors,
            colorscale=colorscale,
            showscale=show_colorbar,
            colorbar=dict(title="Value") if show_colorbar else None,
            opacity=opacity,
            flatshading=True,
            name="Fibonacci Grid",
        )

    def _render_scatter_fallback(
        self,
        data: np.ndarray | None,
        colorscale: str,
        opacity: float,
        show_colorbar: bool,
    ) -> go.Scatter3d:
        """Fallback scatter plot for unsupported grid types."""
        grid_df = self.grid.grid
        theta = grid_df["theta"].to_numpy()
        phi = grid_df["phi"].to_numpy()

        # Convert to 3D Cartesian coordinates
        x = np.sin(theta) * np.sin(phi)
        y = np.sin(theta) * np.cos(phi)
        z = np.cos(theta)

        # Prepare data values
        if data is None:
            values = np.ones(self.grid.ncells) * 0.5
        else:
            values = data

        # Filter hemisphere only
        hemisphere_mask = theta <= np.pi / 2
        x = x[hemisphere_mask]
        y = y[hemisphere_mask]
        z = z[hemisphere_mask]
        values = values[hemisphere_mask]

        return go.Scatter3d(
            x=x,
            y=y,
            z=z,
            mode="markers",
            marker=dict(
                size=6,
                color=values,
                colorscale=colorscale,
                opacity=opacity,
                colorbar=dict(title="Value") if show_colorbar else None,
                cmin=np.nanmin(values),
                cmax=np.nanmax(values),
            ),
            text=[f"Cell {i}<br>Value: {v:.3f}" for i, v in enumerate(values)],
            hoverinfo="text",
        )

    def plot_hemisphere_scatter(
        self,
        data: np.ndarray | None = None,
        title: str | None = None,
        colorscale: str = "Viridis",
        marker_size: int | np.ndarray = 6,
        opacity: float = 0.8,
        width: int = 800,
        height: int = 600,
    ) -> go.Figure:
        """Create 3D scatter plot of cell centers.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, optional
            Plot title
        colorscale : str, default 'Viridis'
            Plotly colorscale name
        marker_size : int or np.ndarray, default 6
            Marker size (constant or per-point array)
        opacity : float, default 0.8
            Marker opacity
        width : int, default 800
            Figure width
        height : int, default 600
            Figure height

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive scatter plot

        """
        # Note: Now renders as mesh, not scatter points
        # Marker size parameter is ignored
        fig = self.plot_hemisphere_surface(
            data=data,
            title=title,
            colorscale=colorscale,
            opacity=opacity,
            width=width,
            height=height,
        )

        return fig

    def plot_cell_mesh(
        self,
        data: np.ndarray | None = None,
        title: str | None = None,
        colorscale: str = "Viridis",
        opacity: float = 0.7,
        show_edges: bool = True,
        width: int = 800,
        height: int = 600,
    ) -> go.Figure:
        """Create 3D mesh plot showing cell boundaries.

        Parameters
        ----------
        data : np.ndarray, optional
            Data values per cell
        title : str, optional
            Plot title
        colorscale : str, default 'Viridis'
            Plotly colorscale name
        opacity : float, default 0.7
            Mesh opacity
        show_edges : bool, default True
            Show cell edges
        width : int, default 800
            Figure width
        height : int, default 600
            Figure height

        Returns
        -------
        plotly.graph_objects.Figure
            Interactive mesh plot

        Notes
        -----
        This method requires grid cells with vertex information.
        Currently supports HTM and geodesic grids.

        """
        traces = []

        # Prepare data
        if data is None:
            values = np.ones(self.grid.ncells) * 0.5
        else:
            values = data

        grid_df = self.grid.grid
        grid_type = self.grid.grid_type.lower()

        # Check if grid supports mesh rendering
        if grid_type == "htm" and "htm_vertex_0" in grid_df.columns:
            # HTM triangular mesh
            for idx, row in enumerate(grid_df.iter_rows(named=True)):
                try:
                    v0 = np.array(row["htm_vertex_0"], dtype=float)
                    v1 = np.array(row["htm_vertex_1"], dtype=float)
                    v2 = np.array(row["htm_vertex_2"], dtype=float)

                    vertices = np.array([v0, v1, v2])

                    # Check hemisphere
                    z_coords = vertices[:, 2]
                    if np.all(z_coords < 0):
                        continue

                    # Normalize color value
                    color_val = (values[idx] - np.nanmin(values)) / (
                        np.nanmax(values) - np.nanmin(values)
                    )
                    color_rgb = sample_colorscale(colorscale, [color_val])[0]

                    # Create triangle mesh
                    trace = go.Mesh3d(
                        x=vertices[:, 0],
                        y=vertices[:, 1],
                        z=vertices[:, 2],
                        i=[0],
                        j=[1],
                        k=[2],
                        color=color_rgb,
                        opacity=opacity,
                        flatshading=True,
                        showscale=False,
                        hoverinfo="skip",
                    )
                    traces.append(trace)
                except KeyError, TypeError, ValueError:
                    continue
        else:
            # Not supported for this grid type
            raise NotImplementedError(
                f"Cell mesh rendering not implemented for {grid_type} grids"
            )

        # Sample colorscale (no longer needed above, but kept for compatibility)

        # Add colorbar trace
        if values is not None and len(traces) > 0:
            dummy_trace = go.Scatter3d(
                x=[None],
                y=[None],
                z=[None],
                mode="markers",
                marker=dict(
                    size=0.1,
                    color=[np.nanmin(values), np.nanmax(values)],
                    colorscale=colorscale,
                    colorbar=dict(title="Value"),
                ),
                showlegend=False,
                hoverinfo="skip",
            )
            traces.append(dummy_trace)

        fig = go.Figure(data=traces)

        # Update layout
        fig.update_layout(
            title=title or "Hemisphere Mesh 3D",
            scene=dict(
                aspectmode="data",
                xaxis=dict(title="East", showbackground=False),
                yaxis=dict(title="North", showbackground=False),
                zaxis=dict(title="Up", showbackground=False),
            ),
            width=width,
            height=height,
            margin=dict(l=0, r=0, b=0, t=40),
        )

        return fig

    def add_spherical_overlays(
        self,
        fig: go.Figure,
        elevation_rings: list[int] | None = None,
        meridians_deg: list[int] | None = None,
        overlay_color: str = "lightgray",
        line_width: float = 1,
    ) -> go.Figure:
        """Add elevation rings and meridians to 3D plot.

        Parameters
        ----------
        fig : plotly.graph_objects.Figure
            Existing 3D figure to add overlays to
        elevation_rings : list of int, optional
            Elevation angles in degrees. Default: [15, 30, 45, 60, 75, 90]
        meridians_deg : list of int, optional
            Meridian angles in degrees. Default: [0, 45, 90, 135, 180, 225, 270, 315]
        overlay_color : str, default 'lightgray'
            Color for overlay lines
        line_width : float, default 1
            Width of overlay lines

        Returns
        -------
        fig : plotly.graph_objects.Figure
            Modified figure with overlays

        """
        if elevation_rings is None:
            elevation_rings = [15, 30, 45, 60, 75, 90]
        if meridians_deg is None:
            meridians_deg = list(range(0, 360, 45))

        # Elevation rings
        for theta_deg in elevation_rings:
            theta = np.radians(theta_deg)
            phi = np.linspace(0, 2 * np.pi, 200)
            x = np.sin(theta) * np.sin(phi)
            y = np.sin(theta) * np.cos(phi)
            z = np.full_like(phi, np.cos(theta))
            fig.add_trace(
                go.Scatter3d(
                    x=x,
                    y=y,
                    z=z,
                    mode="lines",
                    line=dict(color=overlay_color, width=line_width),
                    hoverinfo="skip",
                    showlegend=False,
                )
            )

        # Meridians
        for phi_deg in meridians_deg:
            phi = np.radians(phi_deg)
            theta = np.linspace(0, np.pi / 2, 100)
            x = np.sin(theta) * np.sin(phi)
            y = np.sin(theta) * np.cos(phi)
            z = np.cos(theta)
            fig.add_trace(
                go.Scatter3d(
                    x=x,
                    y=y,
                    z=z,
                    mode="lines",
                    line=dict(color=overlay_color, width=line_width),
                    hoverinfo="skip",
                    showlegend=False,
                )
            )

        return fig

    def add_custom_axes(
        self,
        fig: go.Figure,
        axis_length: float = 1.2,
        axis_color: str = "black",
        show_labels: bool = True,
    ) -> go.Figure:
        """Add custom coordinate axes with labels.

        Parameters
        ----------
        fig : plotly.graph_objects.Figure
            Existing 3D figure
        axis_length : float, default 1.2
            Length of axis lines
        axis_color : str, default 'black'
            Color for axes
        show_labels : bool, default True
            Show axis labels (E, N, Z)

        Returns
        -------
        fig : plotly.graph_objects.Figure
            Modified figure with custom axes

        """
        # Axis lines
        axes_lines = [
            dict(x=[0, axis_length], y=[0, 0], z=[0, 0]),  # East
            dict(x=[0, 0], y=[0, axis_length], z=[0, 0]),  # North
            dict(x=[0, 0], y=[0, 0], z=[0, axis_length]),  # Up
        ]

        for axis in axes_lines:
            fig.add_trace(
                go.Scatter3d(
                    x=axis["x"],
                    y=axis["y"],
                    z=axis["z"],
                    mode="lines",
                    line=dict(color=axis_color, width=6),
                    hoverinfo="skip",
                    showlegend=False,
                )
            )

        # Arrowheads
        arrow_tip = axis_length + 0.05
        arrow_size = 0.1
        for pos, direction in zip(
            [[arrow_tip, 0, 0], [0, arrow_tip, 0], [0, 0, arrow_tip]],
            [[arrow_size, 0, 0], [0, arrow_size, 0], [0, 0, arrow_size]],
        ):
            fig.add_trace(
                go.Cone(
                    x=[pos[0]],
                    y=[pos[1]],
                    z=[pos[2]],
                    u=[direction[0]],
                    v=[direction[1]],
                    w=[direction[2]],
                    sizemode="absolute",
                    sizeref=arrow_size,
                    anchor="tip",
                    showscale=False,
                    colorscale=[[0, axis_color], [1, axis_color]],
                    hoverinfo="skip",
                    showlegend=False,
                )
            )

        # Labels
        if show_labels:
            label_offset = axis_length + 0.15
            for label in [
                dict(x=label_offset, y=0, z=0, text="E"),
                dict(x=0, y=label_offset, z=0, text="N"),
                dict(x=0, y=0, z=label_offset, text="Z"),
            ]:
                fig.add_trace(
                    go.Scatter3d(
                        x=[label["x"]],
                        y=[label["y"]],
                        z=[label["z"]],
                        mode="text",
                        text=[label["text"]],
                        textfont=dict(size=16, color=axis_color),
                        hoverinfo="skip",
                        showlegend=False,
                    )
                )

        return fig