#83: Add support for compound coordinate systems (#84)

.github/workflows/build.yml

@@ -17,7 +17,7 @@ jobs:
       fail-fast: false
       matrix:
         python-version: [3.9]
-        os: [ubuntu-latest, macos-latest, windows-latest]
+        os: [ubuntu-latest, windows-latest]
     runs-on: ${{ matrix.os }}
     steps:
       - uses: actions/checkout@v2

README.md

@@ -52,14 +52,14 @@ Currently, this tool supports the conversion of the following hydrodynamical mod
 ## Installation
 To install the converter, follow these steps:
 
-1. Install dependency manager [Poetry](https://python-poetry.org/docs/):
+1. Install dependency manager [Poetry](https://python-poetry.org/docs/)
 
 **Windows (PowerShell)**
 ```
 (Invoke-WebRequest -Uri https://install.python-poetry.org -UseBasicParsing).Content | python -
 ```
 
-**Linux, macOS, Windows (WSL)**
+**Linux, Windows (WSL)**
 ```
 curl -sSL https://install.python-poetry.org | python3 -
 ```
@@ -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:
+            crs1, crs2 = value.split("+", maxsplit=1)
+            return create_compound_crs(int(crs1), int(crs2))
+
+        return create_crs(int(value))
 
 class CrsShifting(BaseModel):
     """The configuration settings for shifting the coordinates."""

netcdf_to_gltf_converter/netcdf/parser.py

@@ -127,8 +127,8 @@ def _get_shift_values(shift_config: Union[ShiftType, CrsShifting], dataset: Data
 
         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))

netcdf_to_gltf_converter/preprocessing/crs.py

@@ -1,23 +1,55 @@
 import pyproj
 import pyproj.network
+from pyproj import CRS
+from pyproj.crs import CompoundCRS
 
 
-def create_crs_transformer(source_epsg: int, target_epsg: int) -> pyproj.Transformer:
+def create_crs_transformer(source_crs: CRS, target_crs: CRS) -> pyproj.Transformer:
     """Create a coordinate transformer that transforms the coordinates from a source coordinate system to a target coordinate system.
-    
+
     Args:
-        source_epsg (int): The EPSG code of the source coordinate system.
-        target_epsg (int): The EPSG code of the target coordinate system.
+        source_crs (CRS): The source coordinate system.
+        target_crs (CRS): The target coordinate system.
 
     Returns:
         pyproj.Transformer: The coordinate transformer
     """
     pyproj.network.set_network_enabled(True)
-
-    source_crs = pyproj.CRS.from_epsg(source_epsg)
-    target_crs = pyproj.CRS.from_epsg(target_epsg)
-
+
     transformer = pyproj.Transformer.from_crs(
         crs_from=source_crs, crs_to=target_crs, always_xy=True
     )
-    return transformer
+    return transformer
+
+
+def create_crs(epsg: int) -> CRS:
+    """Create a coordinate system from the given EPSG code.
+
+    Args:
+        epsg (int): The EPSG code to create the coordinate system from.
+
+    Returns:
+        CRS: The created coordinate system.
+    """
+    return CRS.from_epsg(epsg)
+
+
+def create_compound_crs(epsg_horizontal: int, epsg_vertical: int) -> CompoundCRS:
+    """Create a compound coordinate system from the given EPSG codes.
+
+    Args:
+        epsg_horizontal (int): EPSG code of the horizontal coordinate system..
+        epsg_vertical (int): EPSG code of the vertical coordinate system.
+
+    Returns:
+        CompoundCRS: The created compound coordinate system.
+    """
+    crs_horizontal = create_crs(epsg_horizontal)
+    crs_vertical = create_crs(epsg_vertical)
+
+    compound_crs = CompoundCRS(
+        name=f"{crs_horizontal.name} + {crs_vertical.name}",
+        components=[crs_horizontal, crs_vertical],
+    )
+
+    return compound_crs

tests/preprocessing/test_crs.py

@@ -1,10 +1,44 @@
-from netcdf_to_gltf_converter.preprocessing.crs import create_crs_transformer
+from pyproj import CRS
+from pyproj.crs import CompoundCRS
+
+from netcdf_to_gltf_converter.preprocessing.crs import (create_compound_crs,
+                                                        create_crs,
+                                                        create_crs_transformer)
 
 
 def test_create_crs_transformer():
-    transformer = create_crs_transformer(source_epsg=4979, target_epsg=7415)
-
+    source_epsg = 4979
+    target_epsg = 7415
+
+    source_crs = create_crs(source_epsg)
+    target_crs = create_crs(target_epsg)
+
+    transformer = create_crs_transformer(source_crs, target_crs)
+
     assert transformer is not None
     assert transformer.is_network_enabled == True
-    assert transformer.source_crs.name == 'WGS 84 (with axis order normalized for visualization)'
-    assert transformer.target_crs.name == 'Amersfoort / RD New + NAP height'
+    assert (
+        transformer.source_crs.name
+        == "WGS 84 (with axis order normalized for visualization)"
+    )
+    assert transformer.target_crs.name == "Amersfoort / RD New + NAP height"
+
+
+def test_create_crs():
+    epsg = 7415
+
+    crs = create_crs(epsg)
+
+    assert isinstance(crs, CRS)
+    assert crs.name == "Amersfoort / RD New + NAP height"
+
+
+def test_create_compound_crs():
+    epsg_horizontal = 32617
+    epsg_vertical = 5703
+
+    crs = create_compound_crs(epsg_horizontal, epsg_vertical)
+
+    assert isinstance(crs, CompoundCRS)
+    assert crs.name == "WGS 84 / UTM zone 17N + NAVD88 height"
+    assert crs.is_compound == True

tests/resources/d-hydro/westerschelde_config.json

@@ -5,8 +5,8 @@
 	"times_per_frame": 1,
 	"shift_coordinates": {
 		"crs_transformation": {
-			"source_epsg": 4979,
-			"target_epsg": 7415
+			"source_crs": 4979,
+			"target_crs": 7415
 		},
 		"shift_x": 3.9842655391511617,
 		"shift_y": 51.40756782573971,

tests/test_config.py

@@ -0,0 +1,43 @@
+from pyproj.crs import CompoundCRS
+
+from netcdf_to_gltf_converter.config import CrsTransformation
+
+
+class TestCrsTransformation:
+    def test_construction_with_epsg_integers(self):
+        source_epsg = 4979
+        target_epsg = 7415
+
+        crs_transformation = CrsTransformation(
+            source_crs=source_epsg, target_crs=target_epsg
+        )
+
+        assert crs_transformation.source_crs.name == "WGS 84"
+        assert crs_transformation.target_crs.name == "Amersfoort / RD New + NAP height"
+
+    def test_construction_with_epsg_strings(self):
+        source_epsg = "4979"
+        target_epsg = "7415"
+
+        crs_transformation = CrsTransformation(
+            source_crs=source_epsg, target_crs=target_epsg
+        )
+
+        assert crs_transformation.source_crs.name == "WGS 84"
+        assert crs_transformation.target_crs.name == "Amersfoort / RD New + NAP height"
+
+    def test_construction_with_compound_epsg_strings(self):
+        source_epsg = "4326+5773"
+        target_epsg = "32617+5703"
+
+        crs_transformation = CrsTransformation(
+            source_crs=source_epsg, target_crs=target_epsg
+        )
+
+        assert crs_transformation.source_crs.name == "WGS 84 + EGM96 height"
+        assert isinstance(crs_transformation.source_crs, CompoundCRS)
+        assert (
+            crs_transformation.target_crs.name
+            == "WGS 84 / UTM zone 17N + NAVD88 height"
+        )
+        assert isinstance(crs_transformation.target_crs, CompoundCRS)