The following is a ‘snaphot’ taken May 2012. In particular, the numerical stats used May 10, 2012 as a sample, but the values are similar from day to day. The features discussed here are modified over time and later updates will be needed.
The Q2 system obtains HADS gauges from NCEP via a ftp site that allows anonymous public login at ftpprd.ncep.noaa.gov. Upon entering the site, a change of directory is made to pub/data1/nccf/com/ingest/prod. At this point, there will be another layer of daily sub-directories with the name pattern ‘shf.YYYYMMDD’. Typically, these directories exist for the past 30 days. The ‘shf’ in the directory name refers to the data being in the SHEF format.
Each daily sub-directory contains a handful of files, but the easiest to track is the one named ‘hrly.prcp.day.MMDD’. Note that the time stamp does not contain the year. Since this archive only goes back one month, there is no confusion about what is the corresponding year. This file gets written in 3 hr intervals a total of 10 times. Each time the file gets larger as more data is appended – while data from former version is always still present. The file is first written about 03:15 UTC on the date in question. The initial file runs about 1.7 Mb, and contains 1 hr gauge accumulations for times ending 00:00 to about 03:08 UTC. The gauge reports are always 1 hr accumulations, but the ending time can be anytime during the day. Most of them (currently ~80%) are reports ending on the hour, so that first file write has many reports for hours ending 00:00, 01:00, 02:00, and 03:00 UTC. But every minute in between has some reports and the most recent report in the 03:15 UTC file write will only be 5 -10 minutes old.
This file is over-written 3 hours later at about 06:15 UTC with a new version. This file is about twice as big (3.5 Mb). It will contain all the previous 00:00-03:08 UTC reports contained in the first version, plus a significant number of new reports. Most of these are for the most recent three hours (reports ending 04:00 – 06:08), but there will also be new reports from earlier in the day that first show up here. This pattern continues as the file is over-written 8 more times at three hour intervals; at 09:15, 12:15, 15:15, 18:15, 21:15 UTC that same day and 00:15 and 03:15 and 06:15 the following day. By the end the file is up to about 14 Mb.
The following figure shows the pattern of new reports. For simplicity, only the gauges that report on the top of the hour are included here. This is about 80% of all the gauge reports in the file. Again this is a snapshot from 5/10/2012. Each new run contains mostly new reports for the last 3-4 hours, but there is also a small number of older reports (up to 18 hrs old). The numbers in the main table in black are the total number of gauge reports in the file up to that point, while the number in red are just the count of those that are new for this file update.
In the blog post the figure text is fairly small since the image is sized to fit the region, but your browser should provide a simple method for viewing a full-size version of the image. For instance, in Google Chrome, right click on the image and select ‘Open image in new tab’ for a much better view.
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.
Spherical earth with radius = 6371 km (this is the IUGG recommended average earth radius obtained via a weighting of 2/3 R
and 1/3 R
Ellipsoid: Inverse of ‘Puissant’s Coast and Geodetic Survey’ method. Uses two constant parameters (semi-major axis = 6378.14 km and eccentricity = 0.081894)
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.
An Introduction to Geodesy
, by Ewing and Mitchell, American Elsevier Pub Co, New York,1970]
A classic iterative ellipsoid calculation by Vincenty can be used to get arbitrarily high accuracy. See reference from our parent organization:
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.
The Q2 blog is now online. See the ‘About’ section for a more complete introduction and mission statement. The format will change for awhile and more menu options will appear as things are posted and more organization is required (and the administrators figure out how all this works).