#78: Add option to transform data values to another coordinate systyem

README.md

@@ -99,8 +99,8 @@ These steps will ensure that the converter is installed within a virtual environ
 - `shift_coordinates` (optional): A value indicating how to shift the coordinates of the data during conversion. When set to `min`, the converter will shift the coordinates such that the smallest x and y become the origin (0,0); variable values remain unchanged. It is also possible to provide custom shift values for the x- and y-coordinates and the variable values (z-coordinates):
 
   - `crs_transformation` (optional): The configuration settings for transforming the provided shift values from one coordinate system to another. The target coordinate system should be the coordinate system of the model.
-    - `source_epsg`: EPSG code of the source coordinate system.
-    - `target_epsg`: EPSG code of the target coordinate system.
+    - `source_crs`: EPSG code of the source coordinate system. Can also be a compound coordinate system, e.g. '32617+5703'
+    - `target_crs`: EPSG code of the target coordinate system.'Can also be a compound coordinate system, e.g. '32617+5703'
   - `shift_x`: A floating value containing the value that should be subtracted from all x-coordinates.
    - `shift_y`: A floating value containing the value that should be subtracted from all y-coordinates.
    - `shift_z`: A floating value containing the value that should be subtracted from all variable values (z-coordinates).

netcdf_to_gltf_converter/config.py

@@ -6,8 +6,10 @@
 from packaging.version import Version
 from pydantic import BaseModel as PydanticBaseModel
 from pydantic import Extra, root_validator, validator
+from pyproj import CRS
 from strenum import StrEnum
 
+from netcdf_to_gltf_converter.preprocessing.crs import create_compound_crs, create_crs
 from netcdf_to_gltf_converter.utils.validation import in_range
 
 Color = List[float]
@@ -108,11 +110,29 @@ class ShiftType(StrEnum):
 class CrsTransformation(BaseModel):
     """The configuration settings for transforming the coordinates."""
 
-    source_epsg: int
-    """int: EPSG code of the source coordinate system."""
+    source_crs: CRS
+    """CRS: The source coordinate system."""
 
-    target_epsg: int
-    """int: EPSG code of the target coordinate system."""
+    target_crs: CRS
+    """CRS: The target coordinate system."""
+
+    @validator("*", pre=True)
+    def validate_epsg(cls, value) -> CRS:
+        if isinstance(value, int):
+            return create_crs(value)
+
+        if isinstance(value, str):
+            return CrsTransformation._create_crs_from_string(value)
+
+        return value
+
+    @staticmethod
+    def _create_crs_from_string(value: str) -> CRS:
+        if "+" in value:
+            parts = value.split("+", maxsplit=1)
+            return create_compound_crs(int(parts[0]), int(parts[1]))
+
+        return create_crs(int(value))
 
 class CrsShifting(BaseModel):
     """The configuration settings for shifting the coordinates."""
@@ -210,6 +230,9 @@ class Config(AbstractJsonConfigFile, AbstractFileVersionFile):
     times_per_frame: int
     """int: The number of time steps per animation frame."""
 
+    vertical_crs_transformation: Optional[CrsTransformation]
+    """Optional[CrsTransformation]: Vertical coordinate transformation from one coordinate system to another."""
+
     shift_coordinates: Optional[Union[ShiftType, CrsShifting]]
     """Optional[Union[ShiftType, CrsShifting]]: The options how to shift the x-, y- and optionally z-coordinates. Typically used to create a reference point, such that an x-, y- and z-coordinate become the origin (0,0,0)."""
 

netcdf_to_gltf_converter/netcdf/netcdf_data.py

@@ -3,6 +3,7 @@
 from typing import List, Tuple
 
 import numpy as np
+from pyproj import CRS, Transformer
 import xarray as xr
 import xugrid.ugrid.connectivity as connectivity
 from xugrid.ugrid.conventions import X_STANDARD_NAMES, Y_STANDARD_NAMES
@@ -22,8 +23,7 @@ def get_coordinate_variables(data, standard_names: tuple) -> List[xr.DataArray]:
     return coord_vars
 
 
-
-class DataVariable():
+class DataVariable:
     """Class that serves as a wrapper object for an xarray.DataArray.
     The wrapper allows for easier retrieval of relevant data.
     """
@@ -47,7 +47,7 @@ def x_coords(self) -> np.ndarray:
             np.ndarray: A 1D np.ndarray of floats with shape (n, 1) where each row contains an x-coordinate.
         """
         return self._x_coords_var.values.flatten()
-    
+
     @property
     def y_coords(self) -> np.ndarray:
         """Get the y-coordinates for this data variable.
@@ -56,7 +56,7 @@ def y_coords(self) -> np.ndarray:
             np.ndarray: A 1D np.ndarray of floats with shape (n, 1) where each row contains a y-coordinate.
         """
         return self._y_coords_var.values.flatten()
-    
+
     @property
     def coordinates(self) -> np.ndarray:
         """Get the coordinates for this data variable.
@@ -84,27 +84,72 @@ def get_data_at_time(self, time_index: int) -> np.ndarray:
         Returns:
             np.ndarray: A 1D np.ndarray of floats.
         """
-        time_filter = {self._time_var.name : time_index}
+        time_filter = {self._time_var.name: time_index}
         return self._data_array.isel(**time_filter).values.flatten()
-
+
+    def _get_data_at_coordinate_index(self, coord_index: int) -> np.ndarray:
+        """Get the data values for this variable at the specified coordinate index.
+        Number of values is equal to the number of time steps.
+
+        Args:
+            coord_index (int): The coordinate index.
+
+        Returns:
+            np.ndarray: The values at the specified coordinate index.
+        """
+        filter = self._get_data_on_coordinate_filter(coord_index)
+        return self._data_array.isel(filter).values
+
+    def _set_data_at_coordinate_index(
+        self, coord_index: int, values: np.ndarray
+    ) -> None:
+        """Set the data values for this variable at the specified coordinate index.
+        Number of values should be equal to the number of time steps.
+
+        Args:
+            coord_index (int): The coordinate index.
+            values (np.ndarray): The new values.
+        """
+        filter = self._get_data_on_coordinate_filter(coord_index)
+        self._data_array.loc[filter] = values
+
+    def transform_data_values(self, crs_transformer: Transformer) -> None:
+        """Transform all data values (z-coordinates) for each xy-coordinate and for each time step.
+
+        Args:
+            crs_transformer (Transformer): The transformer to handle the coordinates' transformation.
+        """
+        for coord_index in range(0, len(self.coordinates)):                
+            z_coords = self._get_data_at_coordinate_index(coord_index)
+            x_coords = np.full(z_coords.size, self.x_coords[coord_index])
+            y_coords = np.full(z_coords.size, self.y_coords[coord_index])
+
+            _, _, z_coords_transformed = crs_transformer.transform(x_coords, y_coords, z_coords)
+            self._set_data_at_coordinate_index(coord_index, z_coords_transformed)
+
     @property
     def min(self) -> float:
         """Get the minimum value for this variable across all dimensions.
-        
+
         Returns:
             float: The minimum variable value.
         """
         return self._data_array.min().values
-    
+
     @property
     def max(self) -> float:
         """Get the maximum value for this variable across all dimensions.
-        
+
         Returns:
             float: The maximum variable value.
         """
         return self._data_array.max().values
-
+
+    def _get_data_on_coordinate_filter(self, coord_index: int):
+        filter = {self._x_coords_var.dims[0]: coord_index}
+        return filter
+
+
 class DatasetBase(ABC):
     """Class that serves as a wrapper object for an xarray.Dataset.
     The wrapper allows for easier retrieval of relevant data.
@@ -183,7 +228,7 @@ def get_variable(self, variable_name: str) -> DataVariable:
     def _raise_if_not_in_dataset(self, name: str):
         if name not in self._dataset:
             raise ValueError(f"Variable with name {name} does not exist in dataset.")
-    
+
     @property
     @abstractmethod
     def face_node_connectivity(self) -> np.ndarray:
@@ -222,7 +267,20 @@ def fill_value(self) -> int:
             int: Integer with the fill value.
         """
         pass
-
+
+    @abstractmethod
+    def transform_vertical_coordinate_system(
+        self, source_crs: CRS, target_crs: CRS, variables: List[str]
+    ) -> None:
+        """Transform the vertical coordinates to another coordinate system.
+
+        Args:
+            source_crs (CRS): The source coordinate system.
+            target_crs (CRS): The target coordinate system.
+            variables (List[str]): The names of the variables for which to transform the values.
+        """
+        pass
+
     @abstractmethod
     def shift_coordinates(self, shift: Vec3, variables: List[str]) -> None:
         """
@@ -236,9 +294,11 @@ def shift_coordinates(self, shift: Vec3, variables: List[str]) -> None:
             variables (List[str]): The names of the variables for which to shift the values.
         """
         pass
-    
+
     @abstractmethod
-    def scale_coordinates(self, scale_horizontal: float, scale_vertical: float, variables: List[str]) -> None:
+    def scale_coordinates(
+        self, scale_horizontal: float, scale_vertical: float, variables: List[str]
+    ) -> None:
         """
         Scale the x- and y-coordinates and the data values, with the scaling factors that are specified.
         The original data set is updated with the new coordinates.
@@ -249,7 +309,6 @@ def scale_coordinates(self, scale_horizontal: float, scale_vertical: float, vari
             variables (List[str]): The names of the variables to scale.
         """
         pass
-
 
     def triangulate(self):
         """Triangulate the provided grid.
@@ -262,6 +321,8 @@ def triangulate(self):
             self.face_node_connectivity, self.fill_value
         )
         self.set_face_node_connectivity(face_node_connectivity)
-        
+
     def _log_grid_bounds(self, bounds: Tuple[float, float, float, float]):
-        logging.info(f"Grid bounds: {bounds[0]} (min x), {bounds[1]} (min y), {bounds[2]} (max x), {bounds[3]} (max y)")
+        logging.info(
+            f"Grid bounds: {bounds[0]} (min x), {bounds[1]} (min y), {bounds[2]} (max x), {bounds[3]} (max y)"
+        )

netcdf_to_gltf_converter/netcdf/parser.py

@@ -103,8 +103,15 @@ def _transform_grid(config: Config, dataset: DatasetBase):
         variables = [var.name for var in config.variables]
         Parser._log_variable_values(dataset, variables)
 
+        if config.vertical_crs_transformation:
+            logging.info(f"TRANSFORM data values from {config.vertical_crs_transformation.source_crs.name} to {config.vertical_crs_transformation.target_crs.name}")
+            dataset.transform_vertical_coordinate_system(config.vertical_crs_transformation.source_crs, 
+                                                         config.vertical_crs_transformation.target_crs, 
+                                                         variables)
+            Parser._log_variable_values(dataset, variables)
+
         if config.shift_coordinates:
-            shift = Parser._get_shift_values(config.shift_coordinates, dataset)
+            shift = Parser._get_shift_values(config, dataset)
 
             logging.info(f"SHIFT model coordinates with: {shift.x} (x), {shift.y} (y), {shift.z} (z)")
             dataset.shift_coordinates(shift, variables)
@@ -121,17 +128,24 @@ def _log_variable_values(dataset: DatasetBase, variables: List[str]):
             logging.info(f"Variable values for '{variable_name}': {variable.min} (min), {variable.max} (max)")
 
     @staticmethod
-    def _get_shift_values(shift_config: Union[ShiftType, CrsShifting], dataset: DatasetBase) -> Vec3:
+    def _get_shift_values(config: Config, dataset: DatasetBase) -> Vec3:
+        shift_config = config.shift_coordinates
         if shift_config == ShiftType.MIN:
             return Vec3(dataset.min_x, dataset.min_y, 0.0)
 
         if isinstance(shift_config, CrsShifting):
             if shift_config.crs_transformation:
-                transformer = create_crs_transformer(shift_config.crs_transformation.source_epsg,
-                                                     shift_config.crs_transformation.target_epsg)
+                transformer = create_crs_transformer(shift_config.crs_transformation.source_crs,
+                                                     shift_config.crs_transformation.target_crs)
 
                 logging.info(f"Transforming shift values from {transformer.source_crs.name} to {transformer.target_crs.name}")
-                return Vec3(*transformer.transform(shift_config.shift_x, shift_config.shift_y, shift_config.shift_z))
+                shift_x, shift_y, shift_z = transformer.transform(shift_config.shift_x, shift_config.shift_y, shift_config.shift_z)
+
+                vert_crs_transformation = config.vertical_crs_transformation
+                if vert_crs_transformation and vert_crs_transformation.target_crs.equals(shift_config.crs_transformation.source_crs):
+                    shift_z = shift_config.shift_z
+
+                return Vec3(shift_x, shift_y, shift_z)
 
             return Vec3(shift_config.shift_x, shift_config.shift_y, shift_config.shift_z)
 

netcdf_to_gltf_converter/netcdf/ugrid/ugrid_data.py

@@ -1,11 +1,13 @@
 from typing import List
 
 import numpy as np
+from pyproj import CRS
 import xarray as xr
 import xugrid as xu
 
 from netcdf_to_gltf_converter.data.vector import Vec3
 from netcdf_to_gltf_converter.netcdf.netcdf_data import DatasetBase
+from netcdf_to_gltf_converter.preprocessing.crs import create_crs_transformer
 
 
 class UgridDataset(DatasetBase):
@@ -22,7 +24,7 @@ def __init__(self, dataset: xr.Dataset) -> None:
         super().__init__(dataset)
         self._ugrid_data_set = xu.UgridDataset(dataset)  
         self._grid = self._get_ugrid2d()
-        super()._log_grid_bounds(self._grid.bounds)
+        self._update()
 
     @property
     def min_x(self) -> float:
@@ -78,6 +80,19 @@ def fill_value(self) -> int:
             int: Integer with the fill value.
         """
         return self._grid.fill_value
+
+    def transform_vertical_coordinate_system(self, source_crs: CRS, target_crs: CRS, variables: List[str]) -> None:
+        """Transform the vertical coordinates to another coordinate system.
+
+        Args:
+            source_crs (CRS): The source coordinate system.
+            target_crs (CRS): The target coordinate system.
+            variables (List[str]): The names of the variables for which to transform the values.
+        """
+        transformer = create_crs_transformer(source_crs, target_crs)
+        for variable_name in variables:
+            variable = self.get_variable(variable_name)
+            variable.transform_data_values(transformer)
 
     def shift_coordinates(self, shift: Vec3, variables: List[str]) -> None:
         """

netcdf_to_gltf_converter/netcdf/xbeach/xbeach_data.py

@@ -1,6 +1,7 @@
 from typing import List, Tuple
 
 import numpy as np
+from pyproj import CRS
 import xarray as xr
 from xugrid.ugrid.conventions import X_STANDARD_NAMES, Y_STANDARD_NAMES
 
@@ -118,7 +119,17 @@ def fill_value(self) -> int:
             int: Integer with the fill value.
         """
         return -1
-
+
+    def transform_vertical_coordinate_system(self, source_crs: CRS, target_crs: CRS, variables: List[str]) -> None:
+        """Transform the vertical coordinates to another coordinate system.
+
+        Args:
+            source_crs (CRS): The source coordinate system.
+            target_crs (CRS): The target coordinate system.
+            variables (List[str]): The names of the variables for which to transform the values.
+        """
+        raise NotImplementedError()
+
     def shift_coordinates(self, shift: Vec3, variables: List[str]) -> None:
         """
         Shift the x- and y-coordinates and the variable values in the data set with the provided shift values.