Set ENU as the default coordinates for cubes LOS and along-track unit vectors (#772)

* add new parameter to set the EPSG for LOS and along-track unit vectors at C++ level

* updating python bindings after adding new parameter to set the EPSG for LOS and along-track unit vectors

* update python unitest

* update C++ unitest

* Update interface with functions that create metadata data cubes with the new parameter to set the EPSG for LOS and along-track unit vectors

* add message to unitest to inform which EPSG is being processed

* move new parameter to the end for backwards compatibility

* move definition of projection for LOS and along-track unit vectors outside the writeVectorDerivedCubes() function to improve efficiency

* set default epsg_los_and_along_track_vectors to zero in the header; update docstrings

Author: gshiroma

cxx/isce3/geometry/metadataCubes.cpp

@@ -64,7 +64,7 @@ static inline void writeVectorDerivedCubes(const int array_pos_i,
         const isce3::core::Vec3& target_llh,
         const isce3::core::Vec3& target_proj, const isce3::core::Orbit& orbit,
         const isce3::core::Ellipsoid& ellipsoid,
-        const isce3::core::ProjectionBase* proj,
+        const isce3::core::ProjectionBase* proj_los_and_along_track_vectors,
         isce3::io::Raster* incidence_angle_raster,
         isce3::core::Matrix<float>& incidence_angle_array,
         isce3::io::Raster* los_unit_vector_x_raster,
@@ -132,8 +132,8 @@ static inline void writeVectorDerivedCubes(const int array_pos_i,
     }
 
 
-    // If lat/lon, compute vectors in ENU coordinates around target
-    if (proj->code() == 4326) {
+    // If null, compute vectors in ENU coordinates around target
+    if (proj_los_and_along_track_vectors == nullptr) {
 
         // LOS unit vector X (ENU)
         if (los_unit_vector_x_raster != nullptr) {
@@ -173,7 +173,7 @@ static inline void writeVectorDerivedCubes(const int array_pos_i,
         const isce3::core::Vec3 target_to_sat_next_llh = 
             ellipsoid.xyzToLonLat(target_to_sat_next_xyz);
         isce3::core::Vec3 target_to_sat_next_proj = 
-            proj->forward(target_to_sat_next_llh);
+            proj_los_and_along_track_vectors->forward(target_to_sat_next_llh);
         const isce3::core::Vec3 look_vector_proj =
                 (target_to_sat_next_proj - target_proj).normalized();
 
@@ -233,11 +233,13 @@ static inline void writeVectorDerivedCubes(const int array_pos_i,
         // Compute velocity vector (proj)
         const isce3::core::Vec3 sat_next_llh =
                 ellipsoid.xyzToLonLat(sat_next_xyz);
-        const isce3::core::Vec3 sat_next_proj = proj->forward(sat_next_llh);
+        const isce3::core::Vec3 sat_next_proj = 
+            proj_los_and_along_track_vectors->forward(sat_next_llh);
 
         const isce3::core::Vec3 sat_previous_llh =
                 ellipsoid.xyzToLonLat(sat_xyz_previous);
-        const isce3::core::Vec3 sat_previous_proj = proj->forward(sat_previous_llh);
+        const isce3::core::Vec3 sat_previous_proj = 
+            proj_los_and_along_track_vectors->forward(sat_previous_llh);
 
         const isce3::core::Vec3 along_track_unit_vector =
                 (sat_next_proj - sat_previous_proj).normalized();
@@ -262,6 +264,7 @@ void makeRadarGridCubes(const isce3::product::RadarGridParameters& radar_grid,
                         const isce3::core::Orbit& orbit,
                         const isce3::core::LUT2d<double>& native_doppler,
                         const isce3::core::LUT2d<double>& grid_doppler,
+                        const int epsg_los_and_along_track_vectors,
                         isce3::io::Raster* slant_range_raster,
                         isce3::io::Raster* azimuth_time_raster,
                         isce3::io::Raster* incidence_angle_raster,
@@ -286,6 +289,12 @@ void makeRadarGridCubes(const isce3::product::RadarGridParameters& radar_grid,
     for (int height_count = 0; height_count < heights.size(); ++height_count) {
 
         auto proj = isce3::core::makeProjection(geogrid.epsg());
+        std::unique_ptr<ProjectionBase>
+            proj_los_and_along_track_vectors =
+                (epsg_los_and_along_track_vectors == 0 or
+                 epsg_los_and_along_track_vectors == 4326) ?
+                 nullptr :
+                 isce3::core::makeProjection(epsg_los_and_along_track_vectors);
 
         const isce3::core::Ellipsoid& ellipsoid = proj->ellipsoid();
 
@@ -377,10 +386,10 @@ void makeRadarGridCubes(const isce3::product::RadarGridParameters& radar_grid,
 
                 isce3::geometry::writeVectorDerivedCubes(
                         i, j, native_azimuth_time, target_llh, target_proj,
-                        orbit, ellipsoid, proj.get(), incidence_angle_raster,
-                        incidence_angle_array, los_unit_vector_x_raster,
-                        los_unit_vector_x_array, los_unit_vector_y_raster,
-                        los_unit_vector_y_array,
+                        orbit, ellipsoid, proj_los_and_along_track_vectors.get(), 
+                        incidence_angle_raster, incidence_angle_array, 
+                        los_unit_vector_x_raster, los_unit_vector_x_array, 
+                        los_unit_vector_y_raster, los_unit_vector_y_array,
                         along_track_unit_vector_x_raster,
                         along_track_unit_vector_x_array,
                         along_track_unit_vector_y_raster,
@@ -448,6 +457,7 @@ void makeGeolocationGridCubes(
         const std::vector<double>& heights, const isce3::core::Orbit& orbit,
         const isce3::core::LUT2d<double>& native_doppler,
         const isce3::core::LUT2d<double>& grid_doppler, const int epsg,
+        const int epsg_los_and_along_track_vectors,
         isce3::io::Raster* coordinate_x_raster,
         isce3::io::Raster* coordinate_y_raster,
         isce3::io::Raster* incidence_angle_raster,
@@ -471,9 +481,15 @@ void makeGeolocationGridCubes(
     for (int height_count = 0; height_count < heights.size(); ++height_count) {
 
         auto proj = isce3::core::makeProjection(epsg);
-                
         const isce3::core::Ellipsoid& ellipsoid = proj->ellipsoid();
 
+        std::unique_ptr<ProjectionBase>
+            proj_los_and_along_track_vectors =
+                (epsg_los_and_along_track_vectors == 0 or
+                 epsg_los_and_along_track_vectors == 4326) ?
+                 nullptr :
+                 isce3::core::makeProjection(epsg_los_and_along_track_vectors);
+
         auto coordinate_x_array = 
                 getNanArrayRadarGrid<double>(coordinate_x_raster, radar_grid);
         auto coordinate_y_array =
@@ -561,7 +577,7 @@ void makeGeolocationGridCubes(
 
                 writeVectorDerivedCubes(
                         i, j, native_azimuth_time, target_llh, target_proj,
-                        orbit, ellipsoid, proj.get(),
+                        orbit, ellipsoid, proj_los_and_along_track_vectors.get(),
                         incidence_angle_raster, incidence_angle_array,
                         los_unit_vector_x_raster, los_unit_vector_x_array,
                         los_unit_vector_y_raster, los_unit_vector_y_array,

cxx/isce3/geometry/metadataCubes.h

@@ -19,7 +19,8 @@ namespace isce3 { namespace geometry {
  * projection system
  * @param[in]  orbit                            Reference orbit
  * @param[in]  ellipsoid                        Reference ellipsoid
- * @param[in]  proj                             Output projection
+ * @param[in]  proj_los_and_along_track_vectors ProjectionBase object for 
+ * representing LOS and along-track unit vectors
  * @param[out] incidence_angle_raster           Incidence angle cube raster
  * @param[out] incidence_angle_array            Incidence angle cube array
  * @param[out] los_unit_vector_x_raster         LOS (target-to-sensor) unit 
@@ -42,7 +43,7 @@ static inline void writeVectorDerivedCubes(
         const double native_azimuth_time, const isce3::core::Vec3& target_llh,
         const isce3::core::Vec3& target_proj, const isce3::core::Orbit& orbit,
         const isce3::core::Ellipsoid& ellipsoid,
-        const isce3::core::ProjectionBase* proj,
+        const isce3::core::ProjectionBase* proj_los_and_along_track_vectors,
         isce3::io::Raster* incidence_angle_raster,
         isce3::core::Matrix<float>& incidence_angle_array,
         isce3::io::Raster* los_unit_vector_x_raster,
@@ -78,16 +79,19 @@ static inline void writeVectorDerivedCubes(
  * high-resolution maps of the corresponding radar geometry variable.
  * 
  * The line-of-sight (LOS) and along-track unit vectors are referenced to
- * the projection defined by the geogrid EPSG code. In case the EPSG code
- * is equal to 4326, ENU coordinates wrt targets are used instead.
+ * the projection defined by the epsg_los_and_along_track_vectors code. 
+ * In the case of ENU, i.e. epsg_los_and_along_track_vectors equals to 0 
+ * or 4326, ENU coordinates are computed wrt targets.
  *
  * @param[in]  radar_grid                  Radar grid
  * @param[in]  geogrid                     Geogrid to be
  *  used as reference for output cubes
  * @param[in]  heights                     Cube heights [m]
  * @param[in]  orbit                       Reference orbit
  * @param[in]  native_doppler              Native image Doppler
- * @param[in]  grid_doppler                Grid Doppler.
+ * @param[in]  grid_doppler                Grid Doppler
+ * @param[in]  epsg_los_and_along_track_vectors EPSG code for LOS and 
+ * along-track unit vectors (0 or 4326 for ENU coordinates)
  * @param[out] slant_range_raster          Slant-range (in meters) 
  * cube raster
  * @param[out] azimuth_time_raster         Azimuth time (in seconds) 
@@ -115,6 +119,7 @@ void makeRadarGridCubes(
         const std::vector<double>& heights, const isce3::core::Orbit& orbit,
         const isce3::core::LUT2d<double>& native_doppler,
         const isce3::core::LUT2d<double>& grid_doppler,
+        const int epsg_los_and_along_track_vectors = 0,
         isce3::io::Raster* slant_range_raster = nullptr,
         isce3::io::Raster* azimuth_time_raster = nullptr,
         isce3::io::Raster* incidence_angle_raster = nullptr,
@@ -147,15 +152,18 @@ void makeRadarGridCubes(
  * high-resolution maps of the corresponding radar geometry variable.
  * 
  * The line-of-sight (LOS) and along-track unit vectors are referenced to
- * the projection defined by the EPSG code. In case the EPSG code
- * is equal to 4326, ENU coordinates wrt targets are used instead.
+ * the projection defined by the epsg_los_and_along_track_vectors code. 
+ * In the case of ENU, i.e. epsg_los_and_along_track_vectors equals to 0 
+ * or 4326, ENU coordinates are computed wrt targets.
  *
  * @param[in]  radar_grid                Cube radar grid
  * @param[in]  heights                   Cube heights
  * @param[in]  orbit                     Orbit
  * @param[in]  native_doppler            Native image Doppler
  * @param[in]  grid_doppler              Grid Doppler.
  * @param[in]  epsg                      Output geolocation EPSG
+ * @param[in]  epsg_los_and_along_track_vectors EPSG code for LOS and 
+ * along-track unit vectors (0 or 4326 for ENU coordinates)
  * @param[out] coordinate_x_raster       Geolocation coordinate X raster
  * @param[out] coordinate_y_raster       Geolocation coordinage Y raster
  * @param[out] incidence_angle_raster    Incidence angle (in degrees wrt 
@@ -180,7 +188,7 @@ void makeGeolocationGridCubes(
         const std::vector<double>& heights, const isce3::core::Orbit& orbit,
         const isce3::core::LUT2d<double>& native_doppler,
         const isce3::core::LUT2d<double>& grid_doppler,
-        const int epsg,
+        const int epsg, const int epsg_los_and_along_track_vectors = 0,
         isce3::io::Raster* coordinate_x_raster = nullptr,
         isce3::io::Raster* coordinate_y_raster = nullptr,
         isce3::io::Raster* incidence_angle_raster = nullptr,

python/extensions/pybind_isce3/geometry/metadataCubes.cpp

@@ -14,6 +14,7 @@ void addbinding_metadata_cubes(py::module & m)
     m.def("make_radar_grid_cubes", &isce3::geometry::makeRadarGridCubes,
           py::arg("radar_grid"), py::arg("geogrid"), py::arg("heights"),
           py::arg("orbit"), py::arg("native_doppler"), py::arg("grid_doppler"),
+          py::arg("epsg_los_and_along_track_vectors") = 0,
           py::arg("slant_range_raster") = nullptr,
           py::arg("azimuth_time_raster") = nullptr,
           py::arg("incidence_angle_raster") = nullptr,
@@ -45,8 +46,9 @@ void addbinding_metadata_cubes(py::module & m)
              high-resolution maps of the corresponding radar geometry variable.
 
              The line-of-sight (LOS) and along-track unit vectors are referenced to
-             the projection defined by the geogrid EPSG code. In case the EPSG code
-             is equal to 4326, ENU coordinates wrt targets are used instead.
+             the projection defined by the epsg_los_and_along_track_vectors code. 
+             In the case of ENU, i.e. epsg_los_and_along_track_vectors equals to 0 
+             or 4326, ENU coordinates are computed wrt targets.
 
         Parameters
         ----------
@@ -62,6 +64,9 @@ void addbinding_metadata_cubes(py::module & m)
                   Native image Doppler
               grid_doppler : isce3.core.LUT2d
                   Grid Doppler
+              epsg_los_and_along_track_vectors : int, optional
+                  EPSG code for LOS and along-track unit vectors 
+                  (0 or 4326 for ENU coordinates)
               slant_range_raster : isce3.io.Raster, optional
                   Slant-range (in meters) cube raster
               azimuth_time_raster : isce3.io.Raster, optional
@@ -95,6 +100,7 @@ void addbinding_metadata_cubes(py::module & m)
         py::arg("native_doppler"),
         py::arg("grid_doppler"),
         py::arg("epsg"),
+        py::arg("epsg_los_and_along_track_vectors") = 0, 
         py::arg("coordinate_x_raster") = nullptr,
         py::arg("coordinate_y_raster") = nullptr,
         py::arg("incidence_angle_raster") = nullptr,
@@ -127,8 +133,9 @@ void addbinding_metadata_cubes(py::module & m)
            high-resolution maps of the corresponding radar geometry variable.
 
            The line-of-sight (LOS) and along-track unit vectors are referenced to
-           the projection defined by the EPSG code. In case the EPSG code
-           is equal to 4326, ENU coordinates wrt targets are used instead.
+           the projection defined by the epsg_los_and_along_track_vectors code. 
+           In the case of ENU, i.e. epsg_los_and_along_track_vectors equals to 0 
+           or 4326, ENU coordinates are computed wrt targets.
 
         Parameters
         ----------
@@ -144,6 +151,9 @@ void addbinding_metadata_cubes(py::module & m)
                   Grid Doppler
               epsg : int
                   Output geolocation EPSG
+              epsg_los_and_along_track_vectors : int, optional
+                  EPSG code for LOS and along-track unit vectors 
+                  (0 or 4326 for ENU coordinates)
               coordinate_x_raster : isce3.io.Raster, optional
                  Geolocation coordinate X raster
               coordinate_y_raster : isce3.io.Raster, optional

python/packages/pybind_nisar/workflows/h5_prep.py

@@ -1092,11 +1092,12 @@ def set_get_geo_info(hdf5_obj, root_ds, geo_grid, z_vect=None, flag_cube=False):
     return yds, xds
 
 
-def add_radar_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, geogrid, heights,
-                                 radar_grid,
-                                 orbit, native_doppler, grid_doppler,
+def add_radar_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, geogrid, 
+                                 heights, radar_grid, orbit, 
+                                 native_doppler, grid_doppler,
                                  threshold_geo2rdr=1e-8,
-                                 numiter_geo2rdr=100, delta_range=1e-8):
+                                 numiter_geo2rdr=100, delta_range=1e-8,
+                                 epsg_los_and_along_track_vectors=0):
     if cube_group_name not in hdf5_obj:
         cube_group = hdf5_obj.create_group(cube_group_name)
     else:
@@ -1105,8 +1106,7 @@ def add_radar_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, geogrid, heights,
     cube_shape = [len(heights), geogrid.length, geogrid.width]
 
     zds, yds, xds = set_get_geo_info(hdf5_obj, cube_group_name, geogrid,
-                                     z_vect=heights,
-                                     flag_cube=True)
+                                     z_vect=heights, flag_cube=True)
 
     # seconds since ref epoch
     ref_epoch = radar_grid.ref_epoch
@@ -1168,6 +1168,7 @@ def add_radar_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, geogrid, heights,
                                          orbit,
                                          native_doppler,
                                          grid_doppler,
+                                         epsg_los_and_along_track_vectors,
                                          slant_range_raster,
                                          azimuth_time_raster,
                                          incidence_angle_raster,
@@ -1231,10 +1232,11 @@ def _get_raster_from_hdf5_ds(group, ds_name, dtype, shape,
 
 
 def add_geolocation_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, radar_grid,
-                                       heights,
-                                       orbit, native_doppler, grid_doppler,
-                                       epsg, threshold_geo2rdr=1e-8,
-                                       numiter_geo2rdr=100, delta_range=1e-8):
+                                       heights, orbit, native_doppler, 
+                                       grid_doppler, epsg,
+                                       threshold_geo2rdr=1e-8,
+                                       numiter_geo2rdr=100, delta_range=1e-8, 
+                                       epsg_los_and_along_track_vectors=0):
     if cube_group_name not in hdf5_obj:
         cube_group = hdf5_obj.create_group(cube_group_name)
     else:
@@ -1307,6 +1309,7 @@ def add_geolocation_grid_cubes_to_hdf5(hdf5_obj, cube_group_name, radar_grid,
                                           native_doppler,
                                           grid_doppler,
                                           epsg,
+                                          epsg_los_and_along_track_vectors,
                                           coordinate_x_raster,
                                           coordinate_y_raster,
                                           incidence_angle_raster,

tests/cxx/isce3/geometry/metadata_cubes/metadata_cubes.cpp

@@ -454,9 +454,14 @@ TEST(radarGridCubeTest, testRadarGridCube)
         isce3::product::GeoGridParameters geogrid(x0, y0, dx, dy, width, length,
                                                   epsg);
 
+        const int epsg_los_and_along_track_vectors = epsg;
+
         // Make cubes
+        std::cout << "calling makeRadarGridCubes() with geogrid EPSG:" 
+                  << epsg << std::endl;
         isce3::geometry::makeRadarGridCubes(
-                radar_grid, geogrid, heights, orbit, zero_doppler, zero_doppler,
+                radar_grid, geogrid, heights, orbit, zero_doppler, 
+                zero_doppler, epsg_los_and_along_track_vectors,
                 &slant_range_raster, &azimuth_time_raster,
                 &incidence_angle_raster, &los_unit_vector_x_raster,
                 &los_unit_vector_y_raster, &along_track_unit_vector_x_raster,
@@ -698,6 +703,7 @@ TEST(metadataCubesTest, testMetadataCubes) {
     int numiter_geo2rdr = 25;
     double delta_range = 1e-6;
     int epsg = 4326;
+    int epsg_los_and_along_track_vectors = epsg;
 
     std::vector<double> heights = {0.0, 1000.0};
 
@@ -730,7 +736,8 @@ TEST(metadataCubesTest, testMetadataCubes) {
     // Make cubes
     isce3::geometry::makeGeolocationGridCubes(
             radar_grid, heights, orbit, zero_doppler, zero_doppler, epsg,
-            &coordinate_x_raster, &coordinate_y_raster, &incidence_angle_raster,
+            epsg_los_and_along_track_vectors, &coordinate_x_raster, 
+            &coordinate_y_raster, &incidence_angle_raster,
             &los_unit_vector_x_raster, &los_unit_vector_y_raster,
             &along_track_unit_vector_x_raster, &along_track_unit_vector_y_raster,
             &elevation_angle_raster, threshold_geo2rdr, numiter_geo2rdr,