diff --git a/docs/user-guide/grid-formats.rst b/docs/user-guide/grid-formats.rst
index 7152b1f5e..b252ce764 100644
--- a/docs/user-guide/grid-formats.rst
+++ b/docs/user-guide/grid-formats.rst
@@ -17,6 +17,7 @@ NetCDF file format. As of the most recent release, the following grid formats ar
* SCRIP
* EXODUS
* ESMF
+* GEOS CS
While each of these formats can be encoded in the UGRID conventions, the amount of information that is parsed from them
varies. The following sections describes how each format is converted into the UGRID conventions and what variables
@@ -111,6 +112,17 @@ References
* https://earthsystemmodeling.org/about/
* https://earthsystemmodeling.org/docs/release/ESMF_8_0_1/ESMF_refdoc/node3.html#SECTION03028200000000000000
+GEOS CS
+=======
+
+The Goddard Earth Observing System (GEOS) Cube Sphere (CS) grid format is equidistant gnomonic cubed-sphere grid with
+6 identical faces that wrap a sphere, with some number of grid cells per face. For example, a C720 GEOS-CS grid has
+6 faces, each with 720x720 elements.
+
+References
+----------
+* https://gmao.gsfc.nasa.gov/gmaoftp/ops/GEOSIT_sample/doc/CS_Description_c180_v1.pdf
+
Parsed Variables
================
@@ -157,6 +169,7 @@ Coordinates
| SCRIP |
EXODUS |
ESMF |
+ GEOS-CS |
| node_latlon |
@@ -165,6 +178,7 @@ Coordinates
Yes |
No |
Yes |
+ Yes |
| edge_latlon |
@@ -173,6 +187,8 @@ Coordinates
No |
No |
No |
+ No |
+
| face_latlon |
@@ -181,6 +197,7 @@ Coordinates
Yes |
No |
Yes |
+ Yes |
| node_xyz |
@@ -189,6 +206,7 @@ Coordinates
No |
Yes |
No |
+ No |
| edge_xyz |
@@ -197,6 +215,7 @@ Coordinates
No |
Yes |
No |
+ No |
| face_xyz |
@@ -205,6 +224,7 @@ Coordinates
No |
Yes |
No |
+ No |
@@ -242,6 +262,7 @@ Connectivity
SCRIP |
EXODUS |
ESMF |
+ GEOS-CS |
| face_node |
@@ -250,6 +271,7 @@ Connectivity
Yes |
Yes |
Yes |
+ Yes |
| face_edge |
@@ -258,6 +280,7 @@ Connectivity
No |
No |
No |
+ No |
| face_face |
@@ -266,6 +289,7 @@ Connectivity
No |
No |
No |
+ No |
| edge_node |
@@ -274,6 +298,7 @@ Connectivity
No |
No |
No |
+ No |
| edge_edge |
@@ -282,6 +307,7 @@ Connectivity
No |
No |
No |
+ No |
| edge_face |
@@ -290,6 +316,7 @@ Connectivity
No |
No |
No |
+ No |
| node_node |
@@ -298,6 +325,7 @@ Connectivity
No |
No |
No |
+ No |
| node_edge |
@@ -306,6 +334,7 @@ Connectivity
No |
No |
No |
+ No |
| node_face |
@@ -314,5 +343,6 @@ Connectivity
No |
No |
No |
+ No |
diff --git a/test/meshfiles/geos-cs/c12/test-c12.native.nc4 b/test/meshfiles/geos-cs/c12/test-c12.native.nc4
new file mode 100644
index 000000000..9263df454
Binary files /dev/null and b/test/meshfiles/geos-cs/c12/test-c12.native.nc4 differ
diff --git a/test/test_geos.py b/test/test_geos.py
new file mode 100644
index 000000000..95ef721d5
--- /dev/null
+++ b/test/test_geos.py
@@ -0,0 +1,31 @@
+import uxarray as ux
+
+import os
+from pathlib import Path
+
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+gridfile_geos_cs = current_path / "meshfiles" / "geos-cs" / "c12" / "test-c12.native.nc4"
+
+
+
+def test_read_geos_cs_grid():
+ """Tests the conversion of a CS12 GEOS-CS Grid to the UGRID conventions.
+
+ A CS12 grid has 6 faces, each with 12x12 faces and 13x13 nodes each.
+ """
+
+ uxgrid = ux.open_grid(gridfile_geos_cs)
+
+ n_face = 6 * 12 * 12
+ n_node = 6 * 13 * 13
+
+ assert uxgrid.n_face == n_face
+ assert uxgrid.n_node == n_node
+
+
+def test_read_geos_cs_uxds():
+ """Tests the creating of a UxDataset from a CS12 GEOS-CS Grid."""
+ uxds = ux.open_dataset(gridfile_geos_cs, gridfile_geos_cs)
+
+ assert uxds['T'].shape[-1] == uxds.uxgrid.n_face
diff --git a/uxarray/core/utils.py b/uxarray/core/utils.py
index 14d569094..550dc1012 100644
--- a/uxarray/core/utils.py
+++ b/uxarray/core/utils.py
@@ -7,28 +7,36 @@ def _map_dims_to_ugrid(
remaps the original dimension name to match the UGRID conventions (i.e.
"nCell": "n_face")"""
- keys_to_drop = []
- for key in _source_dims_dict.keys():
- # obtain all dimensions not present in the original dataset
- if key not in ds.dims:
- keys_to_drop.append(key)
+ if grid.source_grid_spec == "GEOS-CS":
+ # stack dimensions to flatten them to map to nodes or faces
+ for var_name in list(ds.coords) + list(ds.data_vars):
+ if all(key in ds[var_name].sizes for key in ["nf", "Ydim", "Xdim"]):
+ ds[var_name] = ds[var_name].stack(n_face=["nf", "Ydim", "Xdim"])
+ if all(key in ds[var_name].sizes for key in ["nf", "YCdim", "XCdim"]):
+ ds[var_name] = ds[var_name].stack(n_node=["nf", "YCdim", "XCdim"])
+ else:
+ keys_to_drop = []
+ for key in _source_dims_dict.keys():
+ # obtain all dimensions not present in the original dataset
+ if key not in ds.dims:
+ keys_to_drop.append(key)
- for key in keys_to_drop:
- # drop dimensions not present in the original dataset
- _source_dims_dict.pop(key)
+ for key in keys_to_drop:
+ # drop dimensions not present in the original dataset
+ _source_dims_dict.pop(key)
- for dim in set(ds.dims) ^ _source_dims_dict.keys():
- # obtain dimensions that were not parsed source_dims_dict and attempt to match to a grid element
- if ds.sizes[dim] == grid.n_face:
- _source_dims_dict[dim] = "n_face"
- elif ds.sizes[dim] == grid.n_node:
- _source_dims_dict[dim] = "n_node"
- elif ds.sizes[dim] == grid.n_edge:
- _source_dims_dict[dim] = "n_edge"
+ for dim in set(ds.dims) ^ _source_dims_dict.keys():
+ # obtain dimensions that were not parsed source_dims_dict and attempt to match to a grid element
+ if ds.sizes[dim] == grid.n_face:
+ _source_dims_dict[dim] = "n_face"
+ elif ds.sizes[dim] == grid.n_node:
+ _source_dims_dict[dim] = "n_node"
+ elif ds.sizes[dim] == grid.n_edge:
+ _source_dims_dict[dim] = "n_edge"
- # Possible Issue: https://github.com/UXARRAY/uxarray/issues/610
+ # Possible Issue: https://github.com/UXARRAY/uxarray/issues/610
- # rename dimensions to follow the UGRID conventions
- ds = ds.swap_dims(_source_dims_dict)
+ # rename dimensions to follow the UGRID conventions
+ ds = ds.swap_dims(_source_dims_dict)
return ds
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index efaee93f4..f80ebf101 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -21,6 +21,7 @@
from uxarray.io._esmf import _read_esmf
from uxarray.io._vertices import _read_face_vertices
from uxarray.io._topology import _read_topology
+from uxarray.io._geos import _read_geos_cs
from uxarray.io.utils import _parse_grid_type
from uxarray.grid.area import get_all_face_area_from_coords
@@ -198,6 +199,8 @@ def from_dataset(
grid_ds, source_dims_dict = _read_mpas(dataset, use_dual=use_dual)
elif source_grid_spec == "ESMF":
grid_ds, source_dims_dict = _read_esmf(dataset)
+ elif source_grid_spec == "GEOS-CS":
+ grid_ds, source_dims_dict = _read_geos_cs(dataset)
elif source_grid_spec == "Shapefile":
raise ValueError("Shapefiles not yet supported")
else:
diff --git a/uxarray/io/_geos.py b/uxarray/io/_geos.py
new file mode 100644
index 000000000..af0de78ee
--- /dev/null
+++ b/uxarray/io/_geos.py
@@ -0,0 +1,59 @@
+import xarray as xr
+import numpy as np
+
+from uxarray.constants import INT_DTYPE
+from uxarray.conventions import ugrid
+
+
+def _read_geos_cs(in_ds: xr.Dataset):
+ """Reads and encodes a GEOS Cube-Sphere grid into the UGRID conventions.
+
+ https://gmao.gsfc.nasa.gov/gmaoftp/ops/GEOSIT_sample/doc/CS_Description_c180_v1.pdf
+ """
+ out_ds = xr.Dataset()
+
+ node_lon = in_ds["corner_lons"].values.ravel()
+ node_lat = in_ds["corner_lats"].values.ravel()
+
+ out_ds["node_lon"] = xr.DataArray(
+ data=node_lon, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LON_ATTRS
+ )
+
+ out_ds["node_lat"] = xr.DataArray(
+ data=node_lat, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LAT_ATTRS
+ )
+
+ if "lons" in in_ds:
+ face_lon = in_ds["lons"].values.ravel()
+ face_lat = in_ds["lats"].values.ravel()
+
+ out_ds["face_lon"] = xr.DataArray(
+ data=face_lon, dims=ugrid.FACE_DIM, attrs=ugrid.FACE_LON_ATTRS
+ )
+
+ out_ds["face_lat"] = xr.DataArray(
+ data=face_lat, dims=ugrid.FACE_DIM, attrs=ugrid.FACE_LAT_ATTRS
+ )
+
+ nf, nx, ny = in_ds["corner_lons"].shape
+
+ # generate indices for all corner nodes
+ idx = np.arange(nx * ny * nf, dtype=INT_DTYPE).reshape(nf, nx, ny)
+
+ # calculate indices of corner nodes for each face
+ tl = idx[:, :-1, :-1].reshape(-1)
+ tr = idx[:, :-1, 1:].reshape(-1)
+ bl = idx[:, 1:, :-1].reshape(-1)
+ br = idx[:, 1:, 1:].reshape(-1)
+
+ # Concatenate corner node indices for all faces
+ face_node_connectivity = np.column_stack((br, bl, tl, tr))
+
+ out_ds["face_node_connectivity"] = xr.DataArray(
+ data=face_node_connectivity,
+ dims=ugrid.FACE_NODE_CONNECTIVITY_DIMS,
+ attrs=ugrid.FACE_NODE_CONNECTIVITY_ATTRS,
+ )
+
+ # GEOS-CS does not return a source_dims_dict
+ return out_ds, None
diff --git a/uxarray/io/utils.py b/uxarray/io/utils.py
index 59846443b..9c836d47a 100644
--- a/uxarray/io/utils.py
+++ b/uxarray/io/utils.py
@@ -37,6 +37,9 @@ def _parse_grid_type(dataset):
mesh_type = "MPAS"
elif "maxNodePElement" in dataset.dims:
mesh_type = "ESMF"
+ elif all(key in dataset.sizes for key in ["nf", "YCdim", "XCdim"]):
+ # expected dimensions for a GEOS cube sphere grid
+ mesh_type = "GEOS-CS"
else:
raise RuntimeError("Could not recognize dataset format.")
return mesh_type