The role of the model of the earth in gauge and radar bin co-location
Brian Kaney, Feb 2012
Given a gauge at a specific latitude and longitude, disagreement has been noted for the radial and gate that should be used to compare a single radar QPE product with that gauge. In particular, values from the real-time and case study QVS, an in-house comparison tool used by the ROC, the AWIPS system and an archive displayer from NCDC have all been considered.
Inverse Geodetic Problem
This problem has been much studied and is a classic in the field of geodesy. Given two arbitrary points on the globe referenced by latitude and longitude, how does one calculate the range and azimuth angle from one point to the other? This is known as the ‘inverse geodetic problem’. If approximating the earth as a sphere is sufficiently accurate for a given purpose, then an exact solution is available from standard trigonometry. The problem becomes worthy of a name only when one treats the earth more accurately as an ellipsoid. For two latitude and longitude referenced points on an ellipsoid there is no closed form solution to the problem. In other words, one cannot write down a single formula in standard functions that solves the problem. There are many approximation schemes and iterative procedures to get answers with different levels of accuracy.
Four different models for the earth’s shape were considered to at least some degree. One is spherical and three are ellipsoidal. They are listed below and will hereafter be referred to by the numbers.
1.) Spherical earth with radius = 6371 km (this is the IUGG recommended average earth radius obtained via a weighting of 2/3 Requatorial and 1/3 Rpolar).
2.) Ellipsoid: Inverse of ‘Puissant’s Coast and Geodetic Survey’ method. Uses two constant parameters (semi-major axis = 6378.14 km and eccentricity = 0.081894)
3.) Code from the RPG library. This was code sent to HMRG from the ROC back in May or June 2011. This code looks to be an ellipsoidal model and contains two constant parameters: 6380 km and 135 km. An attempt by the HMRG to obtain from the ROC the meaning of the two constants or an overall reference for the technique was not successful. The code looks similar to some techniques HMRG found in the literature, but there was never an exact match.
[Ref: An Introduction to Geodesy, by Ewing and Mitchell, American Elsevier Pub Co, New York,1970]
4.) A classic iterative ellipsoid calculation by Vincenty can be used to get arbitrarily high accuracy. See reference from our parent organization: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf.
Some sources regard this as the current state of the art.
At the time of this writing, the real-time QVS has always used method 1. In the summer of 2011 some dual pol assessment case studies were developed and methods 1, 2, and 3 were all used at different times. The HMRG team was trying to match gauge/QPE comparison statistics that the ROC already had for these events from their own in-house tool. Method 1 was tried first and discrepancies were found. At that point a simple ellipsoid model was sought, namely method 2, and that gave somewhat better agreement overall, but there were still errors. The improvement was not uniform with a significant number of gauges being moved farther from agreement. In both cases, typical errors would be 0.1-0.6 degrees in azimuth and 100 m – 500 m in range.
Then the RPG library code was obtained from the ROC. That method was copied and yielded essentially perfect agreement when ran (apart from rounding differences in the last decimal place, perhaps). Of course, this was not a test of the methods true accuracy as it was just a comparison of the method with itself. In other words, the exercise was just a test that the HMRG could reliably copy formulas into new code.
Since we have no references for method 3, the case studies were set back to using method 1. The HMRG feels it is important to be able to document and reference the techniques and scientific decisions that go into its systems. The fact that other displays, AWIPS or NCDC give errors of a similar size (0.4 degree azimuth and 300 m range) is worth noting, but without access to their source code there is little practical value beyond that. At some point a better earth model might be applied within HMRG, but that path is not clear at the present time. This is also intertwined with other issues. For instance, this entire process still assumes the well known standard solution to the radar propagation equations for an ideal simplified atmosphere. This standard model has radar beams propagating in a large circular arc of 4/3 times the radius of the earth. As the beam propagates it is bent earthward due to refraction as it climbs into thinner and cooler layers of the atmosphere. But this model assumes a spherical earth. Does it even make sense to use this beam model and then suddenly ‘shift’ the ground underneath to an ellipsoid? The beam is going to bend downward toward an ellipsoidal surface and not a spherical surface.
Method 4 was never coded up as it is significantly more complex than methods 2 or 3 and the concern becomes the computational overhead on the QVS to be able to generate images of polar products on the fly.
Comparison of methods
The coding of methods 1, 2 and 3 allowed for some comparisons of the techniques to be made. A map was drawn of the 230 km range region around some sample radar. If the map is in a Mercator projection then each point on the image corresponds to a point of latitude and longitude in a simple Cartesian fashion. For each point, the azimuth and range to that point were calculated via the various methods being considered. The results were then subtracted and color-coded and shown in the following figure.
In the azimuth, both methods have discrepancies compared to spherical that are near 0 degrees along a north south line through the radar. The extremes are about +0.7 degrees in the far west and -0.7 degrees in the far east. Method 3 and method 2 are not shown in direct comparison to each other, but looking at their similar behavior with method 1; they must be quite close to each other.
In the range, method 3 is always higher than method 1. The color scale used makes it difficult to tell how high this difference gets, but it is about 250 m halfway from the radar to the edge, so likely 500 m or more near the outer edge. Method 2 agrees best with a spherical earth along a NE-SW and a NW-SE line. The error is positive east and west of those lines and negative north and south of those lines. Again the extremes cannot be determined, but would appear to be at least +500 m and -500 m.
In the range, methods 2 and 3 are quite different in pattern. A difference between them would be especially large far north and far south of the radar – perhaps 1000 m or more. Although not a well substantiated assertion, the gate difference for method 3 does look suspicious to this author. One feature of an ellipsoidal earth is that the local curvature along the ground is slightly different along a north/south line as opposed to an east/west line. So if one ‘splits the difference’ and uses an average, one would expect the best agreement to be along a NW-SE line and a NE-SW line where the average matches the ground best and then progressively large (but opposite in sign) errors as one rotated to more of an east/west line or a north/south line. That is exactly the pattern seen in method 2, which is the one ellipsoidal method considered with references.