From 4cd3505903a1686bb3f6eb728b0e3c6e701623cd Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Mon, 25 Sep 2023 11:59:39 -0700 Subject: [PATCH] Aztimuth and slant range corrections for CUDA geocode (#1423) * az and srg corrections for CUDA geocode * corrections incorporated into masking geogrid radar indices * update geocode function call * workflow CUDA geocode object init consolidated inside gpu_geocode_rasters * more comprehensive unit tests * output to temporary directory * check number of invalids in tests * removed redundant lines_per_block --- cxx/isce3/cuda/geocode/Geocode.cu | 150 +++++-- cxx/isce3/cuda/geocode/Geocode.h | 49 ++- cxx/isce3/geocode/geocodeSlc.cpp | 6 +- .../pybind_isce3/cuda/geocode/cuGeocode.cpp | 12 + .../packages/nisar/workflows/geocode_insar.py | 42 +- .../extensions/pybind/cuda/geocode/geocode.py | 408 ++++++++++++++---- .../packages/isce3/geocode/geocode_slc.py | 1 - 7 files changed, 514 insertions(+), 154 deletions(-) diff --git a/cxx/isce3/cuda/geocode/Geocode.cu b/cxx/isce3/cuda/geocode/Geocode.cu index 6e7868c8d..a26863b8c 100644 --- a/cxx/isce3/cuda/geocode/Geocode.cu +++ b/cxx/isce3/cuda/geocode/Geocode.cu @@ -12,7 +12,6 @@ #include #include #include -#include #include #include #include @@ -45,21 +44,24 @@ namespace isce3::cuda::geocode { /** Coverts a batch of rows from input geogrid into radar coordinates. Outputs * the pixel-space coordinates (x,y) of each resulting (range, azimuth) pair * with respect to specified radargrid, as well as mask invalid pixels (out - * of bounds or failed to converge in geo2rdr). + * of bounds, failed to converge in geo2rdr, or not contained in + * azTimeCorrections, sRangeCorrections, or nativeDoppler). * * \param[out] rdr_x pointer to device_vector of computed radar grid * x / slant range indices * \param[out] rdr_y pointer to device_vector of computed radar grid * y / azimuth time indices - * \param[out] mask pointer to device_vector mask of valid - * pixels. Mask follows numpy mask_array - * convention where True is masked + * \param[out] mask pointer to device_vector mask of valid pixels + * Mask follows numpy mask_array convention where + * True is masked. * \param[in] ellipsoid Ellipsoid based on output geogrid coordinate * system * \param[in] orbit Orbit associated with radar data * \param[in] dem DEM interpolator. Maybe be of different * coordinate system than output geogrid. - * \param[in] doppler doppler + * \param[in] nativeDoppler doppler centroid of data associated with radar + * grid, in Hz, as a function of azimuth and range + * \param[in] imageDoppler image grid doppler * \param[in] wvl wavelength * \param[in] side look side * \param[in] geo2rdr_params geo2rdr params @@ -69,15 +71,23 @@ namespace isce3::cuda::geocode { * \param[in] block_size Number of elements in a block * \param[in] proj Projection used to covert geogrid XYZ to LLH * of output coordinate system + * \param[in] azTimeCorrection geo2rdr azimuth additive correction, in + * seconds, as a function of azimuth and range + * \param[in] sRangeCorrection geo2rdr slant range additive correction, in + * seconds, as a function of azimuth and range */ __global__ void geoToRdrIndices(double* rdr_x, double* rdr_y, bool* mask, const isce3::core::Ellipsoid ellipsoid, const DeviceOrbitView orbit, isce3::cuda::geometry::gpuDEMInterpolator dem, - const DeviceLUT2d doppler, const double wvl, + const DeviceLUT2d nativeDoppler, + const DeviceLUT2d imageDoppler, + const double wvl, const isce3::core::LookSide side, const Geo2RdrParams geo2rdr_params, const isce3::product::GeoGridParameters geogrid, const RadarGridParamsLite radargrid, const size_t line_start, - const size_t block_size, isce3::cuda::core::ProjectionBase** proj) + const size_t block_size, isce3::cuda::core::ProjectionBase** proj, + const DeviceLUT2d azTimeCorrection, + const DeviceLUT2d sRangeCorrection) { // thread index (1d grid of 1d blocks) const auto tid = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; @@ -97,25 +107,47 @@ __global__ void geoToRdrIndices(double* rdr_x, double* rdr_y, bool* mask, (*proj)->inverse(xyz, llh); llh[2] = dem.interpolateLonLat(llh[0], llh[1]); - // geo2rdr slant range - double r; - // geo2rdr azimuth time initial value - double t = radargrid.sensing_mid; + // slant range to be computed by geo2rdr + double srange; + // initial guess for azimuth time to be updated by geo2rdr + double aztime = radargrid.sensing_mid; // returns 0 if geo2rdr converges else 1 int converged = isce3::cuda::geometry::geo2rdr(llh, ellipsoid, orbit, - doppler, &t, &r, wvl, side, geo2rdr_params.threshold, - geo2rdr_params.maxiter, geo2rdr_params.delta_range); + imageDoppler, &aztime, &srange, wvl, side, + geo2rdr_params.threshold, geo2rdr_params.maxiter, + geo2rdr_params.delta_range); + + // Bool to track if aztime and slant range is not contained in any of the + // following LUT2d's: azTimeCorrection, sRangeCorrection, and nativeDoppler. + // When set to true, this geogrid point will be marked as invalid + bool not_in_lut = false; + + // Apply timing corrections + if (azTimeCorrection.contains(aztime, srange)) { + const auto aztimeCor = azTimeCorrection.eval(aztime, srange); + aztime += aztimeCor; + } else + not_in_lut = true; + + if (sRangeCorrection.contains(aztime, srange)) { + const auto srangeCor = sRangeCorrection.eval(aztime, srange); + srange += srangeCor; + } else + not_in_lut = true; + + if (!nativeDoppler.contains(aztime, srange)) + not_in_lut = true; // convert aztime and range to indices - double y = (t - radargrid.sensing_start) * radargrid.prf; - double x = (r - radargrid.starting_range) / radargrid.range_pixel_spacing; + double y = (aztime - radargrid.sensing_start) * radargrid.prf; + double x = (srange - radargrid.starting_range) / radargrid.range_pixel_spacing; - // check if indinces in bounds and set accordingly + // check if indinces in bounds of radar grid const bool not_in_rdr_grid = y < 0 || y >= radargrid.length || x < 0 || x >= radargrid.width; - const bool invalid_index = not_in_rdr_grid || converged == 0; + const bool invalid_index = not_in_rdr_grid || converged == 0 || not_in_lut; rdr_y[tid] = invalid_index ? 0.0 : y; rdr_x[tid] = invalid_index ? 0.0 : x; mask[tid] = invalid_index; @@ -287,8 +319,7 @@ void invalidValueCheck(const int dtype, const double invalid_value) __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid, const isce3::container::RadarGeometry& rdr_geom, const size_t lines_per_block) : - _lines_per_block(lines_per_block), - _n_blocks((geogrid.length() + _lines_per_block - 1) / _lines_per_block), + _n_blocks((geogrid.length() + lines_per_block - 1) / lines_per_block), _geogrid(geogrid), _rdr_geom(rdr_geom), _ellipsoid(isce3::core::makeProjection(_geogrid.epsg())->ellipsoid()), @@ -298,23 +329,23 @@ __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid, // correctly _az_first_line(_rdr_geom.radarGrid().length() - 1), _az_last_line(0), - _range_first_pixel(_rdr_geom.radarGrid().width() - 1), + _range_first_pixel(rdr_geom.radarGrid().width() - 1), _range_last_pixel(0), - _geo_block_length(_lines_per_block), + _geo_block_length(lines_per_block), _proj_handle(geogrid.epsg()) { // init light weight radar grid - _radar_grid.sensing_start = _rdr_geom.radarGrid().sensingStart(); - _radar_grid.sensing_mid = _rdr_geom.radarGrid().sensingMid(); - _radar_grid.prf = _rdr_geom.radarGrid().prf(); - _radar_grid.starting_range = _rdr_geom.radarGrid().startingRange(); - _radar_grid.range_pixel_spacing = _rdr_geom.radarGrid().rangePixelSpacing(); - _radar_grid.length = _rdr_geom.gridLength(); - _radar_grid.width = _rdr_geom.gridWidth(); + _radar_grid.sensing_start = rdr_geom.radarGrid().sensingStart(); + _radar_grid.sensing_mid = rdr_geom.radarGrid().sensingMid(); + _radar_grid.prf = rdr_geom.radarGrid().prf(); + _radar_grid.starting_range = rdr_geom.radarGrid().startingRange(); + _radar_grid.range_pixel_spacing = rdr_geom.radarGrid().rangePixelSpacing(); + _radar_grid.length = rdr_geom.gridLength(); + _radar_grid.width = rdr_geom.gridWidth(); // Determine max number of elements per block. Last block to be processed // may not contain the max number of elements. - auto n_elem = _lines_per_block * _geogrid.width(); + auto n_elem = _geo_block_length * _geogrid.width(); // Resize all device vectors to max block size. _radar_x.resize(n_elem); @@ -324,8 +355,12 @@ __host__ Geocode::Geocode(const isce3::product::GeoGridParameters& geogrid, void Geocode::setBlockRdrCoordGrid(const size_t block_number, Raster& dem_raster, - const isce3::core::dataInterpMethod _dem_interp_method, - const isce3::geometry::detail::Geo2RdrParams _geo2rdr_params) + const isce3::core::dataInterpMethod dem_interp_method, + const isce3::cuda::container::RadarGeometry& dev_rdr_geom, + const isce3::geometry::detail::Geo2RdrParams geo2rdr_params, + const DeviceLUT2d& nativeDoppler, + const DeviceLUT2d& azTimeCorrection, + const DeviceLUT2d& sRangeCorrection) { // make sure block index does not exceed actual number of blocks if (block_number >= _n_blocks) { @@ -334,11 +369,10 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number, } // Get block extents (of the geocoded grid) - _line_start = block_number * _lines_per_block; + _line_start = block_number * _geo_block_length; // Set block sizes for everything but last block size_t block_size = _geo_block_length * _geogrid.width(); - _geo_block_length = _lines_per_block; // Adjust for last block sized assuming it is sized differently than others if (block_number == (_n_blocks - 1)) { @@ -356,12 +390,9 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number, isce3::geometry::DEMInterpolator host_dem_interp = isce3::geometry::DEMRasterToInterpolator(dem_raster, _geogrid, _line_start, _geo_block_length, _geogrid.width(), - dem_margin_in_pixels, _dem_interp_method); + dem_margin_in_pixels, dem_interp_method); isce3::cuda::geometry::gpuDEMInterpolator dev_dem_interp(host_dem_interp); - // copy RadarGeometry to device - isce3::cuda::container::RadarGeometry dev_rdr_geom(_rdr_geom); - // Create geogrid on device { const unsigned threads_per_block = 256; @@ -370,16 +401,20 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number, geoToRdrIndices<<>>(_radar_x.data().get(), _radar_y.data().get(), _mask.data().get(), _ellipsoid, - dev_rdr_geom.orbit(), dev_dem_interp, dev_rdr_geom.doppler(), + dev_rdr_geom.orbit(), dev_dem_interp, + nativeDoppler, dev_rdr_geom.doppler(), dev_rdr_geom.wavelength(), dev_rdr_geom.lookSide(), - _geo2rdr_params, _geogrid, _radar_grid, _line_start, block_size, - _proj_handle.get_proj()); + geo2rdr_params, _geogrid, _radar_grid, _line_start, block_size, + _proj_handle.get_proj(), azTimeCorrection, sRangeCorrection); checkCudaErrors(cudaPeekAtLastError()); checkCudaErrors(cudaDeviceSynchronize()); } - // find index of min and max values in x/range + // Mask is utilized to _radar_x and _radar_y when determining respective + // min and max values which are then used to determine corresponding radar + // data to read from radar raster + // Find index of unmasked min and max values in x/range const auto [rdr_x_min, rdr_x_max] = masked_minmax(_radar_x, _mask); _range_first_pixel = std::min(_rdr_geom.radarGrid().width() - 1, @@ -388,7 +423,7 @@ void Geocode::setBlockRdrCoordGrid(const size_t block_number, _range_last_pixel = std::max(static_cast(0), static_cast(std::ceil(rdr_x_max) - 1)); - // find index of min and max values in y/azimuth + // Find index of unmasked min and max values in y/azimuth const auto [rdr_y_min, rdr_y_max] = masked_minmax(_radar_y, _mask); _az_first_line = std::min(_rdr_geom.radarGrid().length() - 1, @@ -586,6 +621,9 @@ void Geocode::geocodeRasters( const std::vector& raster_datatypes, const std::vector& invalid_values_double, Raster& dem_raster, + const isce3::core::LUT2d& hostNativeDoppler, + const isce3::core::LUT2d& hostAzTimeCorrection, + const isce3::core::LUT2d& hostSRangeCorrection, const isce3::core::dataInterpMethod dem_interp_method, const double threshold, const int maxiter, const double dr) { @@ -603,7 +641,18 @@ void Geocode::geocodeRasters( std::vector _invalid_values; // geo2rdr params used in radar index calculation - isce3::geometry::detail::Geo2RdrParams _geo2rdr_params{threshold, maxiter, dr}; + isce3::geometry::detail::Geo2RdrParams geo2rdr_params{threshold, maxiter, dr}; + + // 2D LUT Doppler of the SLC image + auto nativeDoppler = DeviceLUT2d(hostNativeDoppler); + + // geo2rdr azimuth additive correction, in seconds, as a function of + // azimuth and range + auto azTimeCorrection = DeviceLUT2d(hostAzTimeCorrection);; + + // geo2rdr slant range additive correction, in seconds, as a function of + // azimuth and range + auto sRangeCorrection = DeviceLUT2d(hostSRangeCorrection); // Prepare interpolator handler and invalid value for each raster for (size_t i_raster = 0; i_raster < n_raster_pairs; ++i_raster) { @@ -680,12 +729,21 @@ void Geocode::geocodeRasters( } } + // copy RadarGeometry to device + const auto dev_rdr_geom = isce3::cuda::container::RadarGeometry(_rdr_geom); + // iterate over blocks for (size_t i_block = 0; i_block < _n_blocks; ++i_block) { // compute then set radar coords for each geocode obj for curret block - setBlockRdrCoordGrid(i_block, dem_raster, dem_interp_method, - _geo2rdr_params); + setBlockRdrCoordGrid(i_block, + dem_raster, + dem_interp_method, + dev_rdr_geom, + geo2rdr_params, + nativeDoppler, + azTimeCorrection, + sRangeCorrection); // iterate over rasters and geocode with associated datatype and // interpolation method diff --git a/cxx/isce3/cuda/geocode/Geocode.h b/cxx/isce3/cuda/geocode/Geocode.h index acad8edcc..63981bd23 100644 --- a/cxx/isce3/cuda/geocode/Geocode.h +++ b/cxx/isce3/cuda/geocode/Geocode.h @@ -11,9 +11,11 @@ #include #include +#include #include #include #include +#include #include #include #include @@ -21,6 +23,9 @@ namespace isce3::cuda::geocode { +template +using DeviceLUT2d = isce3::cuda::core::gpuLUT2d; + /* light weight radar grid container */ struct RadarGridParamsLite { double sensing_start; @@ -83,6 +88,15 @@ class Geocode { * initialize templated interpolators. * \param[in] invalid_values Invalid values for each raster * \param[in] dem_raster DEM used to calculate radar grid indices + * \param[in] nativeDoppler Doppler centroid of data in Hz associated + * radar grid, as a function azimuth and + * range + * \param[in] azTimeCorrection geo2rdr azimuth additive correction, in + * seconds, as a function of azimuth and + * range + * \param[in] sRangeCorrection geo2rdr slant range additive + * correction, in seconds, as a function + * of azimuth and range * \param[in] dem_interp_method DEMinterpolation method * \param[in] threshold Convergence threshold for geo2rdr * \param[in] maxiter Maximum iterations for geo2rdr @@ -96,15 +110,16 @@ class Geocode { const std::vector& raster_datatypes, const std::vector& invalid_values, Raster& dem_raster, + const isce3::core::LUT2d& hostNativeDoppler = {}, + const isce3::core::LUT2d& hostAzTimeCorrection = {}, + const isce3::core::LUT2d& hostSRangeCorrection = {}, const isce3::core::dataInterpMethod dem_interp_method = isce3::core::BIQUINTIC_METHOD, - const double threshold = 1e-8, const int maxiter = 50, + const double threshold = 1e-8, + const int maxiter = 50, const double dr = 10); private: - // number of lines to be processed in a block - size_t _lines_per_block; - // total number of blocks necessary to geocoding a provided geogrid size_t _n_blocks; @@ -156,7 +171,8 @@ class Geocode { * \param[in] raster_datatypes GDAL type of each raster represented by * equivalent int. Necessary to correctly * initialize templated interpolators. - * \param[in] invalid_values Invalid values for each raster + * \param[in] invalid_values_double Invalid values for each raster as double, to + be converted later to the correct type */ void ensureRasterConsistency( const std::vector>& output_rasters, @@ -173,14 +189,29 @@ class Geocode { * \param[in] block_number Index of block where radar grid coordinates * coordinates will be calculated and set * \param[in] dem_raster DEM used to calculate radar grid indices - * \param[in] _dem_interp_method DEMinterpolation method - * \param[in] _geo2rdr_params Parameters used by geo2rdr to compute + * \param[in] dem_interp_method DEMinterpolation method + * \param[in] dev_rdr_geom Radar geometry describing input raster + * stored on-device + * \param[in] geo2rdr_params Parameters used by geo2rdr to compute * geogrids radar grid coordinates + * \param[in] nativeDoppler Doppler centroid of data in Hz associated + * radar grid, as a function azimuth and + * range + * \param[in] azTimeCorrection geo2rdr azimuth additive correction, in + * seconds, as a function of azimuth and + * range + * \param[in] sRangeCorrection geo2rdr slant range additive + * correction, in seconds, as a function + * of azimuth and range */ void setBlockRdrCoordGrid(const size_t block_number, Raster& dem_raster, - const isce3::core::dataInterpMethod _dem_interp_method, - const isce3::geometry::detail::Geo2RdrParams _geo2rdr_params); + const isce3::core::dataInterpMethod dem_interp_method, + const isce3::cuda::container::RadarGeometry& dev_rdr_geom, + const isce3::geometry::detail::Geo2RdrParams geo2rdr_params, + const DeviceLUT2d& nativeDoppler, + const DeviceLUT2d& azTimeCorrection, + const DeviceLUT2d& sRangeCorrection); /** Geocode a block of raster according to grid last set in * setBlockRdrCoordGrid. Only operates on 1st raster band of input/output diff --git a/cxx/isce3/geocode/geocodeSlc.cpp b/cxx/isce3/geocode/geocodeSlc.cpp index 783fd447e..1a06fc1bf 100644 --- a/cxx/isce3/geocode/geocodeSlc.cpp +++ b/cxx/isce3/geocode/geocodeSlc.cpp @@ -512,9 +512,9 @@ void geocodeSlc( uncorrectedSRange(blockLine, pixel) = srange; // apply timing corrections - // if default LUT2d used for both azimuth time and slant range - // corrections, uncorrected slant range and corrected slant - // range will be the same + // If default LUT2d is passed in for either azimuth time and + // slant range corrections, eval() returns 0 so no corrections + // are applied. if (azTimeCorrection.contains(aztime, srange)) { const auto aztimeCor = azTimeCorrection.eval(aztime, srange); diff --git a/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp b/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp index 5c5f83e13..3d771426c 100644 --- a/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp +++ b/python/extensions/pybind_isce3/cuda/geocode/cuGeocode.cpp @@ -43,6 +43,9 @@ void addbinding(pybind11::class_& pyGeocode) py::arg("raster_datatypes"), py::arg("invalid_values"), py::arg("dem_raster"), + py::arg("native_doppler") = isce3::core::LUT2d(), + py::arg("az_time_correction") = isce3::core::LUT2d(), + py::arg("srange_correction") = isce3::core::LUT2d(), py::arg("dem_interp_method") = isce3::core::BIQUINTIC_METHOD, py::arg("threshold") = g2r_defaults.threshold, @@ -65,6 +68,15 @@ void addbinding(pybind11::class_& pyGeocode) Invalid values for each geocoded raster dem_raster: isce3.io.Raster DEM used to calculate radar grid indices + native_doppler: isce3.core.LUT2d + Doppler centroid of data associated with radar grid, in Hz, as + fuction of azimuth and range + az_time_correction: isce3.core.LUT2d + geo2rdr azimuth additive correction, in seconds, as a function + of azimuth and range + srange_correction: isce3.core.LUT2d + geo2rdr slant range additive correction, in seconds, as a + function of azimuth and range dem_interp_method: isce3.core.DataInterpMethod Interpolation method used by DEM interpolator. Default BIQUINTIC_METHOD diff --git a/python/packages/nisar/workflows/geocode_insar.py b/python/packages/nisar/workflows/geocode_insar.py index 587c77be3..7e69c0c03 100644 --- a/python/packages/nisar/workflows/geocode_insar.py +++ b/python/packages/nisar/workflows/geocode_insar.py @@ -363,7 +363,7 @@ def add_radar_grid_cube(cfg, freq, radar_grid, orbit, dst_h5, input_product_type grid_zero_doppler = isce3.core.LUT2d() # The native-Doppler LUT bounds error is turned off to - # computer cubes values outside radar-grid boundaries + # compute cubes values outside radar-grid boundaries native_doppler.bounds_error = False add_radar_grid_cubes_to_hdf5(dst_h5, cube_group_path, cube_geogrid_param, radar_grid_cubes_heights, @@ -742,15 +742,13 @@ def cpu_run(cfg, input_hdf5, output_hdf5, input_product_type=InputProduct.RUNW): geocode_iono_bool = True input_hdf5_iono = input_hdf5 if is_iono_method_sideband and freq == 'A': - ''' - ionosphere_phase_screen from main_side_band or - main_diff_ms_band are computed on radargrid of frequencyB. - The ionosphere_phase_screen is geocoded on geogrid of - frequencyA. Instead of geocoding ionosphere in the RUNW - standard product (frequencyA), geocode the frequencyB in - scratch/ionosphere/method/RUNW.h5 to avoid additional - interpolation. - ''' + # ionosphere_phase_screen from main_side_band or + # main_diff_ms_band are computed on radargrid of frequencyB. + # The ionosphere_phase_screen is geocoded on geogrid of + # frequencyA. Instead of geocoding ionosphere in the RUNW + # standard product (frequencyA), geocode the frequencyB in + # scratch/ionosphere/method/RUNW.h5 to avoid additional + # interpolation. radar_grid_iono = slc.getRadarGrid('B') iono_sideband_bool = True if az_looks > 1 or rg_looks > 1: @@ -879,12 +877,20 @@ def cpu_run(cfg, input_hdf5, output_hdf5, input_product_type=InputProduct.RUNW): t_all_elapsed = time.time() - t_all info_channel.log(f"Successfully ran geocode in {t_all_elapsed:.3f} seconds") -def gpu_geocode_rasters(geocoded_dataset_flags, desired_geo_dataset_names, - interpolation_methods, invalid_values, - freq, pol_list, - geogrid, rdr_geometry, dem_raster, lines_per_block, - input_hdf5, dst_h5, - offset_layers=None, scratch_path='', +def gpu_geocode_rasters(geocoded_dataset_flags, + desired_geo_dataset_names, + interpolation_methods, + invalid_values, + freq, + pol_list, + geogrid, + rdr_geometry, + dem_raster, + lines_per_block, + input_hdf5, + dst_h5, + offset_layers=None, + scratch_path='', compute_stats=True, input_product_type=InputProduct.RUNW, iono_sideband=False): @@ -1309,9 +1315,7 @@ def gpu_run(cfg, input_hdf5, output_hdf5, input_product_type=InputProduct.RUNW): if __name__ == "__main__": - """ - run geocode from command line - """ + # run geocode from command line # load command line args geocode_insar_parser = YamlArgparse() diff --git a/tests/python/extensions/pybind/cuda/geocode/geocode.py b/tests/python/extensions/pybind/cuda/geocode/geocode.py index a83d75687..13591ddb8 100644 --- a/tests/python/extensions/pybind/cuda/geocode/geocode.py +++ b/tests/python/extensions/pybind/cuda/geocode/geocode.py @@ -1,107 +1,363 @@ #/usr/bin/env python3 import itertools +import json import os +import types import numpy as np from osgeo import gdal -import iscetest import isce3 +from isce3.atmosphere.tec_product import tec_lut2d_from_json +from isce3.geometry import compute_incidence_angle +import iscetest from nisar.products.readers import SLC +import pytest +from scipy import interpolate + -def test_cuda_geocode(): - rslc = SLC(hdf5file=os.path.join(iscetest.data, "envisat.h5")) +@pytest.fixture(scope='module') +def unit_test_params(tmp_path_factory): + ''' + test parameters shared by all geocode tests + ''' + # load h5 for doppler and orbit + params = types.SimpleNamespace() - dem_raster = isce3.io.Raster(os.path.join(iscetest.data, - "geocode/zeroHeightDEM.geo")) + # create temporary output directory for test output + params.out_dir = tmp_path_factory.mktemp('test_output') # define geogrid - epsg = 4326 geogrid = isce3.product.GeoGridParameters(start_x=-115.65, - start_y=34.84, - spacing_x=0.0002, - spacing_y=-8.0e-5, - width=500, - length=500, - epsg=epsg) - - geotrans = [geogrid.start_x, geogrid.spacing_x, 0.0, - geogrid.start_y, 0.0, geogrid.spacing_y] - - # Init RadarGeometry, orbit, and doppler from RSLC - radargrid = isce3.product.RadarGridParameters(os.path.join(iscetest.data, - "envisat.h5")) - orbit = rslc.getOrbit() - doppler = rslc.getDopplerCentroid() - rdr_geometry = isce3.container.RadarGeometry(radargrid, - orbit, - doppler) + start_y=34.84, + spacing_x=0.0002, + spacing_y=-8.0e-5, + width=500, + length=500, + epsg=4326) + + params.geogrid = geogrid + + # define geotransform + params.geotrans = [geogrid.start_x, geogrid.spacing_x, 0.0, + geogrid.start_y, 0.0, geogrid.spacing_y] + + input_h5_path = os.path.join(iscetest.data, "envisat.h5") + + params.radargrid = isce3.product.RadarGridParameters(input_h5_path) + + # init SLC object and extract necessary test params from it + rslc = SLC(hdf5file=input_h5_path) + + params.orbit = rslc.getOrbit() + + # XXX This attribute of NISAR RSLCs is supposed to be the image grid Doppler + # but for the envisat.h5 file, it's commonly used in isce3 unit tests to + # represent native Doppler. + img_doppler = rslc.getDopplerCentroid() + params.img_doppler = img_doppler + + params.center_freq = rslc.getSwathMetadata().processed_center_frequency + + params.native_doppler = isce3.core.LUT2d(img_doppler.x_start, + img_doppler.y_start, img_doppler.x_spacing, + img_doppler.y_spacing, np.zeros((geogrid.length,geogrid.width))) + + # create DEM raster object + params.dem_path = os.path.join(iscetest.data, "geocode/zeroHeightDEM.geo") + params.dem_raster = isce3.io.Raster(params.dem_path) + + # half pixel offset and grid size in radians for validataion + params.x0 = np.radians(geogrid.start_x + geogrid.spacing_x / 2.0) + params.dx = np.radians(geogrid.spacing_x) + params.y0 = np.radians(geogrid.start_y + geogrid.spacing_y / 2.0) + params.dy = np.radians(geogrid.spacing_y) + + # multiplicative factor applied to range pixel spacing and azimuth time + # interval to be added to starting range and azimuth time of radar grid + params.offset_factor = 10.0 + + # TEC JSON containing TEC values that generate range offsets that match the + # fixed range offset used to test range correction + params.tec_json_path = params.out_dir / 'test_tec.json' + make_tec_file(params) + # compute unmasked validation arrays based on the geogrid being geocoded to + # params.grid_lon - array where each column is the longitude of the geogrid + # column + # params.grid_lat - array where each row is the latitude of the geogrid row + nx = geogrid.width + ny = geogrid.length + meshx, meshy = np.meshgrid(np.arange(nx), np.arange(ny)) + params.grid_lon = params.x0 + meshx * params.dx + params.grid_lat = params.y0 + meshy * params.dy + + return params + + +def make_tec_file(unit_test_params): + ''' + create TEC file using radar grid from envisat.h5 that yields a uniform + slant range offset when processed with tec_lut2d_from_json() + We ignore topside TEC and simulate total TEC at near and far ranges such + that the slant range delay at near and far ranges are the same. + solve for sub_orbital_tec from: + delta_r = K * sub_orbital_tec * TECU / center_freq**2 / np.cos(incidence) + + yields: + sub_orbital_tec = delta_r * np.cos(incidence) * center_freq**2 / (TECU * K) + ''' + radar_grid = unit_test_params.radargrid + + # create linspace for radar grid sensing time + t_rdr_grid = np.linspace(radar_grid.sensing_start, + radar_grid.sensing_stop + 1.0 / radar_grid.prf, + radar_grid.length) + + # TEC coefficients + K = 40.31 # its a constant in m3/s2 + TECU = 1e16 # its a constant to convert the TEC product to electrons / m2 + + # set delta_r to value used to test slant range offset correction in + # geocode_test_cases() + offset_factor = unit_test_params.offset_factor + delta_r = offset_factor * radar_grid.range_pixel_spacing + + # compute common TEC coefficient used for both near and far TEC + common_tec_coeff = delta_r * unit_test_params.center_freq**2 / (K * TECU) + + # get TEC times in ISO format + # +/- 50 sec from stop/start of radar grid + margin = 50. + # 10 sec increments - also snap to multiples of 10 sec + snap = 10. + start = np.floor(radar_grid.sensing_start / snap) * snap - margin + stop = np.ceil(radar_grid.sensing_stop / snap) * snap + margin + t_tec = np.arange(start, stop + 1.0, snap) + t_tec_iso_fmt = [(radar_grid.ref_epoch + isce3.core.TimeDelta(t)).isoformat()[:-3] + for t in t_tec] + + # compute total TEC + total_tec = [] + for rdr_grid_range in [radar_grid.starting_range, + radar_grid.end_range]: + inc_angs = [compute_incidence_angle(t, rdr_grid_range, + unit_test_params.orbit, + isce3.core.LUT2d(), + radar_grid, + isce3.geometry.DEMInterpolator(), + isce3.core.Ellipsoid()) + for t in t_rdr_grid] + total_tec_rdr_grid = common_tec_coeff * np.cos(inc_angs) + + # near and far top TEC = 0 to allow sub orbital TEC = total TEC + # create extraplotor/interpolators for near and far + total_tec_interp = interpolate.interp1d(t_rdr_grid, total_tec_rdr_grid, + 'linear', + fill_value="extrapolate") + + # compute near and far total TEC + total_tec.append(total_tec_interp(t_tec)) + total_tec_near, total_tec_far = total_tec + + # load relevant TEC into dict and write to JSON + # top TEC = 0 to allow sub orbital TEC = total TEC + tec_zeros = list(np.zeros(total_tec_near.shape)) + tec_dict ={} + tec_dict['utc'] = t_tec_iso_fmt + tec_dict['totTecNr'] = list(total_tec_near) + tec_dict['topTecNr'] = tec_zeros + tec_dict['totTecFr'] = list(total_tec_far) + tec_dict['topTecFr'] = tec_zeros + with open(str(unit_test_params.tec_json_path), 'w') as fp: + json.dump(tec_dict, fp) + + +def geocode_test_cases(unit_test_params): + ''' + Generator for CUDA geocode test cases + + Given a radar grid, generate correction LUT2ds in range and azimuth + directions. Returns axis, offset mode name, range and azimuth correction + LUT2ds and offset corrected radar grid. + ''' + test_case = types.SimpleNamespace() + + radargrid = unit_test_params.radargrid + offset_factor = unit_test_params.offset_factor + + rg_pxl_spacing = radargrid.range_pixel_spacing + range_offset = offset_factor * rg_pxl_spacing + az_time_interval = 1 / radargrid.prf + azimuth_offset = offset_factor * az_time_interval + + # despite uniform value LUT2d set to interp mode nearest just in case + method = isce3.core.DataInterpMethod.NEAREST + + # array of ones to be multiplied by respective offset value + # shape unchanging; no need to be in loop as only starting values change + ones = np.ones(radargrid.shape) + + axes = 'xy' + test_modes = ['', 'rg', 'az', 'rg_az', 'tec', 'blocked'] + for axis, test_mode in itertools.product(axes, test_modes): + is_block_mode = test_mode == 'blocked' + + # skip blocked if for y-axis. x-axis testing sufficient + if axis == 'y' and is_block_mode: + continue + + # if testing block mode, lines_per_block is chosen to be less than the + # total number of lines in the dataset + test_case.lines_per_block = 126 if is_block_mode else 1000 + test_case.axis = axis + test_case.test_mode = test_mode + + # create radar and apply positive offsets in range and azimuth + test_case.radargrid = radargrid.copy() + + # apply offsets as required by mode + if 'rg' in test_mode or 'tec' == test_mode: + test_case.radargrid.starting_range += range_offset + if 'az' in test_mode: + test_case.radargrid.sensing_start += azimuth_offset + + test_case.rdr_geometry = \ + isce3.container.RadarGeometry(test_case.radargrid, + unit_test_params.orbit, + unit_test_params.img_doppler) + + # slant range vector for LUT2d + srange_vec = np.linspace(test_case.radargrid.starting_range, + test_case.radargrid.end_range, + radargrid.width) + + # azimuth vector for LUT2d + az_time_vec = np.linspace(test_case.radargrid.sensing_start, + test_case.radargrid.sensing_stop, + radargrid.length) + + # corrections LUT2ds will use the negative offsets + # should cancel positive offset applied to radar grid + srange_correction = isce3.core.LUT2d() + if 'rg' in test_mode: + srange_correction = isce3.core.LUT2d(srange_vec, + az_time_vec, + range_offset * ones, + method) + elif 'tec' == test_mode: + srange_correction = \ + tec_lut2d_from_json(unit_test_params.tec_json_path, + unit_test_params.center_freq, + unit_test_params.orbit, + test_case.radargrid, + isce3.core.LUT2d(), + unit_test_params.dem_path) + test_case.srange_correction = srange_correction + + az_time_correction = isce3.core.LUT2d() + if 'az' in test_mode: + az_time_correction = isce3.core.LUT2d(srange_vec, + az_time_vec, + azimuth_offset * ones, + method) + test_case.az_time_correction = az_time_correction + + # prepare input and output paths + test_case.input_path = \ + os.path.join(iscetest.data, f"geocodeslc/{axis}.slc") + + common_output_prefix = f'{axis}_{test_case.test_mode}' + test_case.output_path = \ + unit_test_params.out_dir / f'{common_output_prefix}.geo' + + # assign validation array based on axis + if axis == 'x': + test_case.validation_arr = unit_test_params.grid_lon + else: + test_case.validation_arr = unit_test_params.grid_lat + + yield test_case + + +def run_cuda_geocode(test_case, unit_test_params): # Set interp method interp_method = [isce3.core.DataInterpMethod.SINC] raster_dtype = [isce3.io.gdal.GDT_CFloat32] invalid_value = [np.nan] - # If block mode, i.e. suffix != '', have lines per block smaller than - # raster size. Otherwise set lines_per_block large enough to fit entire - # raster. - for xy, (suffix, lines_per_block) in itertools.product(['x', 'y'], - (['_blocked', 126], - ['', 1000])): - # Init CUDA geocode obj - cu_geocode = isce3.cuda.geocode.Geocode( - geogrid, rdr_geometry, lines_per_block) + # init CUDA geocode obj for given test case + cu_geocode = \ + isce3.cuda.geocode.Geocode(unit_test_params.geogrid, + test_case.rdr_geometry, + lines_per_block=test_case.lines_per_block) - # Put all input parameters in iterable list. geocode_rasters only - # works with lists of the same length. - output_raster = [isce3.io.Raster(f"{xy}{suffix}.geo", geogrid.width, - geogrid.length, 1, - gdal.GDT_CFloat32, "ENVI")] + # output and input as arrays to match binding signature + output_raster = [isce3.io.Raster(str(test_case.output_path), + unit_test_params.geogrid.width, + unit_test_params.geogrid.length, 1, + gdal.GDT_CFloat32, "ENVI")] + input_raster = [isce3.io.Raster(test_case.input_path)] - input_raster = [isce3.io.Raster(os.path.join(iscetest.data, - f"geocodeslc/{xy}.slc"))] + # geocode raster for given test case + cu_geocode.geocode_rasters(output_raster, + input_raster, + interp_method, + raster_dtype, + invalid_value, + unit_test_params.dem_raster, + native_doppler=unit_test_params.native_doppler, + az_time_correction=test_case.az_time_correction, + srange_correction=test_case.srange_correction) - cu_geocode.geocode_rasters(output_raster, input_raster, interp_method, - raster_dtype, invalid_value, dem_raster) + output_raster[0].set_geotransform(unit_test_params.geotrans) - output_raster[0].set_geotransform(geotrans) -def test_validate(): +def test_cuda_geocode(unit_test_params): + ''' + run CUDA geocode to make sure it does not crash + ''' + # run array mode for all test cases + for test_case in geocode_test_cases(unit_test_params): + run_cuda_geocode(test_case, unit_test_params) + + +def _get_raster_array_masked(raster_path, array_op=None): + # open raster as dataset, convert to angle as needed, and mask invalids + ds = gdal.Open(raster_path, gdal.GA_ReadOnly) + test_arr = ds.GetRasterBand(1).ReadAsArray() + if array_op is not None: + test_arr = array_op(test_arr) + + # mask with NaN since NaN is used to mark invalid pixels + test_arr = np.ma.masked_invalid(test_arr) + + return test_arr + + +def test_validate(unit_test_params): ''' validate test outputs ''' + # check values of geocoded rasters + for test_case in geocode_test_cases(unit_test_params): + # get masked array phase of complex test data and corresponding mask of + # pixels outside geocoded radar grid + test_phase_arr = \ + _get_raster_array_masked(str(test_case.output_path), np.angle) - # check values of 2 geocoded rasters - for xy, suffix in itertools.product(['x', 'y'], ['', '_blocked']): - # get phase of complex test data and mask zeros - test_raster = f"{xy + suffix}.geo" - ds = gdal.Open(test_raster, gdal.GA_ReadOnly) - test_arr = np.angle(ds.GetRasterBand(1).ReadAsArray()) - test_mask = test_arr == 0.0 - test_arr = np.ma.masked_array(test_arr, mask=test_mask) - ds = None - - # get geotransform from test data - geo_trans = isce3.io.Raster(test_raster).get_geotransform() - x0 = np.radians(geo_trans[0] + geo_trans[1] / 2.0) - dx = np.radians(geo_trans[1]) - y0 = np.radians(geo_trans[3] + geo_trans[5] / 2.0) - dy = np.radians(geo_trans[5]) - - # use geotransform to make lat/lon mesh - pixels, lines = test_arr.shape - meshx, meshy = np.meshgrid(np.arange(lines), np.arange(pixels)) - grid_lon = np.ma.masked_array(x0 + meshx * dx, mask=test_mask) - grid_lat = np.ma.masked_array(y0 + meshy * dy, mask=test_mask) - - # calculate and check error within bounds - if xy == 'x': - err = np.nanmax(test_arr - grid_lon) - else: - err = np.nanmax(test_arr - grid_lat) - assert(err < 1.0e-5), f'{test_raster} max error fail' + # compute maximum error test phase and validate + err_arr = np.abs(test_phase_arr - test_case.validation_arr) + + err = np.nanmax(err_arr) + + assert(err < 5.0e-7), f'{str(test_case.output_path)} max error fail' + + # count number of invalids (NaN) and validate + percent_invalid = \ + np.count_nonzero(np.ma.getmask(test_phase_arr)) \ + / test_phase_arr.size -if __name__ == '__main__': - test_cuda_geocode() - test_validate() + assert(percent_invalid < 0.6), \ + f'{str(test_case.output_path)} percent invalid fail' diff --git a/tests/python/packages/isce3/geocode/geocode_slc.py b/tests/python/packages/isce3/geocode/geocode_slc.py index 19d306dbf..c5fa804c5 100644 --- a/tests/python/packages/isce3/geocode/geocode_slc.py +++ b/tests/python/packages/isce3/geocode/geocode_slc.py @@ -532,7 +532,6 @@ def validate_slc_raster(unit_test_params, mode, raster_layer=1): test_phase_arr, test_mask = \ _get_raster_array_and_mask(output_path, raster_layer, np.angle) - # use geotransform to make lat/lon mesh ny, nx = test_phase_arr.shape meshx, meshy = np.meshgrid(np.arange(nx), np.arange(ny))