RPCoptim: direct_loc_dtm with ESPG!=4326

refactor: WIP use direct_loc_dtm with dtm espg!=4326. RAF RPCoptim, if dtm epsg!=4326 use python los_extrema not the c++ one

refactor: RPCoptim: direct_loc_dtm with ESPG!=4326

CI:desactivate duplicated code on RPCoptim

refactor: few correction

refactor: rm useless attribut in DTMintersection class

CI

Author: August12vaux

shareloc/bindings/bind.cpp

@@ -42,7 +42,7 @@ PYBIND11_MODULE(bindings_cpp, m) {
         .def("eq_plan", &DTMIntersection::eq_plan)
         .def("ter_to_index", &DTMIntersection::ter_to_index)
         .def("index_to_ter", &DTMIntersection::index_to_ter)
-        .def("get_alt_offset", &DTMIntersection::get_alt_offset)
+        .def("get_footprint_corners", &DTMIntersection::get_footprint_corners)
         .def("interpolate", &DTMIntersection::interpolate)
         .def("intersect_dtm_cube", &DTMIntersection::intersect_dtm_cube)
         .def("intersection", &DTMIntersection::intersection)
@@ -124,7 +124,6 @@ PYBIND11_MODULE(bindings_cpp, m) {
                 p.set_transform(t[13].cast<std::array<double,6>>());
                 p.set_nb_rows(t[14].cast<int>());
                 p.set_nb_columns(t[15].cast<int>());
-
                 return p;
                 }
         ));
@@ -234,6 +233,7 @@ PYBIND11_MODULE(bindings_cpp, m) {
         .def("get_den_lon", &RPC::get_den_lon)
         .def("get_num_lat", &RPC::get_num_lat)
         .def("get_den_lat", &RPC::get_den_lat)
+        .def("get_alt_minmax", &RPC::get_alt_minmax)
         .def("get_offset_row", &RPC::get_offset_row)
         .def("get_scale_row", &RPC::get_scale_row)
         .def("get_offset_col", &RPC::get_offset_col)

shareloc/bindings/dtm_intersection.cpp

@@ -131,12 +131,14 @@ array<double, 3> DTMIntersection::index_to_ter(array<double, 3> const& vect_ter)
     return ter;
 }
 
-tuple<double,double> DTMIntersection::get_alt_offset(int epsg)const{
+py::array_t<double> DTMIntersection::get_footprint_corners()const{
 
-    if(epsg!=m_epsg){
-        cout<<"C++ : DTMintersection get_alt_offset not valid epsg"<<endl;
-    }
-    return{0.,0.};
+    vector<double> res = {-0.5,-0.5,\
+                        -0.5,this->get_nb_columns()-0.5,\
+                        this->get_nb_rows()-0.5,this->get_nb_columns()-0.5,\
+                        this->get_nb_rows()-0.5,-0.5};
+    
+    return py::array_t<double>({4, 2}, res.data());
 }
 
 double DTMIntersection::interpolate(double delta_shift_row, double delta_shift_col)const{

shareloc/bindings/dtm_intersection.hpp

@@ -61,8 +61,8 @@ class DTMIntersection
     /**index_to_ter*/
     std::array<double, 3> index_to_ter(std::array<double, 3> const& vect_ter)const;
 
-    /**get_alt_offset*/
-    std::tuple<double,double> get_alt_offset(int epsg)const;
+    /**get_footprint_corners*/
+    pybind11::array_t<double> get_footprint_corners()const;
 
     /**interpolate*/
     double interpolate(double delta_shift_row, double delta_shift_col) const;
@@ -164,10 +164,6 @@ class DTMIntersection
     /**set_nb_columns*/
     void set_nb_columns(int a) noexcept {m_nb_columns = a;};
 
-
-
-
-
 private:
 
         /**alt_data attribut*/

shareloc/bindings/rpc.cpp

@@ -161,6 +161,7 @@ tuple<double,double,double> RPC::direct_loc_dtm(
 {
     if(dtm.get_epsg() != 4326){
         throw runtime_error("C++ : direct_loc_dtm : epsg!=4326 -> Exiting");
+        //cout<<"C++ : direct_loc_dtm : epsg!=4326 -> get_diff_alti_min/max initialized ?"<<endl;
     }
 
     double min_dtm = dtm.get_alt_min() - 1.0;
@@ -204,11 +205,12 @@ tuple<vector<double>,vector<double>,vector<double>> RPC::direct_loc_dtm(
 {
     if(dtm.get_epsg() != 4326){
         throw runtime_error("C++ : direct_loc_dtm : epsg!=4326 -> Exiting");
+        //cout<<"C++ : direct_loc_dtm : epsg!=4326 -> get_diff_alti_min/max initialized ?"<<endl;
     }
 
     double min_dtm = dtm.get_alt_min() - 1.0;
     double max_dtm = dtm.get_alt_max() + 1.0;
-    
+
     size_t nb_points = row.size();
 
     vector<double> lon (2);

shareloc/bindings/rpc.hpp

@@ -207,6 +207,9 @@ class RPC : public GeoModelTemplate
     /**get_den_lat*/
     std::array<double, 20> const& get_den_lat() const noexcept {return m_den_lat;};
 
+    /**get_m_alt_minmax*/
+    std::array<double, 2> const& get_alt_minmax() const noexcept {return m_alt_minmax;};
+
     /**get_offset_row*/
     double get_offset_row() const noexcept {return m_offset_row;};
     /**get_scale_row*/

shareloc/geofunctions/dtm_intersection.py

@@ -31,11 +31,7 @@
 
 # Shareloc imports
 from shareloc.math_utils import interpol_bilin
-from shareloc.proj_utils import (
-    coordinates_conversion,
-    transform_index_to_physical_point,
-    transform_physical_point_to_index,
-)
+from shareloc.proj_utils import transform_index_to_physical_point, transform_physical_point_to_index
 
 
 class DTMIntersection:
@@ -179,27 +175,19 @@ def index_to_ter(self, vect_dtm):
         (vect_ter[1], vect_ter[0]) = transform_index_to_physical_point(self.transform, vect_dtm[0], vect_dtm[1])
         return vect_ter
 
-    def get_alt_offset(self, epsg):
+    def get_footprint_corners(self):
         """
-        returns min/amx altitude offset between dtm coordinates system and another one
+        get_footprint_corners method
 
-        :param epsg: epsg code to compare with
-        :type epsg: int
-        :return: min/max altimetric difference between epsg in parameter minus dtm alti expressed in dtm epsg
-        :rtype: list of float (1x2)
+        :return: footrpint corners of dtm
+        :rtype: np.ndarray(4x2) top_left, top_right, bottom_right, bottom_left
         """
-        if epsg != self.epsg:
-            alti_moy = (self.alt_min + self.alt_max) / 2.0
-            corners = np.zeros([4, 2])
-            corners[:, 0] = [0.0, 0.0, self.nb_rows, self.nb_rows]
-            corners[:, 1] = [0.0, self.nb_columns, self.nb_columns, 0.0]
-            corners -= 0.5  # index to corner
-            ground_corners = np.zeros([3, 4])
-            ground_corners[2, :] = alti_moy
-            ground_corners[1::-1, :] = transform_index_to_physical_point(self.transform, corners[:, 0], corners[:, 1])
-            converted_corners = coordinates_conversion(ground_corners.transpose(), self.epsg, epsg)
-            return [np.min(converted_corners[:, 2]) - alti_moy, np.max(converted_corners[:, 2]) - alti_moy]
-        return [0.0, 0.0]
+
+        corners = np.zeros([4, 2])
+        corners[:, 0] = [0.0, 0.0, self.nb_rows, self.nb_rows]
+        corners[:, 1] = [0.0, self.nb_columns, self.nb_columns, 0.0]
+        corners -= 0.5  # index to corner
+        return corners
 
     def interpolate(self, pos_row, pos_col):
         """
@@ -892,3 +880,6 @@ def get_alt_max(self):  # same api as cpp for direct_loc_dtm
 
     def get_epsg(self):  # same api as cpp for direct_loc_dtm
         return self.epsg
+
+    def get_transform(self):  # same api as cpp for direct_loc_dtm
+        return self.transform

shareloc/geomodels/rpc.py

@@ -27,16 +27,21 @@
 import logging
 import os
 from ast import literal_eval
+from typing import Union
 
 # Third party imports
 import numpy as np
+from affine import Affine
 from numba import config, njit, prange
 
+import bindings_cpp
+from shareloc.geofunctions.dtm_intersection import DTMIntersection
+
 # Shareloc imports
 from shareloc.geomodels.geomodel import GeoModel
 from shareloc.geomodels.geomodel_template import GeoModelTemplate
 from shareloc.geomodels.rpc_readers import rpc_reader
-from shareloc.proj_utils import coordinates_conversion
+from shareloc.proj_utils import coordinates_conversion, transform_index_to_physical_point
 
 # Set numba type of threading layer before parallel target compilation
 config.THREADING_LAYER = "omp"
@@ -315,21 +320,22 @@ def direct_loc_dtm(self, row, col, dtm):
             row = np.array([row])
             col = np.array([col])
 
-        diff_alti_min, diff_alti_max = dtm.get_alt_offset(self.epsg)
+        diff_alti_min, diff_alti_max = self.get_dtm_alt_offset(dtm.get_footprint_corners(), dtm)
+
         # print("min {} max {}".format(dtm.Zmin,dtm.Zmax))
         (min_dtm, max_dtm) = (dtm.get_alt_min() - 1.0 + diff_alti_min, dtm.get_alt_max() + 1.0 + diff_alti_max)
         if min_dtm < self.offset_alt - self.scale_alt:
             logging.debug("minimum dtm value is outside RPC validity domain, extrapolation will be done")
         if max_dtm > self.offset_alt + self.scale_alt:
             logging.debug("maximum dtm value is outside RPC validity domain, extrapolation will be done")
+
         los = self.los_extrema(row, col, min_dtm, max_dtm, epsg=dtm.get_epsg())
 
         # los -> (nb_point,nb_alt,3)
         los = los.T
         los = np.expand_dims(los, axis=0)
         los = np.moveaxis(los, 1, -1)
         los = los.reshape((len(col), 2, 3))
-
         direct_dtm = dtm.intersection_n_los_dtm(los)
 
         return direct_dtm
@@ -559,6 +565,29 @@ def get_alt_min_max(self):
         """
         return [self.offset_alt - self.scale_alt / 2.0, self.offset_alt + self.scale_alt / 2.0]
 
+    def get_dtm_alt_offset(self, corners: np.ndarray, dtm: Union[DTMIntersection, bindings_cpp.DTMIntersection]):
+        """
+        returns min/max altitude offset between dtm coordinates system and RPC one
+
+        :param corners: corners of the DTM's footprint
+        :type corners: np.ndarray (4x2)
+        :param dtm: DTM to get alt offset from
+        :type dtm: DTMIntersection or bindings_cpp.DTMIntersection
+        :return: min/max altimetric difference between RPC's epsg minus dtm alti expressed in dtm epsg
+        :rtype: list of float (1x2)
+        """
+
+        alti_moy = (dtm.get_alt_min() + dtm.get_alt_max()) / 2.0
+
+        ground_corners = np.zeros([3, 4])
+        ground_corners[2, :] = alti_moy
+        transform = dtm.get_transform()
+        if not isinstance(transform, Affine):
+            transform = Affine.from_gdal(*transform)
+        ground_corners[1::-1, :] = transform_index_to_physical_point(transform, corners[:, 0], corners[:, 1])
+        converted_corners = coordinates_conversion(ground_corners.transpose(), dtm.get_epsg(), self.epsg)
+        return [np.min(converted_corners[:, 2]) - alti_moy, np.max(converted_corners[:, 2]) - alti_moy]
+
     def los_extrema(self, row, col, alt_min=None, alt_max=None, fill_nan=False, epsg=None):
         """
         compute los extrema
@@ -603,6 +632,7 @@ def los_extrema(self, row, col, alt_min=None, alt_max=None, fill_nan=False, epsg
             col_array = np.array([col, col])
             alt_array = np.array([los_alt_max, los_alt_min])
         los_edges = self.direct_loc_h(row_array, col_array, alt_array, fill_nan)
+
         if extrapolate:
             diff = los_edges[0::2, :] - los_edges[1::2, :]
             delta_alt = diff[:, 2]
@@ -612,8 +642,10 @@ def los_extrema(self, row, col, alt_min=None, alt_max=None, fill_nan=False, epsg
             coeff_alt_min = np.tile(coeff_alt_min[:, np.newaxis], (1, 3))
             los_edges[0::2, :] = los_edges[1::2, :] + diff * coeff_alt_max
             los_edges[1::2, :] = los_edges[1::2, :] + diff * coeff_alt_min
+
         if epsg is not None and epsg != self.epsg:
             los_edges = coordinates_conversion(los_edges, self.epsg, epsg)
+
         return los_edges