fix: GaussianGrid shape of latitude and longitude should match grid shape

examples/plot_map.py

@@ -63,17 +63,15 @@
     num_y, num_x = data.shape
     grid = OmGrid(reader.get_child_by_name("crs_wkt").read_scalar(), (num_y, num_x))
     lon_grid, lat_grid = grid.get_meshgrid()
-    crs = grid.crs
-    if crs is None:
-        raise ValueError("CRS is None, this should only happen for gaussian grids")
+    crs_name = grid.crs.name if not grid.crs is None else "Gaussian"
 
     # Plot the data
     im = ax.contourf(lon_grid, lat_grid, data, cmap="coolwarm")
     ax.gridlines(draw_labels=True, alpha=0.3)
     ax.set_aspect("equal")
     # ax.set_global()
     plt.colorbar(im, ax=ax, orientation="vertical", pad=0.05, aspect=40, shrink=0.55, label=VARIABLE)
-    plt.title(f"{MODEL_DOMAIN} {VARIABLE} Map\nCRS: {crs.name}", fontsize=12, fontweight="bold", pad=16)
+    plt.title(f"{MODEL_DOMAIN} {VARIABLE} Map\nCRS: {crs_name}", fontsize=12, fontweight="bold", pad=16)
     plt.tight_layout()
 
     output_filename = f"map_{VARIABLE}.png"

examples/regrid_subset_of_projected_domain.py

@@ -79,9 +79,6 @@ def find_boundary_grid_indices(grid: OmGrid, lons: np.ndarray, lats: np.ndarray)
     # Create coordinate arrays
     num_y, num_x = child.shape
     grid = OmGrid(reader.get_child_by_name("crs_wkt").read_scalar(), (num_y, num_x))
-    source_crs = grid.crs
-    if source_crs is None:
-        raise ValueError("CRS is None, this should only happen for gaussian grids")
 
     xmin, xmax, ymin, ymax = find_boundary_grid_indices(grid, TARGET_LONS, TARGET_LATS)
 

python/omfiles/grids/gaussian.py

@@ -740,14 +740,14 @@ def _find_point_xy(self, lat: float, lon: float) -> XYIndex:
 
     @property
     def latitude(self) -> npt.NDArray[np.float64]:
-        """Get 1D array of latitude coordinates for all grid points."""
+        """Get 2D array of latitude coordinates for all grid points."""
         if not hasattr(self, "_latitude"):
             self._compute_coordinates()
         return self._latitude
 
     @property
     def longitude(self) -> npt.NDArray[np.float64]:
-        """Get 1D array of longitude coordinates for all grid points."""
+        """Get 2D array of longitude coordinates for all grid points."""
         if not hasattr(self, "_longitude"):
             self._compute_coordinates()
         return self._longitude
@@ -762,13 +762,5 @@ def _compute_coordinates(self) -> None:
             lats[gridpoint] = lat
             lons[gridpoint] = lon
 
-        self._latitude = lats
-        self._longitude = lons
-
-    def get_meshgrid(self) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
-        """
-        Meshgrids are not meaningful for Gaussian grids, because the grid points are not evenly spaced.
-        """
-        raise NotImplementedError(
-            "Meshgrids are not meaningful for Gaussian grids. Use earthkit.regrid to regrid to a regular grid."
-        )
+        self._latitude = lats.reshape(self.shape)
+        self._longitude = lons.reshape(self.shape)

python/omfiles/grids/om_grid.py

@@ -64,8 +64,15 @@ def get_coordinates(self, x: int, y: int) -> LatLon:
         return self._grid.get_coordinates(x, y)
 
     def get_meshgrid(self) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
-        """Get meshgrid of geographic coordinates."""
-        return self._grid.get_meshgrid()
+        """
+        Get meshgrid of geographic coordinates.
+
+        Useful for plotting with matplotlib/cartopy.
+
+        Returns:
+            (lon_grid, lat_grid) arrays of shape (ny, nx) in geographic coordinates
+        """
+        return (self._grid.longitude, self._grid.latitude)
 
     @property
     def is_gaussian(self) -> bool:

python/omfiles/grids/regular.py

@@ -134,14 +134,3 @@ def get_coordinates(self, x: int, y: int) -> LatLon:
         lon, lat = self.to_wgs84.transform(x_proj, y_proj)
 
         return LatLon(float(lat), float(lon))
-
-    def get_meshgrid(self) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
-        """
-        Get meshgrid of geographic coordinates.
-
-        Useful for plotting with matplotlib/cartopy.
-
-        Returns:
-            (lon_grid, lat_grid) arrays of shape (ny, nx) in geographic coordinates
-        """
-        return (self.longitude, self.latitude)

tests/test_grids.py

@@ -74,6 +74,9 @@ def test_regular_grid(icon_global_grid: OmGrid):
     assert icon_global_grid.find_point_xy(90, -180) == (0, 1440)
     assert icon_global_grid.find_point_xy(90, 179.75) == (2878, 1440)
     assert icon_global_grid.find_point_xy(0, 0) == (1440, 720)
+    assert icon_global_grid.shape == (1441, 2879)
+    assert icon_global_grid.latitude.shape == (1441, 2879)
+    assert icon_global_grid.longitude.shape == (1441, 2879)
 
 
 def test_regular_grid_roundtrip(icon_global_grid: OmGrid):
@@ -319,3 +322,7 @@ def test_ecmwf_grid(ecmwf_ifs_grid: GaussianGrid):
     position = ecmwf_ifs_grid.get_coordinates(flat_grid_coords[1], flat_grid_coords[0])
     assert abs(position[0] - 89.94619) < 0.005
     assert abs(position[1] - 0) < 0.005
+
+    assert ecmwf_ifs_grid.shape == (1, 6599680)
+    assert ecmwf_ifs_grid.latitude.shape == (1, 6599680)
+    assert ecmwf_ifs_grid.longitude.shape == (1, 6599680)