Merge branch '215-n-utiliser-qu-une-seule-unite' into 'master'

Resolve "N'utiliser qu'une seule unité"

Closes #215

See merge request 3d/shareloc!137

shareloc/geofunctions/dtm_intersection.py

@@ -310,8 +310,8 @@ def intersect_dtm_cube(self, los):  # noqa: C901
         :param los:  line of sight
         :type los: numpy.array
         :return: intersection information
-            (True,an intersection has been found ?, (lon,lat) of dtm position, altitude)
-        :rtype: tuple (bool, bool, numpy.array, float)
+            (True,an intersection has been found ?, (lon,lat) of dtm position, altitude, line of sight in index frame)
+        :rtype: tuple (bool, bool, numpy.array, float, numpy.array)
         """
         # los: (n,3):
         point_b = None
@@ -455,54 +455,50 @@ def intersect_dtm_cube(self, los):  # noqa: C901
         # No solution if 0 or 1 single point is found (we have tangent to the cube)
         if nbi < 2:
             b_trouve = False
-            return True, b_trouve, point_b, h_intersect
+            return (True, b_trouve, point_b, h_intersect, los_index)
         # -----------------------------------------------------------------------
         # There are only 2 points left so we cross the cube
         # LAIG-FA-MAJA-2168-CNES: no more filtering on identical points. There may be a number of points > 2
         # Init the current point
         # DTM coordinates
-        point_dtm = np.zeros(3)
-        point_dtm[0] = coord_col_i[0]
-        point_dtm[1] = coord_row_i[0]
-        point_dtm[2] = coord_alt_i[0]
+        point_b = np.zeros(3)
+        point_b[0] = coord_col_i[0]
+        point_b[1] = coord_row_i[0]
+        point_b[2] = coord_alt_i[0]
         # point_dtm is the first intersection with the cube (line, column)
         # -----------------------------------------------------------------------
         # h is gld 3D
         h_intersect = alti_layer_i[0]
         # h_intersect is the h interpolation index (not integer)
         # -----------------------------------------------------------------------
-        # Terrain coordinates
-        point_b = self.index_to_ter(point_dtm)
-        # point_b is the terrain point (lon,lat)
-        # -----------------------------------------------------------------------
         # End, return
         b_trouve = True
-        return (True, b_trouve, point_b, h_intersect)
+
+        return (True, b_trouve, point_b, h_intersect, los_index)
 
     # gitlab issue #56
     # pylint: disable=too-many-locals
     # pylint: disable=too-many-function-args
     # pylint: disable=too-many-nested-blocks
     # pylint: disable=too-many-statements
-    def intersection(self, los, point_b, h_intersect):  # noqa: C901
+    def intersection(self, los_index, point_b, h_intersect):  # noqa: C901
         """
         DTM intersection
 
-        :param los:  line of sight
-        :type los: numpy.array
-        :param point_b:  position of intersection in DTM cube
+        :param los_index:  line of sight in index frame
+        :type los_index: numpy.array
+        :param point_b:  position of intersection in DTM cube in index frame
         :type point_b: numpy.array
         :param h_intersect:  altitude in DTM cube
         :type h_intersect: float
         :return: intersection information (True,an intersection has been found ?, position of intersection)
         :rtype: tuple (bool, bool, numpy.array)
         """
-        los_index = self.ters_to_indexs(los)
-        point_b_dtm = self.ter_to_index(point_b)
+
+        npl = los_index.shape[0]
         point_r = np.zeros(3)
-        (npl, _) = los.shape
         alti = np.arange(npl, -1.0, -1.0)
-        p_1 = point_b_dtm.copy()  # [p_1[0],p_1[1],p_1[2]]
+        p_1 = point_b.copy()  # [p_1[0],p_1[1],p_1[2]]
 
         h_intersect_p1 = h_intersect
 
@@ -528,7 +524,7 @@ def intersection(self, los, point_b, h_intersect):  # noqa: C901
         i_0 = int(np.floor(h_intersect_p1))
 
         #   1.3 - Init the starting point (in p_2)
-        p_2 = point_b_dtm.copy()
+        p_2 = point_b.copy()
         h_intersect_p2 = h_intersect
 
         nb_planes = los_index.shape[0]

shareloc/geomodels/grid.py

@@ -263,9 +263,9 @@ def direct_loc_dtm(self, row, col, dtm):
             row_i = row[point_index]
             col_i = col[point_index]
             los = self.compute_los(row_i, col_i, dtm.epsg)
-            (__, __, position_cube, alti) = dtm.intersect_dtm_cube(los)
+            (__, __, position_cube, alti, los_index) = dtm.intersect_dtm_cube(los)
             if position_cube is not None:
-                (__, __, points_dtm[point_index, :]) = dtm.intersection(los, position_cube, alti)
+                (__, __, points_dtm[point_index, :]) = dtm.intersection(los_index, position_cube, alti)
             else:
                 logging.warning("LOS doesn't instersect DTM cube")
                 points_dtm[point_index, :] = np.full(3, fill_value=np.nan)
@@ -371,9 +371,9 @@ def direct_loc_grid_dtm(self, row0, col0, steprow, stepcol, nbrow, nbcol, dtm):
                 col = col0 + stepcol * j
                 row = row0 + steprow * i
                 los = self.compute_los(row, col, dtm.epsg)
-                (__, __, position_cube, alti) = dtm.intersect_dtm_cube(los)
+                (__, __, position_cube, alti, los_index) = dtm.intersect_dtm_cube(los)
                 if position_cube is not None:
-                    (__, __, point_dtm) = dtm.intersection(los, position_cube, alti)
+                    (__, __, point_dtm) = dtm.intersection(los_index, position_cube, alti)
                 else:
                     point_dtm = np.full(3, fill_value=np.nan)
                 # conversion of all tab

shareloc/geomodels/rpc.py

@@ -695,9 +695,9 @@ def direct_loc_dtm(self, row, col, dtm):
         los = self.los_extrema(row, col, min_dtm, max_dtm, epsg=dtm.epsg)
         for i in range(points_nb):
             los_i = los[2 * i : 2 * i + 2, :]
-            (__, __, position_cube, alti) = dtm.intersect_dtm_cube(los_i)
+            (__, __, position_cube, alti, los_index) = dtm.intersect_dtm_cube(los_i)
             if position_cube is not None:
-                (__, __, position) = dtm.intersection(los_i, position_cube, alti)
+                (__, __, position) = dtm.intersection(los_index, position_cube, alti)
                 direct_dtm[i, :] = position
             else:
                 position = np.full(3, fill_value=np.nan)

tests/geofunctions/test_localization.py

@@ -170,13 +170,14 @@ def test_loc_dir_instersect_cube_dtm(col, row, valid_lon, valid_lat, valid_alt):
     """
     dtmbsq, gri = prepare_loc()
 
-    visee = np.zeros((3, gri.nbalt))
+    los = np.zeros((3, gri.nbalt))
     vislonlat = gri.interpolate_grid_in_plani(row, col)
-    visee[0, :] = vislonlat[0]
-    visee[1, :] = vislonlat[1]
-    visee[2, :] = gri.alts_down
-    visee = visee.T
-    (__, __, position, __) = dtmbsq.intersect_dtm_cube(visee)
+    los[0, :] = vislonlat[0]
+    los[1, :] = vislonlat[1]
+    los[2, :] = gri.alts_down
+    los = los.T
+    (__, __, position, __, __) = dtmbsq.intersect_dtm_cube(los)
+    position = dtmbsq.index_to_ter(position)
     assert position[0] == pytest.approx(valid_lon, abs=1e-12)
     assert position[1] == pytest.approx(valid_lat, abs=1e-12)
     assert position[2] == pytest.approx(valid_alt, abs=1e-12)
@@ -190,9 +191,9 @@ def test_loc_dir_interp_visee_unitaire_gld(row, col, valid_lon, valid_lat):
     Test los interpolation
     """
     ___, gri = prepare_loc()
-    visee = gri.interpolate_grid_in_plani(row, col)
-    assert visee[0][1] == pytest.approx(valid_lon, abs=1e-12)
-    assert visee[1][1] == pytest.approx(valid_lat, abs=1e-12)
+    los = gri.interpolate_grid_in_plani(row, col)
+    assert los[0][1] == pytest.approx(valid_lon, abs=1e-12)
+    assert los[1][1] == pytest.approx(valid_lat, abs=1e-12)
 
 
 @pytest.mark.parametrize("col,row,h", [(50.5, 100.5, 100.0)])
@@ -457,14 +458,14 @@ def test_loc_intersection(row, col, valid_lon, valid_lat, valid_alt):
     """
     dtmbsq, gri = prepare_loc()
 
-    visee = np.zeros((3, gri.nbalt))
+    los = np.zeros((3, gri.nbalt))
     vislonlat = gri.interpolate_grid_in_plani(row, col)
-    visee[0, :] = vislonlat[0]
-    visee[1, :] = vislonlat[1]
-    visee[2, :] = gri.alts_down
-    visee = visee.T
-    (__, __, point_b, alti) = dtmbsq.intersect_dtm_cube(visee)
-    (__, __, point_dtm) = dtmbsq.intersection(visee, point_b, alti)
+    los[0, :] = vislonlat[0]
+    los[1, :] = vislonlat[1]
+    los[2, :] = gri.alts_down
+    los = los.T
+    (__, __, point_b, alti, los_index) = dtmbsq.intersect_dtm_cube(los)
+    (__, __, point_dtm) = dtmbsq.intersection(los_index, point_b, alti)
     assert point_dtm[0] == pytest.approx(valid_lon, abs=1e-12)
     assert point_dtm[1] == pytest.approx(valid_lat, abs=1e-12)
     assert point_dtm[2] == pytest.approx(valid_alt, abs=1e-10)