use CRS.to_2d() for transforming with compound CRS's

GeodetischeInfrastructuur · Nov 19, 2024 · 9624758 · 9624758
1 parent 475ea44
commit 9624758
Show file tree

Hide file tree

Showing 10 changed files with 91 additions and 299 deletions.
diff --git a/src/coordinate_transformation_api/assets/config.yaml b/src/coordinate_transformation_api/assets/config.yaml
@@ -469,6 +469,10 @@ EPSG:7415:
   uri: http://www.opengis.net/def/crs/EPSG/0/7415
 NSGI:Bonaire_DPnet_KADpeil:
   exclude-transformations:
+    - EPSG:7789
+    - EPSG:7912
+    - EPSG:4979
+    - OGC:CRS84h
     - EPSG:28992
     - EPSG:3034
     - EPSG:3035
@@ -924,6 +928,8 @@ NSGI:St_Eustatius2020_GEOGRAPHIC_3D:
 EPSG:7912:
   exclude-transformations:
     - EPSG:9289
+    - NSGI:Bonaire_DPnet_KADpeil
+    - NSGI:Bonaire_DPnet
   uri: http://www.opengis.net/def/crs/EPSG/0/7912
 EPSG:4979:
   exclude-transformations:
@@ -1056,4 +1062,6 @@ NSGI:St_Eustatius2020_GEOCENTRIC:
 EPSG:7789:
   exclude-transformations:
     - EPSG:9289
+    - NSGI:Bonaire_DPnet_KADpeil
+    - NSGI:Bonaire_DPnet
   uri: http://www.opengis.net/def/crs/EPSG/0/7789
diff --git a/src/coordinate_transformation_api/crs_transform.py b/src/coordinate_transformation_api/crs_transform.py
@@ -35,7 +35,6 @@
 )
 from coordinate_transformation_api.models import (
     TransformationNotPossibleError,
-    UnknownCrsError,
 )
 from coordinate_transformation_api.types import CoordinatesType, ShapelyGeometry
 
@@ -196,18 +195,6 @@ def check_axis(s_crs: CRS, t_crs: CRS) -> None:
 
 
 def get_transformer(source_crs: CRS, target_crs: CRS, epoch: float | None) -> Transformer:  # quit
-    check_axis(source_crs, target_crs)
-
-    if exclude_transformation(
-        "{}:{}".format(*source_crs.to_authority()),
-        "{}:{}".format(*target_crs.to_authority()),
-    ):
-        raise TransformationNotPossibleError(
-            "{}:{}".format(*source_crs.to_authority()),
-            "{}:{}".format(*target_crs.to_authority()),
-            "Transformation Excluded",
-        )
-
     # Get available transformer through TransformerGroup
     # TODO check/validate if always_xy=True is correct
     tfg = transformer.TransformerGroup(source_crs, target_crs, allow_ballpark=False, always_xy=True)
@@ -243,18 +230,6 @@ def get_transformer(source_crs: CRS, target_crs: CRS, epoch: float | None) -> Tr
     return tfg.transformers[0]
 
 
-def get_individual_crs_from_compound(compound_crs: CRS) -> tuple[CRS, CRS]:
-    horizontal = compound_crs
-    vertical = compound_crs
-    for crs in compound_crs.sub_crs_list:
-        if len(crs.axis_info) == HORIZONTAL_AXIS_LENGTH:
-            horizontal = crs
-        elif len(crs.axis_info) == VERTICAL_AXIS_LENGTH:
-            vertical = crs
-
-    return (horizontal, vertical)
-
-
 def build_input_coord(coord: CoordinatesType, epoch: float | None) -> CoordinatesType:
     # When 2D input is given with an epoch we need to add a height. So pyproj knows to
     # that the epoch is an epoch and not the height, without this intervention the epoch
@@ -311,45 +286,35 @@ def get_transform_crs_fun(
     if precision is None:
         precision = get_precision(target_crs)
 
-    # We need to do something special for transformation targetting a Compound CRS of 2D coordinates with another height system, like NAP or a LAT height
+    check_axis(source_crs, target_crs)
+    if exclude_transformation(
+        "{}:{}".format(*source_crs.to_authority()),
+        "{}:{}".format(*target_crs.to_authority()),
+    ):
+        raise TransformationNotPossibleError(
+            "{}:{}".format(*source_crs.to_authority()),
+            "{}:{}".format(*target_crs.to_authority()),
+            "Transformation Excluded",
+        )
+
+    # We need to do something special for transformation involving a Compound CRS of 2D coordinates with another height system, like NAP or a LAT height
     # - RD + NAP (EPSG:7415)
     # - ETRS89 + NAP (EPSG:9286)
     # - ETRS89 + LAT-NL (EPSG:9289)
-    # These transformations need to be splitted in a horizontal and vertical transformation.
-    if target_crs is not None and source_crs is not target_crs and target_crs.is_compound:
-        check_axis(source_crs, target_crs)
-
-        target_crs_horizontal, target_crs_vertical = get_individual_crs_from_compound(target_crs)
-
-        if source_crs is not None and source_crs.is_compound:
-            (
-                source_crs_horizontal,
-                source_crs_vertical,
-            ) = get_individual_crs_from_compound(source_crs)
-        else:
-            source_crs_horizontal = source_crs
-            source_crs_vertical = source_crs
-
-        h_transformer = get_transformer(source_crs_horizontal, target_crs_horizontal, epoch)
-
-        # Not all transformation that are possible are defined
-        # When no transformation is found we fall back on the original COMPOUND CRS
-        # Issue is that in some case a transformation is found but not the correct one
-        # These we identify by the laking of a AUTO:CODE, because all our CRS should be
-        # coded. These are also defaulted to the original COMPOUND CRS.
-        try:
-            v_transformer = get_transformer(source_crs_vertical, target_crs_vertical, epoch)
-            if v_transformer.source_crs is not None and v_transformer.source_crs.to_authority() is None:
-                raise UnknownCrsError()  # empty error, we catch it the line below
-        except (TransformationNotPossibleError, UnknownCrsError):
-            v_transformer = get_transformer(source_crs, target_crs, epoch)
+    # These transformations need to be splitted in a horizontal and vertical transformation (vertical transformation actually attempts the 3d transformation).
+    if target_crs is not None and source_crs is not target_crs and (target_crs.is_compound or source_crs.is_compound):
+        target_crs_horizontal = target_crs.to_2d()
+        h_transformer = get_transformer(source_crs, target_crs_horizontal, epoch)
+        v_transformer = get_transformer(
+            source_crs, target_crs, epoch
+        )  # this will do the 3d transformation that might fail, in that case Z/H value is dropped
 
         # note transformers are injected in transform_compound_crs so they are instantiated only once
         _transform_compound_crs = partial(transform_compound_crs, h_transformer, v_transformer, precision, epoch)
         return _transform_compound_crs
     else:
         transformer = get_transformer(source_crs, target_crs, epoch)
-        # note transformer is injected in transform_compound_crs is instantiated once
+        # note transformer is injected in transform_crs is instantiated once
         # creating transformers is expensive
         _transform_crs = partial(transform_crs, transformer, precision, epoch)
         return _transform_crs
@@ -379,7 +344,9 @@ def transform_compound_crs(
 
     output_2d = Position2D(*h[:2])
     output: Position = output_2d
-    if len(v) == THREE_DIMENSIONAL and not math.isinf(v[2]):
+    if len(v) >= THREE_DIMENSIONAL and not math.isinf(
+        v[2]
+    ):  # note len(v) can be larger than three when epoch is supplied
         output = Position3D(*output_2d, v[2])
     else:
         # height coordinate dropped, since v[2] not added
@@ -421,6 +388,9 @@ def transform_crs(
 
     _output = tuple(map(_round_h, transformer.transform(*input)[0:dim]))
 
+    # - 1 van de twee compound
+    # -
+
     output_2d = Position2D(*_output[:2])
     output: Position = output_2d
     if len(_output) >= THREE_DIMENSIONAL:

diff --git a/tests/data/bonaire_validation_data.csv b/tests/data/bonaire_validation_data.csv
@@ -78,16 +78,16 @@ EPSG:7912,EPSG:7789,(inf inf inf inf)
 EPSG:7912,EPSG:7912,(inf inf inf inf)
 EPSG:7912,EPSG:4979,(inf inf inf inf)
 EPSG:7912,OGC:CRS84h,(inf inf inf inf)
-EPSG:4979,NSGI:Bonaire_DPnet_KADpeil,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,NSGI:Bonaire_DPnet,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,NSGI:Bonaire2004_GEOCENTRIC,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,NSGI:Bonaire2004_GEOGRAPHIC_2D,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,NSGI:Bonaire2004_GEOGRAPHIC_3D,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,EPSG:32619,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,EPSG:7789,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,EPSG:7912,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,EPSG:4979,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
-EPSG:4979,OGC:CRS84h,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
+EPSG:4979,NSGI:Bonaire_DPnet_KADpeil,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,NSGI:Bonaire_DPnet,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,NSGI:Bonaire2004_GEOCENTRIC,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,NSGI:Bonaire2004_GEOGRAPHIC_2D,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,NSGI:Bonaire2004_GEOGRAPHIC_3D,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,EPSG:32619,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,EPSG:7789,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,EPSG:7912,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,EPSG:4979,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
+EPSG:4979,OGC:CRS84h,(-68.2537259371291 12.14970719200454 -40.252278505824506 2000.0)
 OGC:CRS84h,NSGI:Bonaire_DPnet_KADpeil,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
 OGC:CRS84h,NSGI:Bonaire_DPnet,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)
 OGC:CRS84h,NSGI:Bonaire2004_GEOCENTRIC,(-68.2537259371291 12.149707191615798 -40.2522831344977 2000.0)