Merge pull request #149 from IGNF/epsg_pyproj

Refactor ProjEngine and Transform class

.github/workflows/test_ubu.yml

@@ -35,6 +35,9 @@ jobs:
     - 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
 
+    - 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 0 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/
 

.github/workflows/test_win.yml

@@ -37,6 +37,10 @@ jobs:
       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
 
+    - 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 0 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}
       run: python ./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/

borea/datastruct/shot.py

@@ -153,9 +153,9 @@ def set_type_z(self, type_z: str) -> None:
             type_z (str): z type height or altitude.
         """
         if type_z == "height":
-            self.pos_shot[2] = ProjEngine().tranform_height(self.pos_shot)
+            self.pos_shot[2] = ProjEngine().tf.tranform_height(self.pos_shot)
         else:
-            self.pos_shot[2] = ProjEngine().tranform_altitude(self.pos_shot)
+            self.pos_shot[2] = ProjEngine().tf.tranform_altitude(self.pos_shot)
 
     def set_linear_alteration(self, linear_alteration: bool) -> None:
         """

borea/datastruct/workdata.py

@@ -59,20 +59,22 @@ def add_shot(self, name_shot: str, pos_shot: np.ndarray,
                                      linear_alteration=linear_alteration,
                                      order_axe=order_axe)
 
-    def set_proj(self, epsg: int, path_geoid: list = None) -> None:
+    def set_proj(self, epsg: list, path_geoid: list = None) -> None:
         """
         Setup a projection system to the worksite.
 
         Args:
-            epsg (int): Code epsg of the porjection ex: 2154.
+            epsg (list): Code epsg of the porjection ex: 2154.
             path_geoid (str): List of GeoTIFF which represents the geoid in grid form.
         """
         ProjEngine.clear()
         try:  # Check if the epsg exist
-            crs = CRS.from_epsg(epsg)
-            del crs
+            for idepsg in epsg:
+                if idepsg:
+                    _ = CRS.from_epsg(idepsg)
+                    del _
         except exceptions.CRSError as e_info:
-            raise exceptions.CRSError(f"Your EPSG:{epsg} doesn't exist") from e_info
+            raise exceptions.CRSError(f"Your EPSG:{epsg} doesn't exist in pyproj.") from e_info
 
         ProjEngine().set_epsg(epsg, path_geoid)
 

borea/format/rpc.py

@@ -111,13 +111,13 @@ 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().carto_to_geog_out(grid_world[0],
-                                                                        grid_world[1],
-                                                                        grid_world[2])
+                x_geog, y_geog, z_geog = ProjEngine().tf.carto_to_geog_out(grid_world[0],
+                                                                           grid_world[1],
+                                                                           grid_world[2])
             else:
-                x_geog, y_geog, z_geog = ProjEngine().carto_to_geog(grid_world[0],
-                                                                    grid_world[1],
-                                                                    grid_world[2])
+                x_geog, y_geog, z_geog = ProjEngine().tf.carto_to_geog(grid_world[0],
+                                                                       grid_world[1],
+                                                                       grid_world[2])
         else:
             x_geog = grid_world[0]*self.fact_rpc
             y_geog = grid_world[1]*self.fact_rpc

borea/geodesy/local_euclidean_proj.py

@@ -37,7 +37,7 @@ def mat_rot_euclidean_local(self, x: float, y: float) -> np.ndarray:
         Returns:
             np.array: Transition matrix.
         """
-        lon, lat = ProjEngine().carto_to_geog(x, y)
+        lon, lat = ProjEngine().tf.carto_to_geog(x, y)
         gamma = ProjEngine().get_meridian_convergence(x, y)
 
         # Matrix for switching to local cartesian coordinates
@@ -92,10 +92,10 @@ def world_to_eucli(self, coor: np.ndarray) -> np.ndarray:
         """
         coor = check_array_transfo(coor[0], coor[1], coor[2])
 
-        coor_geoc = np.array(ProjEngine().carto_to_geoc(coor[0], coor[1], coor[2]))
-        central_geoc = np.array(ProjEngine().carto_to_geoc(self.pt_central[0],
-                                                           self.pt_central[1],
-                                                           self.pt_central[2]))
+        coor_geoc = np.array(ProjEngine().tf.carto_to_geoc(coor[0], coor[1], coor[2]))
+        central_geoc = np.array(ProjEngine().tf.carto_to_geoc(self.pt_central[0],
+                                                              self.pt_central[1],
+                                                              self.pt_central[2]))
         dr = np.vstack([coor_geoc[0] - central_geoc[0],
                         coor_geoc[1] - central_geoc[1],
                         coor_geoc[2] - central_geoc[2]])
@@ -115,13 +115,13 @@ def eucli_to_world(self, coor: np.ndarray) -> np.ndarray:
         """
         coor = np.squeeze(coor)
 
-        central_geoc = np.array(ProjEngine().carto_to_geoc(self.pt_central[0],
-                                                           self.pt_central[1],
-                                                           self.pt_central[2]))
+        central_geoc = np.array(ProjEngine().tf.carto_to_geoc(self.pt_central[0],
+                                                              self.pt_central[1],
+                                                              self.pt_central[2]))
         dr = np.vstack([coor[0] - self.pt_central[0],
                         coor[1] - self.pt_central[1],
                         coor[2] - self.pt_central[2]])
         point_geoc = np.squeeze((self.rot_to_euclidean_local.T @ dr) + np.array([central_geoc]).T)
         x_gc, y_gc, z_gc = check_array_transfo(point_geoc[0], point_geoc[1], point_geoc[2])
-        tup = ProjEngine().geoc_to_carto(x_gc, y_gc, z_gc)
+        tup = ProjEngine().tf.geoc_to_carto(x_gc, y_gc, z_gc)
         return np.array([tup[0], tup[1], tup[2]])

borea/geodesy/proj_engine.py

@@ -11,7 +11,7 @@
 
 # pylint: disable=unpacking-non-sequence
 @dataclass
-class ProjEngine(TransformGeodesy, metaclass=Singleton):
+class ProjEngine(metaclass=Singleton):
     """
     This class provides functions for using a cartographic system.
     """
@@ -21,17 +21,17 @@ class ProjEngine(TransformGeodesy, metaclass=Singleton):
 
     def __post_init__(self) -> None:
         if self.epsg:
-            self.crs = pyproj.CRS.from_epsg(self.epsg)
+            self.crs = pyproj.CRS.from_epsg(self.epsg[0])
             self.proj = pyproj.Proj(self.crs)
-            TransformGeodesy.__tf_init__(self, self.geoid, self.crs)
+            self.tf = TransformGeodesy(self.epsg, self.geoid)
 
-    def set_epsg(self, epsg: int, geoid: list = None) -> None:
+    def set_epsg(self, epsg: list, geoid: list = None) -> None:
         """
         Setter of the class ProjEngine.
         Allows to init the class with data.
 
         Args:
-            epsg (int): Code epsg of the projection ex: 2154.
+            epsg (list): Code epsg of the projection ex: 2154.
             geoid (list): List of geoid to use.
         """
         self.epsg = epsg
@@ -45,7 +45,7 @@ def set_epsg_tf_geog_output(self, epsg_output: int) -> None:
         Args:
             epsg_out (int): Code epsg of the output crs.
         """
-        self.tf_output(self.crs, epsg_output)
+        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:
@@ -60,7 +60,7 @@ def get_meridian_convergence(self, x_carto: Union[np.ndarray, List[float], float
         Returns:
             np.array : Meridian convergence in degree.
         """
-        (x_geog, y_geog) = self.carto_to_geog(x_carto, y_carto)
+        (x_geog, y_geog) = self.tf.carto_to_geog(x_carto, y_carto)
         return -np.array(self.proj.get_factors(x_geog, y_geog).meridian_convergence)
 
     def get_scale_factor(self, x_carto: Union[np.ndarray, List[float], float],
@@ -76,5 +76,5 @@ def get_scale_factor(self, x_carto: Union[np.ndarray, List[float], float],
         Returns:
             np.array: Scale factor and meridian convergence.
         """
-        x_geog, y_geog = self.carto_to_geog(x_carto, y_carto)
+        x_geog, y_geog = self.tf.carto_to_geog(x_carto, y_carto)
         return np.array(self.proj.get_factors(x_geog, y_geog).meridional_scale) - 1

borea/geodesy/transform_geodesy.py

@@ -1,13 +1,11 @@
 """
 Module for class ProjEngine, transform geodesy
 """
-from dataclasses import dataclass
 import pyproj
 import numpy as np
 
 
-# pylint: disable=unsubscriptable-object
-@dataclass
+# pylint: disable=too-many-instance-attributes
 class TransformGeodesy():
     """
     This class provides functions to tranform coordinate system.
@@ -18,49 +16,139 @@ class TransformGeodesy():
         geoid (list): List of geoid to use.
         crs (pyproj): CRS pyproj of the worksite.
     """
-    carto_to_geoc = None
-    geoc_to_carto = None
-    carto_to_geog = None
-    geog_to_carto = None
-    geog_to_geoid = None
-    geoid_to_geog = None
-    carto_to_geog_out = None
-
-    def __tf_init__(self, geoid: list, crs: pyproj) -> None:
-        crs_geoc = pyproj.crs.GeocentricCRS(name=crs.name, datum=crs.datum.name)
-        crs_geog = pyproj.crs.GeographicCRS(name=crs.name, datum=crs.datum.name)
-        # Transform cartographic coordinates to geographic coordinates
-        self.carto_to_geog = pyproj.Transformer.from_crs(crs, crs_geog).transform
-        # Transform geographic coordinates to cartographic coordinates
-        self.geog_to_carto = pyproj.Transformer.from_crs(crs_geog, crs).transform
-        # Transform cartographic coordinates to geocentric coordinates
-        self.carto_to_geoc = pyproj.Transformer.from_crs(crs, crs_geoc).transform
-        # Transform geocentric coordinates to cartographic coordinates
-        self.geoc_to_carto = pyproj.Transformer.from_crs(crs_geoc, crs).transform
-
-        if geoid:
-            self.tf_geoid(geoid)
-
-    def tf_geoid(self, geoid: list) -> None:
-        """
-        Create attribute transform, to transform geographic coordinates to geoide coordinates.
+    def __init__(self, epsg: list, geoid: list) -> None:
+        self.epsg = epsg
+        self.crs = pyproj.CRS.from_epsg(epsg[0])
+        self.geoid = geoid
+        self._carto_to_geoc = None
+        self._geoc_to_carto = None
+        self._carto_to_geog = None
+        self._geog_to_carto = None
+        self._geog_to_geoid = None
+        self._geoid_to_geog = None
+        self.carto_to_geog_out = None
+
+    @property
+    def carto_to_geog(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
 
-        Args:
-            geoid (list): List of geoid to use.
+        Returns:
+            pyproj.Transformer : carto_to_geog
         """
-        try:
-            # Transform geoide coordinates to geographic coordinates
-            self.geoid_to_geog = pyproj.Transformer.from_pipeline(f"+proj=vgridshift "
-                                                                  f"+grids={','.join(geoid)} "
-                                                                  "+multiplier=1").transform
-            # Transform geographic coordinates to geoide coordinates
-            self.geog_to_geoid = pyproj.Transformer.from_pipeline(f"+proj=vgridshift "
-                                                                  f"+grids={','.join(geoid)} "
-                                                                  "+multiplier=-1").transform
-        except pyproj.exceptions.ProjError as e:
-            raise pyproj.exceptions.ProjError(f"{geoid} The name or path of geotif is incorrect or "
-                                              "does not exist in "
-                                              f"{pyproj.datadir.get_data_dir()}!!!{e}") from e
+        if not self._carto_to_geog:
+            try:
+                crs_geog = pyproj.CRS.from_epsg(self.epsg[1])
+            except (IndexError, pyproj.exceptions.CRSError):
+                crs_geog = pyproj.crs.GeographicCRS(name=self.crs.name, datum=self.crs.datum.name)
+
+            self._carto_to_geog = pyproj.Transformer.from_crs(self.crs, crs_geog).transform
+
+        return self._carto_to_geog
+
+    @property
+    def geog_to_carto(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
+
+        Returns:
+            pyproj.Transformer : geog_to_carto
+        """
+        if not self._geog_to_carto:
+            try:
+                crs_geog = pyproj.CRS.from_epsg(self.epsg[1])
+            except (IndexError, pyproj.exceptions.CRSError):
+                crs_geog = pyproj.crs.GeographicCRS(name=self.crs.name, datum=self.crs.datum.name)
+
+            self._geog_to_carto = pyproj.Transformer.from_crs(crs_geog, self.crs).transform
+
+        return self._geog_to_carto
+
+    @property
+    def carto_to_geoc(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
+
+        Returns:
+            pyproj.Transformer : carto_to_geoc
+        """
+        if not self._carto_to_geoc:
+            try:
+                crs_geoc = pyproj.CRS.from_epsg(self.epsg[2])
+            except (IndexError, pyproj.exceptions.CRSError):
+                crs_geoc = pyproj.crs.GeocentricCRS(name=self.crs.name, datum=self.crs.datum.name)
+
+            self._carto_to_geoc = pyproj.Transformer.from_crs(self.crs, crs_geoc).transform
+
+        return self._carto_to_geoc
+
+    @property
+    def geoc_to_carto(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
+
+        Returns:
+            pyproj.Transformer : geoc_to_carto
+        """
+        if not self._geoc_to_carto:
+            try:
+                crs_geoc = pyproj.CRS.from_epsg(self.epsg[2])
+            except (IndexError, pyproj.exceptions.CRSError):
+                crs_geoc = pyproj.crs.GeocentricCRS(name=self.crs.name, datum=self.crs.datum.name)
+
+            self._geoc_to_carto = pyproj.Transformer.from_crs(crs_geoc, self.crs).transform
+
+        return self._geoc_to_carto
+
+    @property
+    def geoid_to_geog(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
+
+        Returns:
+            pyproj.Transformer : geoid_to_geog
+        """
+        if not self.geoid:
+            raise ValueError("Mistake Geoid path")
+
+        if not self._geoid_to_geog:
+            try:
+                # Transform geoide coordinates to geographic coordinates
+                self._geoid_to_geog = pyproj.Transformer.from_pipeline("+proj=vgridshift "
+                                                                       "+grids="
+                                                                       f"{','.join(self.geoid)} "
+                                                                       "+multiplier=1").transform
+            except pyproj.exceptions.ProjError as e:
+                raise pyproj.exceptions.ProjError(f"{self.geoid} The name or path "
+                                                  "of geotif is incorrect or does not exist in "
+                                                  f"{pyproj.datadir.get_data_dir()}!!!{e}") from e
+
+        return self._geoid_to_geog
+
+    @property
+    def geog_to_geoid(self) -> pyproj.Transformer:
+        """
+        Returns the transformation or instantiates it before returning it.
+
+        Returns:
+            pyproj.Transformer : geog_to_geoid
+        """
+        if not self.geoid:
+            raise ValueError("Mistake Geoid path")
+
+        if not self._geog_to_geoid:
+            try:
+                # Transform geoide coordinates to geographic coordinates
+                self._geog_to_geoid = pyproj.Transformer.from_pipeline("+proj=vgridshift "
+                                                                       "+grids="
+                                                                       f"{','.join(self.geoid)} "
+                                                                       "+multiplier=-1").transform
+            except pyproj.exceptions.ProjError as e:
+                raise pyproj.exceptions.ProjError(f"{self.geoid} The name or path "
+                                                  "of geotif is incorrect or does not exist in "
+                                                  f"{pyproj.datadir.get_data_dir()}!!!{e}") from e
+
+        return self._geog_to_geoid
 
     def tranform_height(self, coor: np.ndarray) -> float:
         """
@@ -114,15 +202,14 @@ def tranform_altitude(self, coor: np.ndarray) -> float:
             raise ValueError("out geoid")
         return new_z
 
-    def tf_output(self, crs: pyproj, epsg_out: int = None) -> None:
+    def tf_output(self, epsg_out: int = None) -> None:
         """
         Create the pyproj Transformer from crs of worksite to crs geographic ask.
 
         Args:
-            crs (pyproj): The crs pyproj of the worksite.
             epsg_out (int): Code epsg of the output crs.
         """
-        if epsg_out and epsg_out != crs.to_epsg():
+        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(crs, crs_geog_out).transform
+            self.carto_to_geog_out = pyproj.Transformer.from_crs(self.crs, crs_geog_out).transform