Projection issue for VERY large mosiac

[scollis]

So this is one for the GIS crowd. Using Jonathan's new gridding engine I am doing a massive CONUS wide mesh.. But as you can see form the plot, unless we have radars off the eastern seaboard the mesh does not line up with the map..

Code: https://github.com/scollis/conus_mosiac/blob/master/scripts/conus_mosiac.py#L76

[tjlang]

Sorry, Scott, can't resist, but I fixed your problem in 7 lines of code. :) On a more serious note, yeah if the map projection is based off a single point you will get huge errors near the edges. I would need to dig up how MRMS fixed this, it's probably in the Zhang et al. (2011) paper or something.

[tjlang]

Yeah, reading the paper it sounds like they map each individual radar to its own cylindrical equidistant projection (0.01-deg resolution), then they merge all those individual grids together onto a national projection with the same resolution, and where radars overlap they use an exponential weighting function to merge the data. A lot of intervening QC and other corrections are done as well.

I did something similar, but far less fancy, for this pre-MRMS/NMQ paper: http://onlinelibrary.wiley.com/doi/10.1029/2009JA014500/full

So the basic gist is map individual radars locally, then merge after the fact.

[scollis]

Thanks @tjlang I admit significant ignorance when it comes to projection issues... It would be nice to resolve this so Py-ART's gridding paridigm and map plotting works.. A fundamental question is are the grid points in the right place? Is the issue in the calculation of the displacement of the radars from the origin of the grid?

[jjhelmus]

My guess is that the issue is the projection, or rather lack of a projection, used in calculating displacement of the radar location from the grid origin. This is done by the corner_to_point function which I believe is not valid over large distances. I'll investigate on Monday.

[scollis]

I concur. Basically we need a better function for working out the x-y displacement given latitudinal and longitudinal differences.. Perhaps GDAL

[tjlang]

Yeah for 150-200 km distance from the radar you should be good with whatever, but beyond that you have to be a lot more rigorous. I remember working out that math using fairly simple lat/lon-Cartesian transforms a while back.

[mvanlierwalq]

Isn't part of the problem that x/y distance depends on whether you take your lat or lon step first? If you take your steps infinitesimally, then.... you're basically calculating sin and cos of the forward azimuth and the distance of the great circle between the points? (I don't know... I'm asking). Is x/y displacement really a well posed question on a sphere?

[gamaanderson]

This can be easily solved using the pyProj, is just a question of using forward transformation, the oposite of add_2d_latlon_axis. I'm increasingly believing basemap should become a mandatory dependence.

[jjhelmus]

I've labeled this issue and will verify that it is caused by the poor use of projections and try to come up with a fix at a later time. If anyone else wants debug this further before then, feel free.

[gamaanderson]

I just realized this today, there is one more issue that you may want to take a look while fixing this. Some projection assumption are inconsistent, I will resume everything:

• line 466 uses function corner_to_point that uses a local approximation of the Sinusoidal projection, although I don't think this was a conscious decision
• line 473 uses function radar_coords_to_cart which assumes Azimuthal equidistant projection
• lines 479 and 480 by summing are doing a euclidean translation and therefore using a flat Earth model (FYI this is not choosing a projection, but actually a connection #fancyMathematics)
• In GridMapDisplay the map is plotted in Miller Cylindrical Projection, but in line 65 the data is considered to be in Lambert conformal conic projection (LCC)

Point 2 is probably your goal so this is correct, point 1 is the bug that created this issue and point 3 is limited to the radar range so its errors are smaller. Point 1 can be solved implementing Azimuthal equidistant transformation, its a big formula, but I have it somewhere. Points 3 would be probably better use pyProj, but anyway it will have some significant lost of time for little or non improvement in the results.

[nguy]

I'm by no means an expert, but what about the GDAL library instead of basemap like wradlib uses?

[scollis]

That is one thing we are considering.. However I am leaning much more to Proj as it packages with basemap. I will get to this one soon..

On 7/23/15 2:10 PM, Nick wrote:

I'm by no means an expert, but what about the GDAL library instead of basemap like wradlib https://bitbucket.org/wradlib/wradlib/src/6c2a093c8a04c011f2c2efb92f248a2f0e9de6d0/wradlib/georef.py?at=default uses?

— Reply to this email directly or view it on GitHub https://github.com/ARM-DOE/pyart/issues/294#issuecomment-124214370.

I ride for Parkinsons research http://www.events.org/sponsorship.aspx?id=51573

[jjhelmus]

This issue isn't with basemap vs GDAL vs Captopy so switching away from basemap does not help solve the problem. The issue is that we are not using any a uniform (or any) projection when gridding nor plotting grids as @gamaanderson's comment explains. Basemap can do projections and is already used in Py-ART so my vote would be for that.

[gamaanderson]

I hate to be the one to point one more problem, but plot_latitude_slice and plot_longitude_slice don't do what they say they do. They plot a slice of the data, but as the projection transformation don't preserve straight lines, this is not a slice in a given latitude (as plot_crosshairs would suggest).

To prove this I made an empty grid with the same size as @scollis in the beginning of this issue and them put data into a slice and plot.

them I superposed that into the original image and (manually) circle the regions in the plot and the slice to show they agree

here is the script

If one were to ask me, I would prefer to maintain the function working as they do now, just change they presentation to make clear that fact, also having as parameters the x and y index instead of lat lon ( just as plot_grid get the z index). To make true "straight" slices we should start supporting grids in latlon (or any other cylindrical) projection in issue #285

[scollis]

Very true.. As often pointed out the care in Py-ART put much care into radar classes.. less into grid classes.. we have some work to do here but fortunately we may have some funding soon to tackle this in a holistic manner..

Thanks for pointing this out and feel free to submit an issue.

On 7/25/15 2:48 PM, Anderson wrote:

I hate to be the one to point one more problem, but |plot_latitude_slice| and |plot_longitude_slice| don't do what they say they do. They plot a slice of the data https://github.com/ARM-DOE/pyart/blob/master/pyart/graph/gridmapdisplay.py#L238, but as the projection transformation don't preserve straight lines, this is not a slice in a given latitude (as |plot_crosshairs| would suggest).

To prove this I made an empty grid with the same size as @scollis https://github.com/scollis in the beginning of this issue and them put data into a slice and plot. y-slice https://raw.githubusercontent.com/gamaanderson/scripts/master/y_slice.png

them I superposed that into the original image and (manually) circle the regions in the plot and the slice to show they agree scollis https://raw.githubusercontent.com/gamaanderson/scripts/master/scollis.png

here is the script https://github.com/gamaanderson/scripts/blob/master/script.py

If one were to ask me, I would prefer to maintain the function working as they do now, just change they presentation to make clear that fact, also having as parameters the x and y index instead of lat lon ( just as |plot_grid| get the z index). To make true "straight" slices we should start supporting grids in latlon (or any other cylindrical) projection in issue #285 https://github.com/ARM-DOE/pyart/issues/285

— Reply to this email directly or view it on GitHub https://github.com/ARM-DOE/pyart/issues/294#issuecomment-124880377.

I ride for Parkinsons research http://www.events.org/sponsorship.aspx?id=51573

[deeplycloudy]

I have code that puts radar data into arbitrary map projections using pyproj here: https://github.com/deeplycloudy/lmatools/blob/master/coordinateSystems.py#L127 The class has methods to convert (r, az, el) to and from earth-centered, earth-fixed (ECEF) cartesian WGS84 coordinates. It includes the standard 4/3R refraction correction.

From ECEF coordinates, you can use the map projection class at https://github.com/deeplycloudy/lmatools/blob/master/coordinateSystems.py#L76 to use any of the proj4 map projections, or get lat/lon/alt of any gate using https://github.com/deeplycloudy/lmatools/blob/master/coordinateSystems.py#L55

While this shows the method for getting everything in a common coordinate system, I'm not sure how you'd want this integrated into Py-ART. If you want to point the way within Py-ART, I'd be glad to help with a PR.

[gamaanderson]

Just a comment, the code above uses pyproj.Geod, that is, geodesics; This is probably the most "correct" one can gets, but it is also much slower than pyproj.Proj. However the difference, if any, from Geod to Proj (projection='aeqd') is minimal and only doe to earth eccentricity, in a spherical earth they are egual.

If someone can confirm this to me, I believe in an elliptical earth Proj will work over the intersection of the ellipsoid with planes passing in the center, but those are not the geodesics.

[deeplycloudy]

Regarding the timing for the geodesics, I did a %%timeit timing test on 10.8 M data points using proj.4 https://gist.github.com/deeplycloudy/6390a32b0b15d70fad7d Note that the above depends on a new method I added to lmatools/coordinateSystems.py

(r,az,el) to (lon,lat,alt) (via geodesics): 5.5 s (lon, lat, alt) to aeqd projection: 3.4 s (lon, lat, alt) to ECEF XYZ: 1.7 s

So a map projection in proj.4 takes about 60% of the time of the more complicated geodesic calculation. And the map projection takes twice the time required to calculate an earth centered, earth fixed cartesian XYZ.

For the purposes of gridding, especially at large distances with significant earth curvature, I think it's important to be careful, and I'd tolerate that kind of speed penalty for exactness. I'm particularly tolerant of slightly slower projections during gridding, when the gridding process itself will take much longer than 5 seconds.

In perusing the issue tracker, I noticed that fixing this issue is also related to #247.

[gamaanderson]

Well, I'm not totally sure those are the relevant test, but you got the idea. However those error are really small, limited to radar range. The biggest error we got are:

1. first and at most: incorrect/inconsistent calculation of radar center
2. not that much, but potentiality perceivable, flat earth assumption while translating the gates

[deeplycloudy]

I think the flat earth assumption becomes a big problem when gridding multiple radars in three dimensions - it's no longer limited to radar range.

One thing @jjhelmus, @scollis, and @tjlang and I talked about at the radar conference was implementing a favorite coordinate transform (sounds like az. equidistant is a good choice?) within pyart to keep dependencies to a minimum, and then providing the proj.4 wrappers for anyone who wants to be more careful. That might be a good way to satisfy diverse tolerances for speed vs. exactness at the same time.

[gamaanderson]

Exactly, multiple radars!

That is a solution, we implement spherical earth aeqd -> latlon and latlon -> aeqd by hand, everything else needs proj4. We just need to decide how to implement that into grid. I would like to see Grid provide this transformation, for instance in Grid.xy_to_latlon(x, y) and Grid.latlon_to_xy(lat, lon), so other functions don't have to take care if it is proj4 or not.

We also have to decide how to save the projection information, I would propose as attributes in a NetCDF like variable. Also, what convention we use for those attributes? proj4 arguments or netcdf-java. Proj4 is part of a big community and used by several softwates, while netcdf-java implement the transformation by it self (I checked in the code) and follow a book as reference, not other softwares/libs. As much as I want to support Unidata format, I just don't trust their softwares, so I would use proj4 arguments and make a Proj4->netcdf-java attributes conversion.

[jjhelmus]

I have been silently following this thread but should probably speak up now since I'm hoping to fix this soon. From the discussion my understanding is that the following would solve the majority of the issue. My understanding of the topic is limited and I'm sure @gamaanderson and @deeplycloudy have more knowledge and experience on the matter, so please let me know if I am incorrect on any of these points.

1. Replace the corner_to_point function with one that uses a projection which can be specified by the user. The default projection (az. equidistant?) will use code built into Py-ART to perform the projection when basemap/pyproj is not available.
2. Leave the radar_coords_to_cart function as is. The transforms here are taken from Doviak and Zrnic and I cannot figure out how they can be modified to allow for a selectable projection.
3. In the gridding routines use corner_to_point to calculate the displacement in x, y and z from the grid center to the location of each radar. In the case of a single radar with grid centered at that radar no calculation is requires.
4. Also in the gridding routines, calculate the x, y and z location from the radar to each gate using radar_coords_to_cart. Add these to the displacement of the radar to the center from the previous step. This is a connection as mentioned in am earlier comment by @gamaanderson.
5. The issue with the GridMapDisplay class have been address in already merged Pull Requests (largely #346, thanks @gamaanderson).
6. Store this projection information in some location of the Grid class and use this when appropiate in the GridMapDisplay.

There may be small errors introduced by points 3 and 4 but these will be limited to the radar range. To remove these errors the latitude and longitude for each gate would need to be calculated from the x, y, and z distances from the radar center calculated using radar_coords_to_cart. Then the x, y and z distances from these latitudes and longitudes to the grid center would need to be calculated using pyProj or something similar.

This procedure would require two transformation: xyz_radar -> latlon -> xyz_grid on all gates. For large mosiacs this could be compuationally expensive when millions of gates are involved. I'm guessing the small errors introduced by 3 and 4 would be acceptable considering the extra complexity and computational requirements that the proper procedure would require.

[gamaanderson]

You are absolutely right @jjhelmus.

Just in point 2 we should consider the Cartesian coordinates in az.equidistant, so if the user choose a different projection, with potential big deformation, correction would be needed. This could be done with a option of performing or not the xyz_radar -> latlon -> xyz_grid transformation.

There are two open questions (that I know you are award): how the user select a different projection? How the handle the internal implementation vs proj4 in a transparent but automatized way?

Also, this is not really necessary, but one could also add projection information for radar, this would (for instance) change the output of radar_coords_to_cart. This is really a choose, there are good reasons for and again implementing it.

[gamaanderson]

just for reference here are manual implementations of the az.equidistant transformations

[jjhelmus]

Good points. Let me think on the topic this weekend. I'll try to put together a PR to implement some of these ideas next week. I think having some code to examine and critique will help guide some of the discussion on the selection and representation of the projections.