I’ve been wrestling with x2sys to examine ICESat-2 cross-overs. Following some examples found on the forum, I have successfully produced cross-overs and have solved for their adjustments.
With ICESat-2 data, it can involve huge matrices (depending on region size), especially when each of the six beams are treated as independent (uncorrelated) measures. When two ICESat-2 tracks cross, a total of 36 cross-overs are formed by the individual beams.
What I haven’t found is a way to report/propagate the least-squares solution errors back through the system to provide an error estimate; i.e., how well determined are the adjusted tracks? I’ve tried forming my own system, based on what was observed and what was computed, but the inversion is unstable.
Bottom line question: is there an additional x2sys command that can form the inverse(A’PA) and scale it by the variance of unit weight?
It has taken quite some effort to work through what I was hoping to do.
After reading Paul’s 2010 paper on x2sys, I came to understand his use of Lagrange Multipliers (eq. 9 of said paper) for matrix stability and insuring independence between clusters. Initially, I had reconstructed the normal equation matrix (N) without components of the matrix for the multipliers and, most of the time, couldn’t obtain an inverse.
Once the extra rows and columns are added to the matrix, to untangle the individual clusters, it inverts just fine. Then, I can grab the sqrt(diagonals) of N and scale by the sigma of unit weight (i.e., sigma of the residuals) to obtain an estimate of how well-determined the cross-over adjusted heights are, on each track.
What made this hard is that I am treating each ICESat-2 beam (of six beams) as individual “tracks.” I had to remove any cross-overs occurring on nearly-aligned tracks (using differences of the headings) to create the COE list, then feed that into x2sys_solve. Kind of ugly, but I suspect that when x2sys was written, Paul wasn’t thinking about the prospect of lots and lots of beams and crossings. Much better when applied to radar altimetry. A new option in x2sys_list could do something like that, to insure that all cross-overs are ascending/descending pairs (or did I miss that option?).
This was what I was after, and I’m guessing that the way I phrased the original post was beyond comprehension (perhaps even wrt my own!).
Paul, I’d like to suggest that one more x2sys command (x2sys_error?) be considered that can provide the formal errors associated with each track, after the adjustment. That would save all the rigmarole I devised (in scripts and in Fortran) to accomplish this.
Seems like we could add a new option to x2sys_cross that would skip crossings where the azimuths at the cross points differ by less than cutoff, like -S excludes crossings when the ship was moving too slowly to be mapping anything. Maybe -Ecutoff would work. The code has to check for line orientation, not azimuth, and then compute difference in orientation and if less than cutoff we skip it. I don’t think there is a good default cutoff value so it needs t one set. John, I am busy with other non-GMT things so do you mind turning this into a feature request on GitHub.
As for the error calculations. Should it not be done by x2sys_solve instead?
I’ll look for it when it is released in mainstream GMT.
Yes, I agree that error calculations seem most appropriately performed within x2sys_solve. Basically you need that N-matrix, including the rows/cols with Lagrange Multipliers, to be inverted, then scale the diagonal terms by the variance of unit weight (i.e., sum of the squares of the residuals, divided by the DoF), followed by a square root.
In terms of implementation, the development team will have to decide how to make it jibe with existing options. It doesn’t seem to be in any way related to the -E option. Whatever way it is done, the output needs to go somewhere, in addition to the general output of the command.