From bd0cae68e45518ff2452d5a69a38267c060f2811 Mon Sep 17 00:00:00 2001 From: mikebedington Date: Mon, 1 Jun 2026 20:06:39 +0200 Subject: [PATCH 1/5] Add reading cartesian coordinate files to FVCOM reader --- pyfvcom2/fvcom_reader.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/pyfvcom2/fvcom_reader.py b/pyfvcom2/fvcom_reader.py index 1fd358e..b036d5a 100644 --- a/pyfvcom2/fvcom_reader.py +++ b/pyfvcom2/fvcom_reader.py @@ -142,9 +142,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[0], epsg_code=coordinate_system[1]) return self._grid @property @@ -459,14 +459,22 @@ 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'][:] + # This assumes the coordinate projection is filled in and is either EPSG:xxxxx or just xxxxxx + coordinate_system = ['cartesian', self._metadata_dataset.CoordinateProjection.split(":")[-1]] + else: + x1 = self._metadata_dataset.variables['lon'][:] + x2 = self._metadata_dataset.variables['lat'][:] + coordinate_system = ['geographic', None] + 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. From bcd8a30f4f0ce0d7a0956558271a641084a25555 Mon Sep 17 00:00:00 2001 From: mikebedington Date: Fri, 5 Jun 2026 11:07:25 +0200 Subject: [PATCH 2/5] Add error handling for missing CoordinateProjection code --- pyfvcom2/fvcom_reader.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyfvcom2/fvcom_reader.py b/pyfvcom2/fvcom_reader.py index b036d5a..21792b7 100644 --- a/pyfvcom2/fvcom_reader.py +++ b/pyfvcom2/fvcom_reader.py @@ -464,7 +464,11 @@ def _extract_mesh_data(self) -> MeshData: x1 = self._metadata_dataset.variables['x'][:] x2 = self._metadata_dataset.variables['y'][:] # This assumes the coordinate projection is filled in and is either EPSG:xxxxx or just xxxxxx - coordinate_system = ['cartesian', self._metadata_dataset.CoordinateProjection.split(":")[-1]] + try: + coordinate_system = ['cartesian', self._metadata_dataset.CoordinateProjection.split(":")[-1]] + except AttributeError: + raise PyFVCOM2ValueError(f"Cartesian coordinates specificed but no CoordinateProjection global attribute provided") + else: x1 = self._metadata_dataset.variables['lon'][:] x2 = self._metadata_dataset.variables['lat'][:] From 6fea5542dc24c8c4c99535685e713436c6c9641b Mon Sep 17 00:00:00 2001 From: mikebedington Date: Sat, 6 Jun 2026 13:36:07 +0200 Subject: [PATCH 3/5] Slightly cleaner way of handling the projection code in fvcom_reader --- pyfvcom2/fvcom_reader.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/pyfvcom2/fvcom_reader.py b/pyfvcom2/fvcom_reader.py index 21792b7..198619d 100644 --- a/pyfvcom2/fvcom_reader.py +++ b/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 = {} @@ -144,7 +151,7 @@ def grid(self) -> Grid: if self._grid is None: mesh_data, coordinate_system = self._extract_mesh_data() sigma_data = self._extract_sigma_data() - self._grid = Grid(mesh_data, sigma_data, coordinate_system[0], epsg_code=coordinate_system[1]) + self._grid = Grid(mesh_data, sigma_data, coordinate_system, epsg_code=self._projection) return self._grid @property @@ -463,16 +470,16 @@ def _extract_mesh_data(self) -> MeshData: if getattr(self._metadata_dataset, 'CoordinateSystem', None) == 'Cartesian': x1 = self._metadata_dataset.variables['x'][:] x2 = self._metadata_dataset.variables['y'][:] - # This assumes the coordinate projection is filled in and is either EPSG:xxxxx or just xxxxxx - try: - coordinate_system = ['cartesian', self._metadata_dataset.CoordinateProjection.split(":")[-1]] - except AttributeError: - raise PyFVCOM2ValueError(f"Cartesian coordinates specificed but no CoordinateProjection global attribute provided") + 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', None] + coordinate_system = 'geographic' x3 = self._return_grid_variable_data('h')[:] open_bdy_node_lists = None From 01a579c25a7a392c7f84f51f6a3c040f4692fae0 Mon Sep 17 00:00:00 2001 From: mikebedington Date: Sat, 6 Jun 2026 13:36:35 +0200 Subject: [PATCH 4/5] Raise error if cartesian coordinates and no projection code available --- pyfvcom2/grid.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyfvcom2/grid.py b/pyfvcom2/grid.py index 714ed5a..94ed844 100644 --- a/pyfvcom2/grid.py +++ b/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 From 0af37e70b0326e64eab2df0c8a478f751df563b7 Mon Sep 17 00:00:00 2001 From: mikebedington Date: Tue, 23 Jun 2026 11:20:54 +0200 Subject: [PATCH 5/5] Bug fix, missing return argument type --- pyfvcom2/fvcom_reader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfvcom2/fvcom_reader.py b/pyfvcom2/fvcom_reader.py index 198619d..b2dcf6f 100644 --- a/pyfvcom2/fvcom_reader.py +++ b/pyfvcom2/fvcom_reader.py @@ -457,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: