Question about routine radar_coords_to_cart in multi-radar gridding

[miaoneng]

Recently I reviewed the procedure of multi-radar gridding and find something might not be accurate in grid_mapper. If my understanding is correct, map_to_grid function in grid_mapper module iterate a list of radars, apply radar_coords_to_cart for each sweep to obtain a (x, y, z) for each gate. However, the radar_coords_to_cart function is a radar-station centric function. This means azimuth==0 indicates the true north thus it is along the Y-axis. It is OK in any single radar mapping but in any multi-radar mapping, as long as two and more radar (especially they're far away)

Consider an unreal corner case like this. If two radar are located near to north pole: Radar A is (88.9N, 90E), Radar B is (88.9N, 90W). If mapping together, it would be a wrong to use radar_coords_to_cart individually and directly mosaick them using calculated (x, y): actually y_a and y_b should on same line but current routine in grid_mapper probably produces two parallel lines for y_a and y_b but put x_a and x_b on a same line.

To correct this "rotating" issue, I got two ways now I could think out to solve this (if it is an issue). One way is plot multi-radar grids in a geographically projected map, not a lat-lon or unprojected planar map, then plus the geodesic convergence angle at the radar station to all azimuth recored; the second way is in radar_coords_to_cart, the finally result is converted back to lat-lon then all final results are plotted in lat-lon coordinates. The previous one needs the a calculation of geodesic convergence angle but I cannot figure out how to do this in proj.4 lib (unless we use a conic projection). The second way needs a lot of conversion from x-y to lat-lon along the great circle (might cause performance issue). Probably an another way would be completely change the way to calculate range gate using the equation given by Doviak and Zrnic in 1993. However, that could be more code to add.

I am not sure if it is considered in other procedure so the radar azimuth is actually corrected by geodesic convergence angle. If so, could you be kind to let me know the branch/file. I find this problem in my other research so I came back to find some references here but failed.

Thanks

[jjhelmus]

Thank you for reporting this @striges, this is most definitely a defect in the multi-radar gridding routines in Py-ART. Some of this issue was previously identified in issue #294.

There seems to be two issues with the way gridding of multiple radars is currently done in Py-ART:

1. Finding the Cartesian distances between two radar stations currently uses the corner_to_point function which performs a sinusoidal projections, see the issue comment by @gamaanderson. This is not appropriate for radars which are far apart. Issue #294 is primarily concerned with this and a solution was proposed but has yet to be implemented.
2. To find the Cartesian coordinate of each radar gate, a sum is taken between the Cartesian displacement of the radar to the grid origin (calculated using corner_to_point) and the Cartesian location of each gate from the radar origin (calculated using radar_coords_to_cart). As pointed out by @gamaanderson this is a connection. Although not exact, I though that the error introduced by this procedure would be minimal. From this issue I've realized that the problem in using a connection is not only that the sum introduces small errors, but also that the two coordinate systems may not be aligned and are rotated. This is the main concern of this issue.

I'm not certain as to the best solution to this rotation issue. Perhaps the only method is to calculate the latitude and longitude of all radar gates and then project these onto a single Cartesian grid? This was proposed in issue #294 but seemed unnecessarily computationally demanding.

Perhaps @gamaanderson or @deeplycloudy can provide some insight into this issue?

[gamaanderson]

You are correct, I have missed that.

A connection defines not only the distortion of a translation, but also the rotation and that can be really big. The problem is that Doviak and Zrnic give us formulas for a radar specific projection, to plot several radar we need to pass that to a single projection. If that one is geographic we just need x,y -> lat,lon; if not, since we don't have a x,y -> x',y', we need to do x,y -> lat,lon -> x',y'. Now it is just the computational problem of where and how we do those conversions, and in this question I don't think I have anything to say that @jjhelmus don't already know.

[miaoneng]

I explored a little more on geodesic convergence angle correction. If we use Lambert Conformal Conic, this can be calculate easily. For a radar station at (lat_r, lon_r) the angle correction is (lon_r - lon_0) * sin(lat_0), here lat_0, lon_0 are origin point(0,0) in the projection. Then we could apply this angle correction on all radar stations.

In terms of the geodesic distance between two radar stations, a simple way is to project two radars in to LCC then use Cartesian distance. Actually if we define standard parallels at 30N and 50N with original latitude at 40N, central meridian at 96W with NAD83 datum, the largest distance distortion will be around 2% in CONUS. It is more accurate than using a constant earth radius (R=6371km) to calculate geodesic distance.

Since LCC is projectable and inverse-projectable, it will be easy to inverse-project any range gate back to lat-lon. When plotting, we could plot using Basemap. Then only concern is, there are so many gates to project and inverse project. This could be a performance bottleneck.

[gamaanderson]

I think a important question is: do we want to implement projection specific formulas? I'm sure for any projection we choose there is a more or less simple formula for the rotation. More important is that we give the user enough control so he can bypass the bottleneck if he is willing to pay the price.

[miaoneng]

Here is a few additional comments after reading comments from @gamaanderson in issue #294. From my point of view, we should always avoid R=6371km at any time. Because pyproj is already available there, it would be not difficult to calculate geodesic distance using functions in pyproj/geodesic.py. In a single radar scenario, Azimuthal Equidistant is definitely a good choice if we put original point at radar station. For multi-radar case, currently I think we can only use a conformal conic projection because it preserve the shape of any great circle.

But I think @gamaanderson is idea is better: do a (x,y) -> (lon,lat), if py-art internally stores all rage gates using (lat, lon), gridding could be done in any error-free way.

[miaoneng]

@gamaanderson You point is also one of my concern: whether a "projection specification" thing is good in Py-ART?. I think the answer is no. But I mentioned LCC because I consider this is an easy way to make such correction and inverse-project range gates back to lat-lon for any further plotting or flatting. Another thought is: in the real application, the error is really minimal. So we can safely ignore this.

[gamaanderson]

Using R=6371km is the price one pays for making manual implementations, I'm also no fan of it. Keep in mind matplotlib and therefore pyproj is still a optional dependence of pyart.

Maybe I'm not getting your point. Are so suggesting to convert all to LCC, make translation and rotation there and then convert back? Considering pyproj available this is just as costly as converting to lat,lon (as long as I know), but in lat,lon there is not translation to be done so no rotation need to be corrected.

[miaoneng]

Sorry for the confusion. Let me give an example. If a radar is at (lat_r, lon_r), given a beam at azimuth alpha. Considering to a beam actually is a great circle, so we get all range gates on that beam presented by (lat, lon) using pyproj by calculating the geodesic distance with geodesic bearing of alpha. We can repeat this operations for all radars in a multi-radar case. But then if a distance-weighted interpolation is needed, we come to the problem that we need reproject (lat, lon) back to in order to get the line distance. I consider this procedure involves too many unnecessary project/inverse project calculations. So I think why not put all radars on LCC, make geodesic convergence corrections on each radar then gridding, interpolating. If a lat-lon grid is desired, we could inverse project them back. If another projection is desired, we project the mosaic to another coordinates. That simplifies calculations. Furthermore, I just feel that LCC is not a bad idea at all.

Is pyproj is a dependency of matplotlib? I don't think matplotlib requires pyproj because it is directly based on either proj4 or gdal C codebase (which one I don't remember)

[jjhelmus]

pyproj is included with the matplotlib Basemap toolkit. This is a separate install from matplotlib. Currently, Py-ART requires matplotlib to be installed but not Basemap. Basemap is quite a large install (~180 MB) as it includes a number of large data files. Therefore we would like to keep Basemap as an optional dependency, which means pyproj is an optional dependency.

The rough plan is to implement one map projection directly in Py-ART (azimuth equidistant seems like the best choice) so that multi-radar gridding is possible without pyproj/Basemap as long as the user accepts this projection. If Basemap/pyproj is available the user would be able to choose any number of map projections for the gridding.

I should not that it is possible to install pyproj independently of Basemap, but Py-ART does not support this yet. prproj by itself is a much smaller install (~1-2 MB).

[miaoneng]

Yes. The azimuth equidistant is a good choice. Simple to implement and small distortion on conformality. The convergence angle is simply the difference of longitudes.

[gamaanderson]

Thank you @jjhelmus for undoing my confusion with matplotlib instead of Basemap. An also the explanation about the size, I didn't knew that.

@striges So I think I did understand your suggestion. It is as valid as passing everything to latlon, just that LCC is preferred for high-latitudes.