Skip to content

Commit 66f5b40

Browse files
gshiromagmgunter
authored and
GitHub Enterprise
committed
Substitutes DEM margin units for geocoding from lat/lon degrees to a multiple of DEM pixels (#783)
* fix default DEM margin for DEMs in projected coordinates * add units to the message to the user informing the DEM block margin * add units to the message to the user informing the DEM block margin (2) * remove DEM margin in degrees from GeocodeCov and GeocodeSLC * remove DEM margin in degrees from GeocodeCov and GeocodeSLC (2) * run "git-clang-format develop" * fix rvalue number type (from int to long) * update block margin * remove dem_block_margin from geocode_insar.py * rename pybind_nisar to nisar * rename pybind_nisar to nisar (2) * merge branch develop into update_geocode_interp_dem_margin (2) * rename pybind_nisar to nisar (3) * remove dem_margin/dem_block_marging from geocode_insar.py * merge develop into geocode_metadata (2) * merge develop into update_geocode_interp_dem_margin (3) * merge develop into update_geocode_interp_dem_margin (4) * git merge develop into geogrid_topo (3) * merge develop into update_geocode_interp_dem_margin (5) * merge develop into update_geocode_interp_dem_margin (6) * Update cxx/isce3/geometry/loadDem.h Co-authored-by: Geoffrey M Gunter <[email protected]> * remove unnecessary placeholder function in gslc_runconfig.py Co-authored-by: Geoffrey M Gunter <[email protected]>
1 parent 0e66909 commit 66f5b40

File tree

48 files changed

+66
-321
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+66
-321
lines changed

‎cxx/isce3/cuda/geocode/Geocode.cu‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -178,7 +178,6 @@ __host__
178178
Geocode::Geocode(const isce3::product::GeoGridParameters & geogrid,
179179
const isce3::container::RadarGeometry & rdr_geom,
180180
const Raster & dem_raster,
181-
const double dem_margin,
182181
const size_t lines_per_block,
183182
const isce3::core::dataInterpMethod data_interp_method,
184183
const isce3::core::dataInterpMethod dem_interp_method,
@@ -195,7 +194,6 @@ Geocode::Geocode(const isce3::product::GeoGridParameters & geogrid,
195194
_range_first_pixel(_rdr_geom.radarGrid().width() - 1),
196195
_range_last_pixel(0),
197196
_dem_raster(dem_raster),
198-
_dem_margin(dem_margin),
199197
_interp_float_handle(data_interp_method),
200198
_interp_cfloat_handle(data_interp_method),
201199
_interp_double_handle(data_interp_method),
@@ -269,9 +267,10 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number)
269267
_mask.resize(block_size);
270268

271269
// prepare device DEMInterpolator
270+
int dem_margin_in_pixels = 50;
272271
isce3::geometry::DEMInterpolator host_dem_interp = isce3::geometry::loadDEM(
273272
_dem_raster, _geogrid, _line_start, _geo_block_length,
274-
_geogrid.width(), _dem_margin, _dem_interp_method);
273+
_geogrid.width(), dem_margin_in_pixels, _dem_interp_method);
275274
isce3::cuda::geometry::gpuDEMInterpolator dev_dem_interp(host_dem_interp);
276275
dev_dem_interp.initProjInterp();
277276

‎cxx/isce3/cuda/geocode/Geocode.h‎

Lines changed: 7 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -64,20 +64,18 @@ class Geocode {
6464
* \param[in] geogrid Geogrid defining output product
6565
* \param[in] rdr_geom Radar geometry describing input rasters
6666
* \param[in] dem_raster DEM used to calculate radar grid indices
67-
* \param[in] dem_margin Extra margin applied to bounding box to
68-
* load DEM. Units are need to match
69-
* geogrid EPSG units. \param[in] lines_per_block Number of lines to
70-
* be processed per block \param[in] data_interp_method Data
71-
* interpolation method \param[in] dem_interp_method DEMinterpolation
72-
* method \param[in] threshold Convergence threshold for geo2rdr
67+
* \param[in] lines_per_block Number of lines to be processed per block
68+
* \param[in] data_interp_method Data interpolation method
69+
* \param[in] dem_interp_method DEMinterpolation method
70+
* \param[in] threshold Convergence threshold for geo2rdr
7371
* \param[in] maxiter Maximum iterations for geo2rdr
7472
* \param[in] dr Step size for numerical gradient for
7573
* geo2rdr
7674
* \param[in] invalid_value Value assigned to invalid geogrid pixels
7775
*/
78-
Geocode(const isce3::product::GeoGridParameters& geogrid,
79-
const isce3::container::RadarGeometry& rdr_geom,
80-
const isce3::io::Raster& dem_raster, const double dem_margin,
76+
Geocode(const isce3::product::GeoGridParameters & geogrid,
77+
const isce3::container::RadarGeometry & rdr_geom,
78+
const isce3::io::Raster & dem_raster,
8179
const size_t lines_per_block = 1000,
8280
const isce3::core::dataInterpMethod data_interp_method =
8381
isce3::core::BILINEAR_METHOD,
@@ -152,9 +150,6 @@ class Geocode {
152150
// DEM used to calculate radar grid indices
153151
isce3::io::Raster _dem_raster;
154152

155-
// Extra margin for the dem relative to the geocoded grid
156-
double _dem_margin;
157-
158153
// radar grid boundaries of block last passed to setBlockRdrCoordGrid,
159154
// not for entire geogrid
160155
size_t _az_first_line;

‎cxx/isce3/geocode/GeocodeCov.cpp‎

Lines changed: 17 additions & 122 deletions
Original file line numberDiff line numberDiff line change
@@ -29,10 +29,10 @@ namespace isce3 { namespace geocode {
2929
inline float _getRadarGridOffset(
3030
isce3::core::Matrix<float>& offset_array, float y, float x)
3131
{
32-
const long y_index = std::min(
33-
std::max((long) y, (long) 0), (long) offset_array.length());
34-
const long x_index =
35-
std::min(std::max((long) x, (long) 0), (long) offset_array.width());
32+
const long y_index = std::min(std::max(static_cast<long>(y), 0L),
33+
static_cast<long>(offset_array.length()));
34+
const long x_index = std::min(std::max(static_cast<long>(x), 0L),
35+
static_cast<long>(offset_array.width()));
3636
return offset_array(y_index, x_index);
3737
}
3838

@@ -266,22 +266,15 @@ void Geocode<T>::geocodeInterp(
266266
info << "flatten phase (0:false, 1:true): " << flatten
267267
<< pyre::journal::newline;
268268

269-
if (std::isnan(_demBlockMargin))
270-
// If _demBlockMargin is NaN, add margin of 50 (DEM) pixels
271-
_demBlockMargin =
272-
50.0 * std::max(geogrid.spacingX(), geogrid.spacingY());
273-
274-
info << "DEM block margin " << _demBlockMargin << pyre::journal::newline;
275-
276-
info << "remove phase screen (0: false, 1: true): "
277-
<< std::to_string(phase_screen_raster != nullptr)
278-
<< pyre::journal::newline;
279-
info << "apply azimuth offset (0: false, 1: true): "
280-
<< std::to_string(offset_az_raster != nullptr)
281-
<< pyre::journal::newline;
282-
info << "apply range offset (0: false, 1: true): "
283-
<< std::to_string(offset_rg_raster != nullptr)
284-
<< pyre::journal::newline;
269+
info << "remove phase screen (0: false, 1: true): "
270+
<< std::to_string(phase_screen_raster != nullptr)
271+
<< pyre::journal::newline;
272+
info << "apply azimuth offset (0: false, 1: true): "
273+
<< std::to_string(offset_az_raster != nullptr)
274+
<< pyre::journal::newline;
275+
info << "apply range offset (0: false, 1: true): "
276+
<< std::to_string(offset_rg_raster != nullptr)
277+
<< pyre::journal::newline;
285278

286279
// number of bands in the input raster
287280
int nbands = inputRaster.numBands();
@@ -449,10 +442,12 @@ void Geocode<T>::geocodeInterp(
449442
out_geo_dem_array.fill(std::numeric_limits<float>::quiet_NaN());
450443
}
451444

452-
// load a block of DEM for the current geocoded grid
445+
// load a block of DEM for the current geocoded grid with a margin of
446+
// 50 DEM pixels
447+
int dem_margin_in_pixels = 50;
453448
isce3::geometry::DEMInterpolator demInterp = isce3::geometry::loadDEM(
454449
demRaster, geogrid, lineStart, geoBlockLength, geogrid.width(),
455-
_demBlockMargin, dem_interp_method);
450+
dem_margin_in_pixels, dem_interp_method);
456451

457452
// X and Y indices (in the radar coordinates) for the
458453
// geocoded pixels (after geo2rdr computation)
@@ -898,106 +893,6 @@ void Geocode<T>::_baseband(isce3::core::Matrix<std::complex<T2>>& data,
898893
}
899894
}
900895

901-
template<class T>
902-
void Geocode<T>::_loadDEM(isce3::io::Raster& demRaster,
903-
isce3::geometry::DEMInterpolator& demInterp,
904-
isce3::core::ProjectionBase* proj, int lineStart,
905-
int blockLength, int blockWidth, double demMargin)
906-
{
907-
// Create projection for DEM
908-
int epsgcode = demRaster.getEPSG();
909-
910-
// Initialize bounds
911-
double minX = -1.0e64;
912-
double maxX = 1.0e64;
913-
double minY = -1.0e64;
914-
double maxY = 1.0e64;
915-
916-
// Projection systems are different
917-
if (epsgcode != proj->code()) {
918-
919-
// Create transformer to match the DEM
920-
std::unique_ptr<isce3::core::ProjectionBase> demproj(
921-
isce3::core::createProj(epsgcode));
922-
923-
// Skip factors
924-
const int askip = std::max(static_cast<int>(blockLength / 10.), 1);
925-
const int rskip = std::max(static_cast<int>(blockWidth / 10.), 1);
926-
927-
// Construct vectors of line/pixel indices to traverse perimeter
928-
std::vector<int> lineInd, pixInd;
929-
930-
// Top edge
931-
for (int j = 0; j < blockWidth; j += rskip) {
932-
lineInd.emplace_back(0);
933-
pixInd.emplace_back(j);
934-
}
935-
936-
// Right edge
937-
for (int i = 0; i < blockLength; i += askip) {
938-
lineInd.emplace_back(i);
939-
pixInd.emplace_back(blockWidth);
940-
}
941-
942-
// Bottom edge
943-
for (int j = blockWidth; j > 0; j -= rskip) {
944-
lineInd.emplace_back(blockLength - 1);
945-
pixInd.emplace_back(j);
946-
}
947-
948-
// Left edge
949-
for (int i = blockLength; i > 0; i -= askip) {
950-
lineInd.emplace_back(i);
951-
pixInd.emplace_back(0);
952-
}
953-
954-
// Loop over the indices
955-
for (int i = 0; i < lineInd.size(); ++i) {
956-
Vec3 outpt = {_geoGridStartX + _geoGridSpacingX * (0.5 + pixInd[i]),
957-
_geoGridStartY +
958-
_geoGridSpacingY * (0.5 + lineInd[i]),
959-
0.0};
960-
961-
Vec3 dempt;
962-
if (!projTransform(proj, demproj.get(), outpt, dempt)) {
963-
minX = std::min(minX, dempt[0]);
964-
maxX = std::max(maxX, dempt[0]);
965-
minY = std::min(minY, dempt[1]);
966-
maxY = std::max(maxY, dempt[1]);
967-
}
968-
}
969-
} else {
970-
// Use the corners directly as the projection system is the same
971-
double Y1 = _geoGridStartY + _geoGridSpacingY * lineStart;
972-
double Y2 =
973-
_geoGridStartY + _geoGridSpacingY * (lineStart + blockLength);
974-
minY = std::min(Y1, Y2);
975-
maxY = std::max(Y1, Y2);
976-
minX = _geoGridStartX;
977-
maxX = _geoGridStartX + _geoGridSpacingX * (blockWidth);
978-
}
979-
980-
// If not LonLat, scale to meters
981-
demMargin = (epsgcode != 4326) ? isce3::core::decimaldeg2meters(demMargin)
982-
: demMargin;
983-
984-
// Account for margins
985-
minX -= demMargin;
986-
maxX += demMargin;
987-
minY -= demMargin;
988-
maxY += demMargin;
989-
990-
// load the DEM for this bounding box
991-
demInterp.loadDEM(demRaster, minX, maxX, minY, maxY);
992-
993-
if (demInterp.width() == 0 || demInterp.length() == 0)
994-
std::cout << "warning there are not enough DEM coverage in the "
995-
"bounding box. "
996-
<< std::endl;
997-
998-
// declare the dem interpolator
999-
demInterp.declare();
1000-
}
1001896

1002897
template<class T>
1003898
int Geocode<T>::_geo2rdr(const isce3::product::RadarGridParameters& radar_grid,

‎cxx/isce3/geocode/GeocodeCov.h‎

Lines changed: 0 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -360,11 +360,6 @@ class Geocode {
360360

361361
void linesPerBlock(size_t linesPerBlock) { _linesPerBlock = linesPerBlock; }
362362

363-
void demBlockMargin(double demBlockMargin)
364-
{
365-
_demBlockMargin = demBlockMargin;
366-
}
367-
368363
void radarBlockMargin(int radarBlockMargin)
369364
{
370365
_radarBlockMargin = radarBlockMargin;
@@ -472,11 +467,6 @@ class Geocode {
472467
bool flag_upsample_radar_grid,
473468
geocodeMemoryMode geocode_memory_mode, pyre::journal::info_t& info);
474469

475-
void _loadDEM(isce3::io::Raster& demRaster,
476-
isce3::geometry::DEMInterpolator& demInterp,
477-
isce3::core::ProjectionBase* _proj, int lineStart,
478-
int blockLength, int blockWidth, double demMargin);
479-
480470
std::string _get_nbytes_str(long nbytes);
481471

482472
int _geo2rdr(const isce3::product::RadarGridParameters& radar_grid,
@@ -591,9 +581,6 @@ class Geocode {
591581
// epsg code for the output geogrid
592582
int _epsgOut = 0;
593583

594-
// margin around a computed bounding box for DEM (in degrees)
595-
double _demBlockMargin;
596-
597584
// margin around the computed bounding box for radar dara (integer number of
598585
// lines/pixels)
599586
int _radarBlockMargin;

‎cxx/isce3/geocode/geocodeSlc.cpp‎

Lines changed: 7 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -263,13 +263,12 @@ void geocodeSlc(
263263
const isce3::core::LUT2d<double>& imageGridDoppler,
264264
const isce3::core::Ellipsoid& ellipsoid, const double& thresholdGeo2rdr,
265265
const int& numiterGeo2rdr, const size_t& linesPerBlock,
266-
const double& demBlockMargin, const bool flatten,
267-
const AzRgFunc& azCarrierPhase, const AzRgFunc& rgCarrierPhase,
268-
const std::complex<float> invalidValue)
266+
const bool flatten, const AzRgFunc& azCarrierPhase,
267+
const AzRgFunc& rgCarrierPhase, const std::complex<float> invalidValue)
269268
{
270269
geocodeSlc(outputRaster, inputRaster, demRaster, radarGrid, radarGrid,
271270
geoGrid, orbit,nativeDoppler, imageGridDoppler, ellipsoid,
272-
thresholdGeo2rdr, numiterGeo2rdr, linesPerBlock, demBlockMargin,
271+
thresholdGeo2rdr, numiterGeo2rdr, linesPerBlock,
273272
flatten, azCarrierPhase, rgCarrierPhase, invalidValue);
274273
}
275274

@@ -326,7 +325,7 @@ void geocodeSlc(
326325
const isce3::core::LUT2d<double>& imageGridDoppler,
327326
const isce3::core::Ellipsoid& ellipsoid, const double& thresholdGeo2rdr,
328327
const int& numiterGeo2rdr, const size_t& linesPerBlock,
329-
const double& demBlockMargin, const bool flatten,
328+
const bool flatten,
330329
const AzRgFunc& azCarrierPhase, const AzRgFunc& rgCarrierPhase,
331330
const std::complex<float> invalidValue)
332331
{
@@ -363,8 +362,7 @@ void geocodeSlc(
363362
// get a DEM interpolator for a block of DEM for the current geocoded
364363
// grid
365364
isce3::geometry::DEMInterpolator demInterp = isce3::geometry::loadDEM(
366-
demRaster, geoGrid, lineStart, geoBlockLength, geoGrid.width(),
367-
demBlockMargin);
365+
demRaster, geoGrid, lineStart, geoBlockLength, geoGrid.width());
368366

369367
// X and Y indices (in the radar coordinates) for the
370368
// geocoded pixels (after geo2rdr computation)
@@ -534,7 +532,7 @@ template void geocodeSlc<AzRgFunc>( \
534532
const isce3::core::Ellipsoid& ellipsoid, \
535533
const double& thresholdGeo2rdr, \
536534
const int& numiterGeo2rdr, const size_t& linesPerBlock, \
537-
const double& demBlockMargin, const bool flatten, \
535+
const bool flatten, \
538536
const AzRgFunc& azCarrierPhase, const AzRgFunc& rgCarrierPhase, \
539537
const std::complex<float> invalidValue); \
540538
template void geocodeSlc<AzRgFunc>( \
@@ -549,7 +547,7 @@ template void geocodeSlc<AzRgFunc>( \
549547
const isce3::core::Ellipsoid& ellipsoid, \
550548
const double& thresholdGeo2rdr, \
551549
const int& numiterGeo2rdr, const size_t& linesPerBlock, \
552-
const double& demBlockMargin, const bool flatten, \
550+
const bool flatten, \
553551
const AzRgFunc& azCarrierPhase, const AzRgFunc& rgCarrierPhase, \
554552
const std::complex<float> invalidValue)
555553

‎cxx/isce3/geocode/geocodeSlc.h‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,6 @@ namespace isce3 { namespace geocode {
2424
* \param[in] thresholdGeo2rdr threshold for geo2rdr computations
2525
* \param[in] numiterGeo2rdr maximum number of iterations for Geo2rdr convergence
2626
* \param[in] linesPerBlock number of lines in each block
27-
* \param[in] demBlockMargin margin of a DEM block in degrees
2827
* \param[in] flatten flag to flatten the geocoded SLC
2928
* \param[in] azCarrier azimuth carrier phase of the SLC data, in radians, as a function of azimuth and range
3029
* \param[in] rgCarrier range carrier phase of the SLC data, in radians, as a function of azimuth and range
@@ -40,7 +39,7 @@ void geocodeSlc(isce3::io::Raster& outputRaster, isce3::io::Raster& inputRaster,
4039
const isce3::core::LUT2d<double>& imageGridDoppler,
4140
const isce3::core::Ellipsoid& ellipsoid,
4241
const double& thresholdGeo2rdr, const int& numiterGeo2rdr,
43-
const size_t& linesPerBlock, const double& demBlockMargin,
42+
const size_t& linesPerBlock,
4443
const bool flatten = true,
4544
const AzRgFunc& azCarrier = AzRgFunc(),
4645
const AzRgFunc& rgCarrier = AzRgFunc(),
@@ -66,7 +65,6 @@ void geocodeSlc(isce3::io::Raster& outputRaster, isce3::io::Raster& inputRaster,
6665
* \param[in] thresholdGeo2rdr threshold for geo2rdr computations
6766
* \param[in] numiterGeo2rdr maximum number of iterations for Geo2rdr convergence
6867
* \param[in] linesPerBlock number of lines in each block
69-
* \param[in] demBlockMargin margin of a DEM block in degrees
7068
* \param[in] flatten flag to flatten the geocoded SLC
7169
* \param[in] azCarrier azimuth carrier phase of the SLC data, in radians, as a function of azimuth and range
7270
* \param[in] rgCarrier range carrier phase of the SLC data, in radians, as a function of azimuth and range
@@ -83,7 +81,7 @@ void geocodeSlc(isce3::io::Raster& outputRaster, isce3::io::Raster& inputRaster,
8381
const isce3::core::LUT2d<double>& imageGridDoppler,
8482
const isce3::core::Ellipsoid& ellipsoid,
8583
const double& thresholdGeo2rdr, const int& numiterGeo2rdr,
86-
const size_t& linesPerBlock, const double& demBlockMargin,
84+
const size_t& linesPerBlock,
8785
const bool flatten = true,
8886
const AzRgFunc& azCarrier = AzRgFunc(),
8987
const AzRgFunc& rgCarrier = AzRgFunc(),

‎cxx/isce3/geometry/loadDem.cpp‎

Lines changed: 5 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ namespace isce3 { namespace geometry {
88
DEMInterpolator loadDEM(
99
isce3::io::Raster& demRaster,
1010
const isce3::product::GeoGridParameters& geoGrid, int lineStart,
11-
int blockLength, int blockWidth, double demMargin,
11+
int blockLength, int blockWidth, const int demMarginInPixels,
1212
isce3::core::dataInterpMethod demInterpMethod)
1313
{
1414
// DEM interpolator
@@ -90,15 +90,11 @@ DEMInterpolator loadDEM(
9090
maxX = geoGrid.startX() + geoGrid.spacingX() * (blockWidth - 1);
9191
}
9292

93-
// If not LonLat, scale to meters
94-
demMargin = (epsgcode != 4326) ? isce3::core::decimaldeg2meters(demMargin)
95-
: demMargin;
96-
9793
// Account for margins
98-
minX -= demMargin;
99-
maxX += demMargin;
100-
minY -= demMargin;
101-
maxY += demMargin;
94+
minX -= demMarginInPixels * demRaster.dx();
95+
maxX += demMarginInPixels * demRaster.dx();
96+
minY -= demMarginInPixels * std::abs(demRaster.dy());
97+
maxY += demMarginInPixels * std::abs(demRaster.dy());
10298

10399
std::cout << minX << " , " << maxX << " , " << minY << ", " << maxY
104100
<< std::endl;

0 commit comments

Comments
 (0)