Fvcom reader to handle non-geographic coordinates

pyfvcom2/fvcom_reader.py

@@ -30,11 +30,13 @@ class FVCOMReader:
     DAYS_PER_MILLISECOND = 1.0 / (1000.0 * 60.0 * 60.0 * 24.0)
 
     def __init__(self,
-                 file_paths: Union[str, List[str]]):
+                 file_paths: Union[str, List[str]],
+                 projection: str=None):
         """Initialize the FVCOMReader with the path to the netCDF file.
 
         Args:
             file_paths (str, list): Path to the FVCOM netCDF file.
+            projection (str, optional): EPSG projection code, as either EPSG:xxxx or xxxx
         """
         # Handle single file path or list of file paths
         if isinstance(file_paths, str):
@@ -45,13 +47,18 @@ def __init__(self,
         # Load only the first file initially for metadata and time-independent data
         print(f'Accessing FVCOM metadata from: {self.file_paths[0]}')
         self._metadata_dataset = Dataset(self.file_paths[0])
-
         self._grid = None  # Lazy initialization
         self._z_level_cache = {}  # Cache for time-dependent z levels
 
         # Build the time index mapping for multiple files
         self._build_time_index_mapping()
 
+        if projection is None and hasattr(self._metadata_dataset,'CoordinateProjection'):
+            self._projection = self._metadata_dataset.CoordinateProjection.split(':')[-1]
+        else:
+            self._projection = projection.split(':')[-1]
+
+
     def _build_time_index_mapping(self):
         """Build a mapping from datetime to (file_path, local_time_index)"""
         self._time_to_file_map = {}
@@ -142,9 +149,9 @@ def grid(self) -> Grid:
             Grid: The grid object containing mesh structure and open boundaries.
         """
         if self._grid is None:
-            mesh_data = self._extract_mesh_data()
+            mesh_data, coordinate_system = self._extract_mesh_data()
             sigma_data = self._extract_sigma_data()
-            self._grid = Grid(mesh_data, sigma_data, "geographic")
+            self._grid = Grid(mesh_data, sigma_data, coordinate_system, epsg_code=self._projection)
         return self._grid
 
     @property
@@ -450,7 +457,7 @@ def get_interpolation_coordinates(self, horizontal_position: str,
                                                        vertical_coordinate_system=vertical_coordinate_system,
                                                        dates=dates)
 
-    def _extract_mesh_data(self) -> MeshData:
+    def _extract_mesh_data(self) -> tuple[MeshData, str]:
         """Extract mesh data from FVCOM output file.
 
         Returns:
@@ -459,14 +466,26 @@ def _extract_mesh_data(self) -> MeshData:
         # Extract basic mesh components
         nodes = np.arange(1, self._metadata_dataset.dimensions['node'].size+1) # TBC zero based indexing kept here.
         triangles = self._metadata_dataset.variables['nv'][:].T - 1  # Convert to 0-based indexing, transpose to (n_elem, 3)
-        x1 = self._metadata_dataset.variables['lon'][:]
-        x2 = self._metadata_dataset.variables['lat'][:]
-        x3 = self._return_grid_variable_data('h')[:]
+
+        if getattr(self._metadata_dataset, 'CoordinateSystem', None) == 'Cartesian':
+            x1 = self._metadata_dataset.variables['x'][:]
+            x2 = self._metadata_dataset.variables['y'][:]
+            if self._projection is not None:
+                # This assumes the coordinate projection in the netCDF is either EPSG:xxxxx or just xxxxxx
+                coordinate_system = 'cartesian'
+            else:
+                raise PyFVCOM2ValueError(f"Cartesian coordinates specificed but no CoordinateProjection provided in file or passed")
+
+        else:
+            x1 = self._metadata_dataset.variables['lon'][:]
+            x2 = self._metadata_dataset.variables['lat'][:]
+            coordinate_system = 'geographic'
 
+        x3 = self._return_grid_variable_data('h')[:]
         open_bdy_node_lists = None
         bdy_types = None
 
-        return MeshData(triangles, nodes, x1, x2, x3, bdy_types, open_bdy_node_lists)
+        return MeshData(triangles, nodes, x1, x2, x3, bdy_types, open_bdy_node_lists), coordinate_system
 
     def _extract_sigma_data(self) -> SigmaData:
         """Extract sigma coordinate data from FVCOM output file.

pyfvcom2/grid.py

@@ -176,7 +176,10 @@ def __init__(
         if coordinate_system == "cartesian":
             self._x = mesh_data.x1
             self._y = mesh_data.x2
-            self.epsg_code = epsg_code
+            if epsg_code is not None:
+                self.epsg_code = epsg_code
+            else:
+                raise PyFVCOM2ValueError("epsg code must be provided for cartesian coordinates")
             self._lon, self._lat = lonlat_from_utm(self._x, self._y, epsg_code)
         elif coordinate_system == "geographic":
             self._lon = mesh_data.x1