search

OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.74k stars 788 forks source link

Azimuthal Equidistant projection is not accurate - Suggestion for fixing #246

Closed proj4-bot closed 9 years ago

proj4-bot commented 9 years ago

Reported by bbauerma on 24 Sep 2014 14:41 UTC Dear Proj.4 team!

We used the GDAL/OGR libraries for the geotiff/shapefile transformation from ellipsoidal WGS84 to planar coordinates, using the Azimuthal Equidistant. We discovered significant location inaccuracies in the peripheries (distance to the centre > 1000km, errors >> 1km). We traced the errors back to the approximating formulae for azimuths and geodesics in the proj4s aeqd.c file, where the formulae are as in Snyder 1987, page 199 + 201.

After some research, we found the work about geodesics of C.F. Karney and we already successfully integrated his Python-library of the GeographicLib into our code for the precise coordinate transformation. More precisely, we used for the forward projection his Geodesic.WGS84.Inverse(lat0, lon0, lat, lon) for the calculation of the geodesic distance c and its azimuth Az. For the inverse projection, we used his Geodesic.WGS84.Direct(lat0, lon0, azi, dist).

Consequently, we changed the code in our local aeqd.c so that it uses the geodesic-functions from the GeographicLib. When using this modified aeqd.c via GDAL/OGR, we match perfectly with the already accuracy-tested results from ArcGIS. Im wondering if the algorithm for the geodesic forward and inverse problem shouldnt be integrated also in the proj4-libraries for projections (Azimuthal Equidistant, Cassini-Soldner, ?) so that GDAL/OGR-community can make use of it too. We are happy to provide our code.

Thanks in advance, Bernhard

Migrated-From: https://trac.osgeo.org/proj/ticket/246

proj4-bot commented 9 years ago

Comment by rouault on 24 Sep 2014 18:53 UTC Attaching your modified code as a patch would be a starting point.

proj4-bot commented 9 years ago

Comment by bbauerma on 2 Dec 2014 11:31 UTC Dear proj.4 team!

Is there any update on this ticket? Do you need further information from our side?

All the best, Bernhard

proj4-bot commented 9 years ago

Comment by rouault on 2 Dec 2014 11:42 UTC Proj 4.9.0 has now geodesic API : http://trac.osgeo.org/proj/browser/trunk/proj/src/geodesic.h

Your patch should be updated against proj trunk to use the new geodesic API that has been integrated. Attaching your changes as a real patch against the subversion head would be easier for review too.

Other point: I didn't see new tests that show the improvements. This would typically involve editing nad/testvarious and nad/tv_out.dist

proj4-bot commented 9 years ago

Comment by bbauerma on 2 Dec 2014 16:59 UTC Hey there! Thanks for the quick response.

I have to say that I do not understand much of software design and so on - I'm just a stupid scientist ;-) I was hoping you and your team can fix this inaccuracy and implement it into proj4 so that others can also profit and we ourselves do not have to apply our patch. So I want to highlight what was really changed by us:

The azimuthal_equidistant_enhancements.7z contains all edits we did so that it worked on our machines. Essentially, the only real changes are in PJ_aeqd.c: Lines 70-119 Lines 172-216 (the lines as in our PJ_aeqd.c) In there we are using the geodesic.c from Karney.

I hope it is clearer now to you. All the best, Bernhard

proj4-bot commented 9 years ago

Comment by qulogic on 21 Feb 2015 09:45 UTC Hi!

I've taken the liberty of going through Bernhard's files, and applying them to SVN. I then cleaned up the indentation, fixed a couple minor bugs (global g -> per-instance g), and optimized it slightly.

I have attached the change as a patch against SVN 2616. I have not added any tests yet, because a) I don't yet know how to run them, and b) I don't know how to get the ''correct'' results.

proj4-bot commented 9 years ago

Comment by qulogic on 21 Feb 2015 09:46 UTC PS, since this patch also changes the inverse transform, it may also fix #211.

proj4-bot commented 9 years ago

Comment by qulogic on 22 Feb 2015 01:39 UTC I added a test based on #211. You can see that the forward calculation result is slightly different from before. But this is not unexpected since the Snyder approximation is only good to about 800km. Additionally, taking that result and putting it into the inverse calculation ''does'' return the original location. So that should confirm that #211 is fixed, at least.

I don't really have a way to provide a completely arbitrary test as I have no other way of calculating this projection.

(Sorry, I forgot to tick the "replace attachment" option, so it's uploaded twice, but they're the same.)

proj4-bot commented 9 years ago

Comment by bbauerma on 26 Feb 2015 12:09 UTC Thank you very much for the fixing!

We tested the algorithm for some random points against a) methods from ArcPy (ArcGIS) b) "hand-made" mathematics according to Karney.

We concluded that in each method a), b), and this algorithm, the inverse- and forward projection conclude without residue.

I can suggest as test:

of test coordinates --> inverse proj. --> forward proj --> coordinates

of test coordinates - coordinates = residue =? 0.000000

How is the procedure to make this fix effective to the latest PROJ4 and GDAL/OGR release?

Thanks again, Bernhard

proj4-bot commented 9 years ago

Comment by qulogic on 1 Mar 2015 06:58 UTC I did make a test of forward->inverse to see if it worked out based on #211, but there's only one.

I guess the process for getting this patch in a release depends on rouault or warmerdam. On a related note, the to-be-released cartopy will include support for this projection, but due to this issue, full global plots end up very distorted, so it would be nice if the fix were in some release that could be pointed to.

proj4-bot commented 9 years ago

Comment by rouault on 1 Mar 2015 08:24 UTC I don't know if Charles has an opinion on the new maths ?

proj4-bot commented 9 years ago

Comment by karney on 2 Mar 2015 14:29 UTC I haven't checked this enhancement at all. However, the basic methodology is sound: rely on proj's new-found ability to do geodesic calculations to define the azimuthal equidistant projection instead of relying on an approximate solution.

This is something to take up after the release of 4.9.1. I will look at it then (and bug me if I don't).

proj4-bot commented 9 years ago

Comment by bbauerma on 24 Apr 2015 11:49 UTC Sorry for inquiring again, but will this enhancement be included in the next proj.4 and/or GDAL release? I'm especially looking towards the GDAL 2.0 release here.

Can I support you guys? What kind of test is required? Thx!

proj4-bot commented 9 years ago

Comment by karney on 24 Apr 2015 11:52 UTC This is something I'll tackle once Howard Butler switches proj.4 to github. (Howard, is this move still on? I lost track of where this discussion ended up.)

proj4-bot commented 9 years ago

Attachment added by bbauerma on 25 Sep 2014 15:19 UTC https://trac.osgeo.org/proj/attachment/ticket/246/azimuthal_equidistant_enhancements.7z

proj4-bot commented 9 years ago

Attachment added by qulogic on 22 Feb 2015 01:35 UTC https://trac.osgeo.org/proj/attachment/ticket/246/proj-246-aeqd-fix.2.diff

proj4-bot commented 9 years ago

Attachment added by qulogic on 22 Feb 2015 01:35 UTC https://trac.osgeo.org/proj/attachment/ticket/246/proj-246-aeqd-fix.diff

QuLogic commented 9 years ago

This one's been fixed by #275/#281.

QuLogic commented 9 years ago

Ping?

rouault commented 9 years ago

@QuLogic What is this ping about : to just close the ticket ?

QuLogic commented 9 years ago

Yep.