Skip to content

Commit cd18cc2

Browse files
gshiromagmgunter
authored and
GitHub Enterprise
committed
Fix topo for DEMs in projected coordinates (#780)
* rename functions starting with GetDemCoords to getDemCoords * substitute topo interpolateXY() with getDemCoords() * update envisat golden dataset * update winnipeg golden dataset * make topo unitest print which files and bands it is checking * use ECEF coordinates rather than Earth radius to calculate height derivatives * remove unnecessary code * bring changes to CUDA topo * bring changes to CUDA topo (2) * bring changes to CUDA topo (3) * add z-coordinate to output projection inverse transform for cases in which the input ellipsoid differs from the output ellipsoid * add z-coordinate to output projection inverse transform for cases in which the input ellipsoid differs from the output ellipsoid (2) * declare variable and assign value in the same line * revert changes: declare variable and assign value in the same line * add docstrings to getDemCoordsSameEpsg() and getDemCoordsDiffEpsg() * Update cxx/isce3/geometry/Topo.cpp Co-authored-by: Geoffrey M Gunter <[email protected]> Co-authored-by: Geoffrey M Gunter <[email protected]>
1 parent a9cd337 commit cd18cc2

File tree

8 files changed

+110
-18
lines changed

8 files changed

+110
-18
lines changed

‎cxx/isce3/cuda/geometry/gpuTopo.cu‎

Lines changed: 32 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -114,16 +114,41 @@ void setOutputTopoLayers(const Vec3& targetLLH,
114114
}
115115
layers.hdg(index, heading);
116116

117+
// Project output coordinates to DEM coordinates
118+
Vec3 input_coords_llh;
119+
(*projOutput)->inverse({x, y, targetLLH[2]}, input_coords_llh);
120+
Vec3 dem_vect;
121+
(*(demInterp.proj()))->forward(input_coords_llh, dem_vect);
122+
117123
// East-west slope using central difference
118-
double aa = demInterp.interpolateXY(x - demInterp.deltaX(), y);
119-
double bb = demInterp.interpolateXY(x + demInterp.deltaX(), y);
120-
double gamma = targetLLH[1];
121-
double alpha = ((bb - aa) * degrees) / (2.0 * ellipsoid.rEast(gamma) * demInterp.deltaX());
124+
double aa = demInterp.interpolateXY(dem_vect[0] - demInterp.deltaX(), dem_vect[1]);
125+
double bb = demInterp.interpolateXY(dem_vect[0] + demInterp.deltaX(), dem_vect[1]);
126+
127+
Vec3 dem_vect_p_dx = {dem_vect[0] + demInterp.deltaX(), dem_vect[1], dem_vect[2]};
128+
Vec3 dem_vect_m_dx = {dem_vect[0] - demInterp.deltaX(), dem_vect[1], dem_vect[2]};
129+
Vec3 input_coords_llh_p_dx, input_coords_llh_m_dx;
130+
(*(demInterp.proj()))->inverse(dem_vect_p_dx, input_coords_llh_p_dx);
131+
(*(demInterp.proj()))->inverse(dem_vect_m_dx, input_coords_llh_m_dx);
132+
const Vec3 input_coords_xyz_p_dx = ellipsoid.lonLatToXyz(input_coords_llh_p_dx);
133+
const Vec3 input_coords_xyz_m_dx = ellipsoid.lonLatToXyz(input_coords_llh_m_dx);
134+
double dx = (input_coords_xyz_p_dx - input_coords_xyz_m_dx).norm();
135+
136+
double alpha = (bb - aa) / dx;
122137

123138
// North-south slope using central difference
124-
aa = demInterp.interpolateXY(x, y - demInterp.deltaY());
125-
bb = demInterp.interpolateXY(x, y + demInterp.deltaY());
126-
double beta = ((bb - aa) * degrees) / (2.0 * ellipsoid.rNorth(gamma) * demInterp.deltaY());
139+
aa = demInterp.interpolateXY(dem_vect[0], dem_vect[1] - demInterp.deltaY());
140+
bb = demInterp.interpolateXY(dem_vect[0], dem_vect[1] + demInterp.deltaY());
141+
142+
Vec3 dem_vect_p_dy = {dem_vect[0], dem_vect[1] + demInterp.deltaY(), dem_vect[2]};
143+
Vec3 dem_vect_m_dy = {dem_vect[0], dem_vect[1] - demInterp.deltaY(), dem_vect[2]};
144+
Vec3 input_coords_llh_p_dy, input_coords_llh_m_dy;
145+
(*(demInterp.proj()))->inverse(dem_vect_p_dy, input_coords_llh_p_dy);
146+
(*(demInterp.proj()))->inverse(dem_vect_m_dy, input_coords_llh_m_dy);
147+
const Vec3 input_coords_xyz_p_dy = ellipsoid.lonLatToXyz(input_coords_llh_p_dy);
148+
const Vec3 input_coords_xyz_m_dy = ellipsoid.lonLatToXyz(input_coords_llh_m_dy);
149+
double dy = (input_coords_xyz_p_dy - input_coords_xyz_m_dy).norm();
150+
151+
double beta = (bb - aa) / dy;
127152

128153
// Compute local incidence angle
129154
enu /= enu.norm();

‎cxx/isce3/geometry/RTC.h‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -329,10 +329,34 @@ void computeRtcAreaProj(isce3::io::Raster& dem,
329329
isce3::core::dataInterpMethod::BIQUINTIC_METHOD,
330330
double threshold = 1e-8, int num_iter = 100, double delta_range = 1e-8);
331331

332+
/*
333+
Interpolate DEM at position (x, y) considering that input_proj and
334+
dem_interp have same coordinate systems. The function is written to
335+
have the same interface of getDemCoordsDiffEpsg()
336+
* @param[in] x X-coordinate in input coordinates
337+
* @param[in] y Y-coordinate in input coordinates
338+
* @param[in] dem_interp DEM interpolation object
339+
* @param[in] input_proj Input projection object
340+
* @returns 3-elements vector containing the x and
341+
* y coordinates over DEM projection coordinates, and interpolated
342+
* DEM value at that position: {x_dem, y_dem, z_dem}
343+
*/
332344
inline isce3::core::Vec3 getDemCoordsSameEpsg(double x, double y,
333345
const DEMInterpolator& dem_interp,
334346
isce3::core::ProjectionBase* input_proj);
335347

348+
/*
349+
Convert x and y coordinates to from input_proj coordinates to
350+
DEM (dem_interp) coordinates and interpolate DEM at that position.
351+
3-elements vector containing the
352+
* @param[in] x X-coordinate in input coordinates
353+
* @param[in] y Y-coordinate in input coordinates
354+
* @param[in] dem_interp DEM interpolation object
355+
* @param[in] input_proj Input projection object
356+
* @returns 3-elements vector containing the x and
357+
* y coordinates over DEM projection coordinates, and interpolated
358+
* DEM value at that position: {x_dem, y_dem, z_dem}
359+
*/
336360
inline isce3::core::Vec3 getDemCoordsDiffEpsg(double x, double y,
337361
const DEMInterpolator& dem_interp,
338362
isce3::core::ProjectionBase* input_proj);

‎cxx/isce3/geometry/Topo.cpp‎

Lines changed: 47 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include <isce3/product/Product.h>
2121

2222
// isce3::geometry
23+
#include <isce3/geometry/RTC.h>
2324
#include "DEMInterpolator.h"
2425
#include "TopoLayers.h"
2526

@@ -143,6 +144,9 @@ topo(Raster & demRaster, TopoLayers & layers)
143144
const double endingRange = _radarGrid.endingRange();
144145
const double midRange = _radarGrid.midRange();
145146

147+
info << "DEM EPSG: " << demRaster.getEPSG() << pyre::journal::newline;
148+
info << "Output EPSG: " << _epsgOut << pyre::journal::endl;
149+
146150
// Loop over blocks
147151
size_t totalconv = 0;
148152
for (size_t block = 0; block < nBlocks; ++block) {
@@ -619,16 +623,39 @@ _setOutputTopoLayers(Vec3 & targetLLH, TopoLayers & layers, size_t line,
619623
}
620624
layers.hdg(line, bin, heading);
621625

626+
// Project output coordinates to DEM coordinates
627+
auto input_coords_llh = _proj->inverse({x, y, targetLLH[2]});
628+
Vec3 dem_vect = demInterp.proj()->forward(input_coords_llh);
629+
622630
// East-west slope using central difference
623-
double aa = demInterp.interpolateXY(x - demInterp.deltaX(), y);
624-
double bb = demInterp.interpolateXY(x + demInterp.deltaX(), y);
625-
double gamma = targetLLH[1];
626-
double alpha = ((bb - aa) * degrees) / (2.0 * _ellipsoid.rEast(gamma) * demInterp.deltaX());
631+
double aa = demInterp.interpolateXY(dem_vect[0] - demInterp.deltaX(), dem_vect[1]);
632+
double bb = demInterp.interpolateXY(dem_vect[0] + demInterp.deltaX(), dem_vect[1]);
633+
634+
Vec3 dem_vect_p_dx = {dem_vect[0] + demInterp.deltaX(), dem_vect[1], dem_vect[2]};
635+
Vec3 dem_vect_m_dx = {dem_vect[0] - demInterp.deltaX(), dem_vect[1], dem_vect[2]};
636+
Vec3 input_coords_llh_p_dx, input_coords_llh_m_dx;
637+
demInterp.proj()->inverse(dem_vect_p_dx, input_coords_llh_p_dx);
638+
demInterp.proj()->inverse(dem_vect_m_dx, input_coords_llh_m_dx);
639+
const Vec3 input_coords_xyz_p_dx = _ellipsoid.lonLatToXyz(input_coords_llh_p_dx);
640+
const Vec3 input_coords_xyz_m_dx = _ellipsoid.lonLatToXyz(input_coords_llh_m_dx);
641+
double dx = (input_coords_xyz_p_dx - input_coords_xyz_m_dx).norm();
642+
643+
double alpha = (bb - aa) / dx;
627644

628645
// North-south slope using central difference
629-
aa = demInterp.interpolateXY(x, y - demInterp.deltaY());
630-
bb = demInterp.interpolateXY(x, y + demInterp.deltaY());
631-
double beta = ((bb - aa) * degrees) / (2.0 * _ellipsoid.rNorth(gamma) * demInterp.deltaY());
646+
aa = demInterp.interpolateXY(dem_vect[0], dem_vect[1] - demInterp.deltaY());
647+
bb = demInterp.interpolateXY(dem_vect[0], dem_vect[1] + demInterp.deltaY());
648+
649+
Vec3 dem_vect_p_dy = {dem_vect[0], dem_vect[1] + demInterp.deltaY(), dem_vect[2]};
650+
Vec3 dem_vect_m_dy = {dem_vect[0], dem_vect[1] - demInterp.deltaY(), dem_vect[2]};
651+
Vec3 input_coords_llh_p_dy, input_coords_llh_m_dy;
652+
demInterp.proj()->inverse(dem_vect_p_dy, input_coords_llh_p_dy);
653+
demInterp.proj()->inverse(dem_vect_m_dy, input_coords_llh_m_dy);
654+
const Vec3 input_coords_xyz_p_dy = _ellipsoid.lonLatToXyz(input_coords_llh_p_dy);
655+
const Vec3 input_coords_xyz_m_dy = _ellipsoid.lonLatToXyz(input_coords_llh_m_dy);
656+
double dy = (input_coords_xyz_p_dy - input_coords_xyz_m_dy).norm();
657+
658+
double beta = (bb - aa) / dy;
632659

633660
// Compute local incidence angle
634661
const Vec3 enunorm = enu.normalized();
@@ -675,6 +702,17 @@ setLayoverShadow(TopoLayers& layers, DEMInterpolator& demInterp,
675702
// Initialize mask to zero for this block
676703
layers.mask() = 0;
677704

705+
// Prepare function getDemCoords() to interpolate DEM
706+
std::function<Vec3(double, double,
707+
const isce3::geometry::DEMInterpolator&,
708+
isce3::core::ProjectionBase*)> getDemCoords;
709+
710+
if (_epsgOut == demInterp.epsgCode()) {
711+
getDemCoords = isce3::geometry::getDemCoordsSameEpsg;
712+
} else {
713+
getDemCoords = isce3::geometry::getDemCoordsDiffEpsg;
714+
}
715+
678716
// Loop over lines in block
679717
#pragma omp parallel for firstprivate(x, y, ctrack, ctrackGrid, \
680718
slantRangeGrid, maskGrid)
@@ -720,12 +758,11 @@ setLayoverShadow(TopoLayers& layers, DEMInterpolator& demInterp,
720758
const double y_grid = y[k] * frac1 + y[k+1] * frac2;
721759

722760
// Interpolate DEM at x/y
723-
const float z_grid = demInterp.interpolateXY(x_grid, y_grid);
761+
Vec3 demXYZ = getDemCoords(x_grid, y_grid, demInterp, _proj);
724762

725763
// Convert DEM XYZ to ECEF XYZ
726764
Vec3 llh, xyz, satToGround;
727-
Vec3 demXYZ{x_grid, y_grid, z_grid};
728-
_proj->inverse(demXYZ, llh);
765+
demInterp.proj()->inverse(demXYZ, llh);
729766
_ellipsoid.lonLatToXyz(llh, xyz);
730767

731768
// Compute and save slant range

‎tests/cxx/isce3/geometry/topo/topo.cpp‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -58,10 +58,13 @@ TEST(TopoTest, RunTopo) {
5858
TEST(TopoTest, CheckResults) {
5959

6060
// Open generated topo raster
61+
std::cout << "test file: ./topo.vrt" << std::endl;
6162
isce3::io::Raster testRaster("topo.vrt");
6263

6364
// Open reference topo raster
64-
isce3::io::Raster refRaster(TESTDATA_DIR "topo/topo.vrt");
65+
std::string ref_filename = TESTDATA_DIR "topo/topo.vrt";
66+
std::cout << "reference file:" << ref_filename << std::endl;
67+
isce3::io::Raster refRaster(ref_filename);
6568

6669
// The associated tolerances
6770
std::vector<double> tols{1.0e-5, 1.0e-5, 0.15, 1.0e-4, 1.0e-4, 0.02, 0.02};
@@ -75,6 +78,9 @@ TEST(TopoTest, CheckResults) {
7578

7679
// Loop over topo bands
7780
for (size_t k = 0; k < refRaster.numBands(); ++k) {
81+
82+
std::cout << "comparing band: " << k + 1 << std::endl;
83+
7884
// Compute sum of absolute error
7985
double error = 0.0;
8086
size_t count = 0;

‎tests/data/topo/localInc.rdr‎

0 Bytes
Binary file not shown.

‎tests/data/topo/localPsi.rdr‎

0 Bytes
Binary file not shown.

‎tests/data/topo_winnipeg/localInc.rdr‎

0 Bytes
Binary file not shown.

‎tests/data/topo_winnipeg/localPsi.rdr‎

0 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)