Merge pull request #154 from IGNF/tf_proj_data

New functionality transform projection of points
IGNF · Sep 30, 2024 · c45ae90 · c45ae90
2 parents eb15890 + dffad7d
commit c45ae90
Show file tree

Hide file tree

Showing 82 changed files with 1,542 additions and 290 deletions.
diff --git a/.github/workflows/test_mac.yml b/.github/workflows/test_mac.yml
@@ -35,11 +35,11 @@ jobs:
 
     - name: Run regression test opk to opk
       shell: bash -el {0}
-      run: python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run: python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test opk to opk with multi epsg
       shell: bash -el {0}
-      run : python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run : python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test control opk
       shell: bash -el {0}
@@ -67,8 +67,12 @@ jobs:
 
     - name: Run regression test spaceresection opk
       shell: bash -el {0}
-      run: python ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC -ou degree -oa True
+      run: python ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC --ou degree --oa True
 
     - name: Run regression test opk to conical
       shell: bash -el {0}
       run: python ./borea_tools/opk_to_conl.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -w ./test/tmp/
+
+    - name: Run transform projextion of points file
+      shell: bash -el {0}
+      run: python ./borea_tools/transform_proj_points.py -g ./dataset/GCP_test.app -l PTXYZ -e 2154 -y ./dataset/fr_ign_RAF20.tif --oe 4326 --oz height -n GCP_4326 -w ./test/tmp
diff --git a/.github/workflows/test_ubu.yml b/.github/workflows/test_ubu.yml
@@ -33,10 +33,10 @@ jobs:
       run: ./ci/test_example.sh
 
     - name: Run regression test opk to opk
-      run: python3 ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run: python3 ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test opk to opk with multi epsg
-      run : python3 ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run : python3 ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test control opk
       run: python3 ./borea_tools/opk_control.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./dataset/terrain_test.mes -g ./dataset/GCP_test.app -d 13 -l PTXYH -p inter -w ./test/tmp/
@@ -57,7 +57,10 @@ jobs:
       run: python3 ./borea_tools/ptfile_image_to_world.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./dataset/terrain_test.mes -p square -n Coor3d_pt_image -w ./test/tmp/
 
     - name: Run regression test spaceresection opk
-      run: python3 ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC -ou degree -oa True
+      run: python3 ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC --ou degree --oa True
 
     - name: Run regression test opk to conical
       run: python3 ./borea_tools/opk_to_conl.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -w ./test/tmp/
+
+    - name: Run transform projextion of points file
+      run: python3 ./borea_tools/transform_proj_points.py -g ./dataset/GCP_test.app -l PTXYZ -e 2154 -y ./dataset/fr_ign_RAF20.tif --oe 4326 --oz height -n GCP_4326 -w ./test/tmp
diff --git a/.github/workflows/test_win.yml b/.github/workflows/test_win.yml
@@ -35,11 +35,11 @@ jobs:
 
     - name: Run regression test conv opk
       shell: bash -el {0}
-      run: python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run: python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test opk to opk with multi epsg
       shell: bash -el {0}
-      run : python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC -ou radian -oa False
+      run : python ./borea_tools/opk_to_opk.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 --geoc 4964 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -n Test -w ./test/tmp/ -o NXYHPOKC --ou radian --oa False
 
     - name: Run regression test control opk
       shell: bash -el {0}
@@ -67,8 +67,12 @@ jobs:
 
     - name: Run regression test spaceresection opk
       shell: bash -el {0}
-      run: python ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC -ou degree -oa True
+      run: python ./borea_tools/spaceresection_opk.py -p 825439 6289034 1500 -d height -c ./dataset/Camera1.txt -e 2154 -y ./dataset/fr_ign_RAF20.tif -m ./dataset/MNT_France_25m_h_crop.tif --fm height -t ./test/data/dataset2/all_liaisons2.mes -g ./test/data/dataset2/all_liaisons2_world.mes -l PXYH -n SpaceResection -w ./test/tmp/ -o NXYZOPKC --ou degree --oa True
 
     - name: Run regression test opk to conical
       shell: bash -el {0}
       run: python ./borea_tools/opk_to_conl.py -r ./dataset/23FD1305_alt_test.OPK -i NXYZOPKC -f 2 -e 2154 -y ./dataset/fr_ign_RAF20.tif -c ./dataset/Camera1.txt -m ./dataset/MNT_France_25m_h_crop.tif --fm height -w ./test/tmp/
+
+    - name: Run transform projextion of points file
+      shell: bash -el {0}
+      run: python ./borea_tools/transform_proj_points.py -g ./dataset/GCP_test.app -l PTXYZ -e 2154 -y ./dataset/fr_ign_RAF20.tif --oe 4326 --oz height -n GCP_4326 -w ./test/tmp
diff --git a/README.md b/README.md
@@ -16,6 +16,7 @@ Why Borea? B for Box and orea is a back slang of aero.
 * Transforms coordinates file terrain from image: [borea_tools/docs_tools/README_ptfile_image_to_world.md](./borea_tools/docs_tools/README_ptfile_image_to_world.md)
 * Transforms coordinates file image from terrain: [borea_tools/docs_tools/README_ptfile_world_to_image.md](./borea_tools/docs_tools/README_ptfile_world_to_image.md)
 * Calculates opk by space resection: [borea_tools/docs_tools/README_spaceresection_opk.md](./borea_tools/docs_tools/README_spaceresection_opk.md)
+* Transform projection of points file: [borea_tools/docs_tools/README_transform_proj_points.md](./borea_tools/docs_tools/README_transform_proj_points.md)
 * Python lib: [README_borea_lib.md](./README_borea_lib.md)
 
 ## Dependency

diff --git a/README_borea_lib.md b/README_borea_lib.md
@@ -43,6 +43,7 @@ Some tools are already implemented in the library:
 * Transforms coordinates file terrain from image: `ptfile-image-to-world -h` [doc](https://github.com/IGNF/Borea/tree/main/borea_tools/docs_tools/README_ptfile_image_to_world.md)
 * Transforms coordinates file image from terrain: `ptfile-world-to-image -h` [doc](https://github.com/IGNF/Borea/tree/main/borea_tools/docs_tools/README_ptfile_world_to_image.md)
 * Calculates opk by space resection: `spaceresection-opk -h` [doc](https://github.com/IGNF/Borea/tree/main/borea_tools/docs_tools/README_spaceresection_opk.md)
+* Transform projection of points file: `transform-proj-points -h` [doc](./borea_tools/docs_tools/README_transform_proj_points.md)
 
 ## Read data and instantiate worksite
 

diff --git a/borea/datastruct/shot.py b/borea/datastruct/shot.py
@@ -185,6 +185,12 @@ def set_order_axe(self, order_axe: str) -> None:
             self.ori_shot = -(R.from_euler("x", np.pi) *
                               R.from_matrix(self.mat_rot)).as_euler(order_xyz, degrees=unitori)
 
+    def set_proj_pos(self) -> None:
+        """
+        Change the projection of the position shot with the epsg_output.
+        """
+        self.pos_shot = np.array(ProjEngine().tf.proj_to_proj_out(*self.pos_shot))
+
     def get_z_remove_scale_factor(self) -> float:
         """
         Return Z after removing the scale factor. The Z of the object is NOT modified.

diff --git a/borea/datastruct/workdata.py b/borea/datastruct/workdata.py
@@ -33,6 +33,7 @@ def __init__(self, name: str) -> None:
         self.type_z_data = None
         self.type_z_shot = None
         self.approxeucli = False
+        self.epsg_output = False
 
     # pylint: disable-next=too-many-arguments too-many-positional-arguments
     def add_shot(self, name_shot: str, pos_shot: np.ndarray,
@@ -59,24 +60,28 @@ def add_shot(self, name_shot: str, pos_shot: np.ndarray,
                                      linear_alteration=linear_alteration,
                                      order_axe=order_axe)
 
-    def set_proj(self, epsg: list, path_geoid: list = None) -> None:
+    def set_proj(self, epsg: list, path_geoid: list = None, epsg_output: int = None) -> None:
         """
         Setup a projection system to the worksite.
 
         Args:
-            epsg (list): Code epsg of the porjection ex: 2154.
+            epsg (list): Code epsg of the projection ex: [2154].
             path_geoid (str): List of GeoTIFF which represents the geoid in grid form.
+            epsg_output (int): Code epsg of the output projection. If you want to change.
         """
         ProjEngine.clear()
-        try:  # Check if the epsg exist
-            for idepsg in epsg:
-                if idepsg:
+        for idepsg in [*epsg, epsg_output]:
+            if idepsg:
+                try:  # Check if the epsg exist
                     _ = CRS.from_epsg(idepsg)
                     del _
-        except exceptions.CRSError as e_info:
-            raise exceptions.CRSError(f"Your EPSG:{epsg} doesn't exist in pyproj.") from e_info
+                except exceptions.CRSError as e_info:
+                    raise exceptions.CRSError(f"Your EPSG:{epsg}"
+                                              " doesn't exist in pyproj.") from e_info
+        if epsg_output:
+            self.epsg_output = True
 
-        ProjEngine().set_epsg(epsg, path_geoid)
+        ProjEngine().set_epsg(epsg, path_geoid, epsg_output)
 
     # pylint: disable-next=too-many-arguments too-many-positional-arguments
     def add_camera(self, name_camera: str, ppax: float, ppay: float,

diff --git a/borea/format/conl.py b/borea/format/conl.py
@@ -33,19 +33,15 @@ class Conl:
     cam: Camera
     proj: str
 
-    def save_conl(self, path_conical: Path, linalt: bool = True) -> None:
+    def save_conl(self, path_conical: Path) -> None:
         """
         Save the shot as light conical file.
 
         Args:
             path_conical (Path): path to the light conical file.
-            linalt (bool): If you want z shot corrected by linear alteration.
         """
         date_now = datetime.now()
 
-        # Scale factor correction
-        self.shot.set_linear_alteration(linalt)
-
         # creation XML
         ori = ET.Element("orientation", {})
         ET.SubElement(ori, "lastmodificationbylibori",

diff --git a/borea/format/rpc.py b/borea/format/rpc.py
@@ -22,7 +22,7 @@ class Rpc:
     def __init__(self) -> None:
         self.param_rpc = {}
         self.fact_rpc = None
-        self.output_epsg = None
+        self.output_epsg = False
 
     @classmethod
     def from_shot(cls, shot: Shot, cam: Camera, param_rpc: dict, unit_data: dict) -> None:
@@ -110,10 +110,9 @@ def normalize_data(self, grid_img: np.ndarray, grid_world: np.ndarray) -> tuple:
         """
         if self.fact_rpc is None:
             if self.output_epsg:
-                ProjEngine().set_epsg_tf_geog_output(self.output_epsg)
-                x_geog, y_geog, z_geog = ProjEngine().tf.carto_to_geog_out(grid_world[0],
-                                                                           grid_world[1],
-                                                                           grid_world[2])
+                x_geog, y_geog, z_geog = ProjEngine().tf.proj_to_proj_out(grid_world[0],
+                                                                          grid_world[1],
+                                                                          grid_world[2])
             else:
                 x_geog, y_geog, z_geog = ProjEngine().tf.carto_to_geog(grid_world[0],
                                                                        grid_world[1],

diff --git a/borea/geodesy/proj_engine.py b/borea/geodesy/proj_engine.py
@@ -23,30 +23,23 @@ def __post_init__(self) -> None:
         if self.epsg:
             self.crs = pyproj.CRS.from_epsg(self.epsg[0])
             self.proj = pyproj.Proj(self.crs)
-            self.tf = TransformGeodesy(self.epsg, self.geoid)
+            self.tf = TransformGeodesy(self.epsg, self.geoid, self.epsg_output)
 
-    def set_epsg(self, epsg: list, geoid: list = None) -> None:
+    def set_epsg(self, epsg: list, geoid: list = None, epsg_output: int = None) -> None:
         """
         Setter of the class ProjEngine.
         Allows to init the class with data.
 
         Args:
-            epsg (list): Code epsg of the projection ex: 2154.
+            epsg (list): Code epsg of the projection ex: [2154].
             geoid (list): List of geoid to use.
+            epsg_output (int): Code epsg of the output projection. If you want to change.
         """
         self.epsg = epsg
         self.geoid = geoid
+        self.epsg_output = epsg_output
         self.__post_init__()
 
-    def set_epsg_tf_geog_output(self, epsg_output: int) -> None:
-        """
-        Create the pyproj Transformer from crs of worksite to crs geographic ask.
-
-        Args:
-            epsg_out (int): Code epsg of the output crs.
-        """
-        self.tf.tf_output(epsg_output)
-
     def get_meridian_convergence(self, x_carto: Union[np.ndarray, List[float], float],
                                  y_carto: Union[np.ndarray, List[float], float]) -> np.ndarray:
         """

diff --git a/borea/geodesy/transform_geodesy.py b/borea/geodesy/transform_geodesy.py
@@ -3,6 +3,7 @@
 """
 import pyproj
 import numpy as np
+import pandas as pd
 
 
 # pylint: disable=too-many-instance-attributes
@@ -14,10 +15,12 @@ class TransformGeodesy():
 
     Args:
         geoid (list): List of geoid to use.
-        crs (pyproj): CRS pyproj of the worksite.
+        epsg (list): Code epsg of the projection ex: [2154].
+        epsg_output (int): Code epsg of the output projection. If you want to change.
     """
-    def __init__(self, epsg: list, geoid: list) -> None:
+    def __init__(self, epsg: list, geoid: list, epsg_output: int) -> None:
         self.epsg = epsg
+        self.epsg_output = epsg_output
         self.crs = pyproj.CRS.from_epsg(epsg[0])
         self.geoid = geoid
         self._carto_to_geoc = None
@@ -26,7 +29,7 @@ def __init__(self, epsg: list, geoid: list) -> None:
         self._geog_to_carto = None
         self._geog_to_geoid = None
         self._geoid_to_geog = None
-        self.carto_to_geog_out = None
+        self._proj_to_proj_out = None
 
     @property
     def carto_to_geog(self) -> pyproj.Transformer:
@@ -150,6 +153,20 @@ def geog_to_geoid(self) -> pyproj.Transformer:
 
         return self._geog_to_geoid
 
+    @property
+    def proj_to_proj_out(self) -> pyproj.Transformer:
+        """
+        Create the pyproj Transformer from crs of worksite to crs geographic ask.
+
+        Returns:
+            pyproj.Transformer : carto_to_geog_out
+        """
+        if not self._proj_to_proj_out:
+            crs_out = pyproj.CRS.from_epsg(self.epsg_output)
+            self._proj_to_proj_out = pyproj.Transformer.from_crs(self.crs, crs_out).transform
+
+        return self._proj_to_proj_out
+
     def tranform_height(self, coor: np.ndarray) -> float:
         """
         Converting z in altitude to z in height of point.
@@ -172,7 +189,7 @@ def tranform_height(self, coor: np.ndarray) -> float:
             raise ValueError("The geoid has not been entered, "
                              "cannot transform z altitude to height.") from info
 
-        if new_z == np.inf:
+        if np.all(new_z == np.inf):
             raise ValueError("out geoid")
         return new_z
 
@@ -202,14 +219,22 @@ def tranform_altitude(self, coor: np.ndarray) -> float:
             raise ValueError("out geoid")
         return new_z
 
-    def tf_output(self, epsg_out: int = None) -> None:
+    def transform_pt_proj(self, df_pt: pd.DataFrame, type_z_input: str = None,
+                          type_z_output: str = None) -> pd.DataFrame:
         """
-        Create the pyproj Transformer from crs of worksite to crs geographic ask.
-
-        Args:
-            epsg_out (int): Code epsg of the output crs.
+        Tranform the input projection to the output projection of points coordinates
         """
-        if epsg_out and epsg_out != self.epsg[0]:
-            crs_out = pyproj.CRS.from_epsg(epsg_out)
-            crs_geog_out = pyproj.crs.GeographicCRS(name=crs_out.name, datum=crs_out.datum.name)
-            self.carto_to_geog_out = pyproj.Transformer.from_crs(self.crs, crs_geog_out).transform
+        if type_z_input and type_z_output:
+            if type_z_input != type_z_output:
+                if type_z_output == "altitude":
+                    df_pt["z"] = self.tranform_altitude(np.array([df_pt['x'],
+                                                                  df_pt['y'],
+                                                                  df_pt['z']]))
+                if type_z_output == "height":
+                    df_pt["z"] = self.tranform_height(np.array([df_pt['x'],
+                                                                df_pt['y'],
+                                                                df_pt['z']]))
+
+        df_pt["x"], df_pt["y"] = self.proj_to_proj_out(df_pt['x'], df_pt['y'])
+
+        return df_pt