Move projection-based axis labels and aspect ratio onto geometry

src/mikeio/dataset/_data_plot.py

@@ -317,7 +317,7 @@ def contour(
         pos = ax.contour(x, y, da.values, **kwargs)
         # fig.colorbar(pos, label=self._label_txt())
         ax.clabel(pos, fmt="%1.2f", inline=1, fontsize=9)
-        self._set_aspect_and_labels(ax, self.da.geometry, y)
+        self._set_aspect_and_labels(ax, self.da.geometry)
         if title is not None:
             ax.set_title(title)
         return ax
@@ -349,7 +349,7 @@ def contourf(
 
         pos = ax.contourf(x, y, da.values, **kwargs)
         fig.colorbar(pos, label=label, pad=0.01)
-        self._set_aspect_and_labels(ax, self.da.geometry, y)
+        self._set_aspect_and_labels(ax, self.da.geometry)
         if title is not None:
             ax.set_title(title)
         return ax
@@ -381,7 +381,7 @@ def pcolormesh(
 
         pos = ax.pcolormesh(xn, yn, da.values, **kwargs)
         fig.colorbar(pos, label=label, pad=0.01)
-        self._set_aspect_and_labels(ax, self.da.geometry, yn)
+        self._set_aspect_and_labels(ax, self.da.geometry)
         if title is not None:
             ax.set_title(title)
         return ax
@@ -397,26 +397,13 @@ def _get_xn_yn(self) -> tuple[np.ndarray, np.ndarray]:
         return xn, yn
 
     @staticmethod
-    def _set_aspect_and_labels(ax: Axes, geometry: Any, y: np.ndarray) -> None:
-        if geometry.is_spectral:
-            ax.set_xlabel("Frequency [Hz]")
-            ax.set_ylabel("Directions [degree]")
-        elif geometry._is_rotated:
-            ax.set_xlabel("[m]")
-            ax.set_ylabel("[m]")
-        elif geometry.projection == "NON-UTM":
-            ax.set_xlabel("[m]")
-            ax.set_ylabel("[m]")
-        elif geometry.is_geo:
-            ax.set_xlabel("Longitude [degrees]")
-            ax.set_ylabel("Latitude [degrees]")
-            mean_lat = np.mean(y)
-            aspect_ratio = 1.0 / np.cos(np.pi * mean_lat / 180)
-            ax.set_aspect(aspect_ratio)
-        else:
-            ax.set_xlabel("Easting [m]")
-            ax.set_ylabel("Northing [m]")
-            ax.set_aspect("equal")
+    def _set_aspect_and_labels(ax: Axes, geometry: Any) -> None:
+        xlabel, ylabel = geometry._axis_labels
+        ax.set_xlabel(xlabel)
+        ax.set_ylabel(ylabel)
+        aspect = geometry._plot_aspect
+        if aspect is not None:
+            ax.set_aspect(aspect)
 
 
 class DataArrayPlotterFM(DataArrayPlotter):

src/mikeio/spatial/_FM_geometry.py

@@ -26,10 +26,9 @@
     _plot_map,
     BoundaryPolygons,
     Polygon,
-    _set_xy_label_by_projection,  # TODO remove
     _to_polygons,  # TODO remove
 )
-from ._geometry import Geometry0D, GeometryPoint2D, _Geometry
+from ._geometry import Geometry0D, GeometryPoint2D, _Geometry, _geographic_aspect
 
 from ._grid_geometry import Grid2D
 from ._distance import xy_to_bbox
@@ -152,7 +151,7 @@ def mesh(
         from matplotlib.collections import PatchCollection  # type: ignore
 
         ax = self._get_ax(ax=ax, figsize=figsize)
-        ax.set_aspect(self._plot_aspect())
+        ax.set_aspect(self.g._plot_aspect)
 
         patches = _to_polygons(self.g.node_coordinates, self.g.element_table)
         fig_obj = PatchCollection(
@@ -162,7 +161,9 @@ def mesh(
         self.outline(ax=ax)
         ax.set_title(title)
         ax = self._set_plot_limits(ax)
-        _set_xy_label_by_projection(ax, self.g.projection)
+        xlabel, ylabel = self.g._axis_labels
+        ax.set_xlabel(xlabel)
+        ax.set_ylabel(ylabel)
         return ax
 
     def outline(
@@ -181,7 +182,7 @@ def outline(
 
         """
         ax = self._get_ax(ax=ax, figsize=figsize)
-        ax.set_aspect(self._plot_aspect())
+        ax.set_aspect(self.g._plot_aspect)
 
         linwid = 1.2
         out_col = "0.4"
@@ -212,7 +213,7 @@ def boundary_nodes(
         import matplotlib.pyplot as plt
 
         ax = self._get_ax(ax=ax, figsize=figsize)
-        ax.set_aspect(self._plot_aspect())
+        ax.set_aspect(self.g._plot_aspect)
 
         nc = self.g.node_coordinates
         c = self.g.codes
@@ -247,13 +248,6 @@ def _set_plot_limits(self, ax: Axes) -> Axes:
         ax.set_ylim(bbox.bottom - xybuf, bbox.top + xybuf)
         return ax
 
-    def _plot_aspect(self) -> Literal["equal"] | float:
-        if self.g.is_geo:
-            mean_lat = np.mean(self.g.node_coordinates[:, 1])
-            return 1.0 / np.cos(np.pi * mean_lat / 180)
-        else:
-            return "equal"
-
 
 class _GeometryFM(_Geometry):
     def __init__(
@@ -292,6 +286,12 @@ def __init__(
         if reindex:
             self._reindex()
 
+    @property
+    def _plot_aspect(self) -> Literal["equal"] | float:
+        if self.is_geo:
+            return _geographic_aspect(self.node_coordinates[:, 1])
+        return "equal"
+
     def _calc_element_coordinates(self) -> NDArray[np.floating]:
         element_table = self.element_table
         node_coords = self.node_coordinates

src/mikeio/spatial/_FM_plot.py

@@ -130,8 +130,10 @@ def _plot_map(
     if ax is None:
         _, ax = plt.subplots(figsize=figsize)
 
-    _set_aspect_ratio(ax, nc, geometry.projection)
-    _set_xy_label_by_projection(ax, geometry.projection)
+    ax.set_aspect(geometry._plot_aspect)
+    xlabel, ylabel = geometry._axis_labels
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel(ylabel)
 
     if plot_type == "outline_only":
         _plot_outline_only(ax, boundary_polylines)
@@ -372,15 +374,6 @@ def _add_colorbar(
     )
 
 
-def _set_aspect_ratio(ax: Axes, nc: np.ndarray, projection: str) -> None:
-    is_geo = projection == "LONG/LAT"
-    if is_geo:
-        mean_lat = np.mean(nc[:, 1])
-        ax.set_aspect(1.0 / np.cos(np.pi * mean_lat / 180))
-    else:
-        ax.set_aspect("equal")
-
-
 def _add_non_tri_mesh(
     ax: Axes, nc: np.ndarray, element_table: np.ndarray, plot_type: str
 ) -> None:
@@ -407,18 +400,6 @@ def _add_outline(ax: Axes, boundary_polylines: list[Polygon]) -> None:
         ax.plot(*line.xy.T, color="0.4", linewidth=1.2)
 
 
-def _set_xy_label_by_projection(ax: Axes, projection: str) -> None:
-    if (not projection) or projection == "NON-UTM":
-        ax.set_xlabel("x [m]")
-        ax.set_ylabel("y [m]")
-    elif projection == "LONG/LAT":
-        ax.set_xlabel("Longitude [degrees]")
-        ax.set_ylabel("Latitude [degrees]")
-    else:
-        ax.set_xlabel("Easting [m]")
-        ax.set_ylabel("Northing [m]")
-
-
 def _is_tri_only(element_table: np.ndarray) -> bool:
     return max([len(el) for el in element_table]) == 3
 

src/mikeio/spatial/_geometry.py

@@ -4,9 +4,16 @@
 from dataclasses import dataclass
 from typing import Any, Sequence
 
+import numpy as np
 from mikecore.Projections import MapProjection
 
 
+def _geographic_aspect(latitudes: np.ndarray) -> float:
+    """Aspect ratio correcting for longitude compression at a given latitude."""
+    mean_lat = float(np.mean(latitudes))
+    return 1.0 / np.cos(np.pi * mean_lat / 180)
+
+
 @dataclass
 class BoundingBox:
     """Bounding box for spatial data."""
@@ -79,6 +86,15 @@ def is_local_coordinates(self) -> bool:
         """Are coordinates relative (NON-UTM)?"""
         return self._projstr == "NON-UTM"
 
+    @property
+    def _axis_labels(self) -> tuple[str, str]:
+        """(x, y) axis labels for plotting, derived from the projection."""
+        if self.is_geo:
+            return "Longitude [degrees]", "Latitude [degrees]"
+        if self.is_local_coordinates:  # NON-UTM
+            return "x [m]", "y [m]"
+        return "Easting [m]", "Northing [m]"
+
     @property
     @abstractmethod
     def dims(self) -> tuple[str, ...]:

src/mikeio/spatial/_grid_geometry.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 from functools import cached_property
 from pathlib import Path
-from typing import Any, Sequence, TYPE_CHECKING, overload
+from typing import Any, Literal, Sequence, TYPE_CHECKING, overload
 from dataclasses import dataclass
 import numpy as np
 
@@ -16,6 +16,7 @@
     GeometryPoint3D,
     GeometryUndefined,
     _Geometry,
+    _geographic_aspect,
 )
 
 from .._interpolation import Interpolant
@@ -375,26 +376,12 @@ def outline(
         return ax
 
     def _set_aspect_and_labels(self, ax: Axes) -> None:
-        g = self.g
-        if g.is_spectral:
-            ax.set_xlabel("Frequency [Hz]")
-            ax.set_ylabel("Directions [degree]")
-        elif g._is_rotated:
-            ax.set_xlabel("[m]")
-            ax.set_ylabel("[m]")
-        elif g.projection == "NON-UTM":
-            ax.set_xlabel("[m]")
-            ax.set_ylabel("[m]")
-        elif g.is_geo:
-            ax.set_xlabel("Longitude [degrees]")
-            ax.set_ylabel("Latitude [degrees]")
-            mean_lat = np.mean(g.y)
-            aspect_ratio = 1.0 / np.cos(np.pi * mean_lat / 180)
-            ax.set_aspect(aspect_ratio)
-        else:
-            ax.set_xlabel("Easting [m]")
-            ax.set_ylabel("Northing [m]")
-            ax.set_aspect("equal")
+        xlabel, ylabel = self.g._axis_labels
+        ax.set_xlabel(xlabel)
+        ax.set_ylabel(ylabel)
+        aspect = self.g._plot_aspect
+        if aspect is not None:
+            ax.set_aspect(aspect)
 
 
 @dataclass
@@ -507,6 +494,22 @@ def dims(self) -> tuple[str, ...]:
     def _is_rotated(self) -> Any:
         return np.abs(self._orientation) > 1e-5
 
+    @property
+    def _axis_labels(self) -> tuple[str, str]:
+        if self.is_spectral:
+            return "Frequency [Hz]", "Directions [degree]"
+        if self._is_rotated:
+            return "[m]", "[m]"
+        return super()._axis_labels
+
+    @property
+    def _plot_aspect(self) -> Literal["equal"] | float | None:
+        if self.is_spectral or self._is_rotated or self.is_local_coordinates:
+            return None
+        if self.is_geo:
+            return _geographic_aspect(self.y)
+        return "equal"
+
     def _create_in_bbox(
         self,
         bbox: BoundingBox | tuple[float, float, float, float] | Sequence[float],