Closed parosen closed 1 month ago
Thank you so much, @parosen , for reporting this issue and for providing all the information we need to debug it.
Indeed, the error is in ISCE3 in the code that creates the zeroDopplerTime
vector for the geolocationGrid
cubes here.
In short, the problem lies with the end time point for np.linspace()
function. We need to remove the subtraction of the d_az
term from az_f
. I already issued an internal pull request (PR) to fix this.
Detailed explanation:
The radar grid convention employed in ISCE3 follows the "pixel center convention", which means that the starting range and az. time refer to the center of the first pixel (first range sample of the first az. line). The zeroDopplerTime
vector also uses the pixel center conventions. The first azimuth point az_0
can be directly retrieved from radar_grid.sensing_start
and the last azimuth point az_f
can be obtained by multiplying the az. pixel spacing (PRI) d_az
by the number of az. lines minus one radar_grid.length - 1
. The "minus one" is due to the pixel center convention. Note that if there was a single pixel in the radar grid, the position of the first and last pixel would be the same (i.e., the center of that single pixel).
When creating the zeroDopplerTime
vector, the current code is calling np.linspace(start, stop, length, dtype=dtype)
as:
az_vect = np.linspace(az_0, az_f - d_az, radar_grid.length, dtype=np.float64)
The end time should be simply az_f
, i.e., we should remove the term d_az
from the end point az_f - d_az
:
az_vect = np.linspace(az_0, az_f, radar_grid.length, dtype=np.float64)
The bug causes an error of -d_az
in the last point, which represents an error of - i *
d_az/ (N - 1)
for each az. time point i
in N
az. lines. Notice that if we have 2 points, the error in the second point (end point) is d_az / 1
(not d_az / 2
). If we have 3 points, the error in the second point (middle point) is d_az / 2
(i.e., not d_az / 3
).
So, the az. time spacing error d_az_current_code
in the current (wrong) code is given by:
d_az_current_code = d_az - d_az / (N - 1) = d_az * (N - 1 - 1) / (N - 1) = d_az * (N - 2) / (N - 1)
If we only have the "wrong time spacing" d_az_current_code
and we want to compute the "correct time spacing" d_az
, we have:
d_az = d_az_current_code * (N-1) / (N-2) = (d_az_current_code - d_az_current_code) + d_az_current_code * (N-1) / (N-2) = d_az_current_code + d_az_current_code * { [N-1 - (N-2)] / (N-2) } = d_az_current_code + d_az_current_code * [1 / (N-2)]
In other words, if we have the wrong time spacing d_az_current_code
, we have to add d_az_current_code * [1 / (N-2)]
to get the "correct time spacing d_az
, which is in agreement with your findings above.
@fastice FYI, please see @gshiroma's response above.
In testing personal code to convert Range/Doppler/height to lat/lon on a few boundary points from the grid:
When the time is zero, using the statevectors, my code gets the near and far range perfectly (see images below).
But as the time increases, the error grows linearly. I can’t see how my code could be the cause because it only uses the time to interpolate the state vectors. (if it were range I would think it was me).
Then I noticed the shift is equivalent to about a time step. So I empirically determined a correction as dT/(N - 2), where N is the number of time steps, which fixed the error perfectly (see second figure with the only difference being I commented out the zeroing of the correction).
I don’t know why N-2 works, I would more expect it to be N-1 (or even N).
But my guess is that somewhere in the code that computes the geolocation grid, the zeroDoppler time step is being computed incorrectly (probably a something like a divide by N instead of N-1, or vice versa).
I can’t entirely rule out an error on my part, but I am 95% sure it's the NISAR code (plus this would probably explain the shortfall in the grid because the time step is a bit too small).