File Exchange

image thumbnail

geographiclib

version 1.50 (115 KB) by Charles Karney
MATLAB implementations of a subset of the C++ library, GeographicLib

25 Downloads

Updated 19 Sep 2019

View License

GeographicLib toolbox
Version 1.50 2019-09-19

This toolbox provides native MATLAB implementations of a subset of the
C++ library, GeographicLib. Key components of this toolbox are

* Geodesics: direct, inverse, area calculations.
* Projections: transverse Mercator, polar stereographic, etc.
* Grid systems: UTM, UPS, MGRS.
* Geoid lookup: EGM84, EGM96, EGM2008 geoids supported.
* Geometric transformations: geocentric, local cartesian.
* Great ellipses: direct, inverse, area calculations.

There is some overlap between this toolbox and MATLAB's Mapping
Toolbox. However, this toolbox offers:

* better accuracy;
* treatment of oblate and prolate ellipsoids;
* guaranteed convergence for geoddistance;
* calculation of area and differential properties of geodesics;
* ellipsoidal versions of the equidistant azimuthal and gnomonic
projections.

Subsets of this package were previously released as:

Geodesics on an ellipsoid of revolution (deprecated)
Geodesic projections for an ellipsoid (withdrawn)
Great ellipses (withdrawn)

Including all of the functionality in a single toolbox allows easier
sharing of code (via a common private directory).

Extensive documentation on the C++ library is available at

https://geographiclib.sourceforge.io/1.50

Geoid lookup requires the installation of one or more geoid models.
Instructions for this are given in

https://geographiclib.sourceforge.io/1.50/geoid.html#geoidinst

A change log for this package is available at

https://geographiclib.sourceforge.io/1.50/changes.html

Cite As

Charles Karney (2019). geographiclib (https://www.mathworks.com/matlabcentral/fileexchange/50605-geographiclib), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (23)

> Thanks for all your efforts on this library and for supporting the
> MATLAB community. Just getting started so I'm not going to rate this
> yet. One potential issue (not bug) to report is that there are
> differences between Mathworks and geographiclib implementations of the
> ecc2flat function. When geographiclib is on the path earlier than the
> mapping toolbox, then (in MATLAB 2018a at least), an error is thrown
> by distance.m from the mapping toolbox when a 5th input (ellipsoid) is
> specified. This happens regardless of whether a [SemiMajor Ecc] or
> referenceEllipsoid is used. I tested with wgs84Ellipsoid not a
> complex/prolate value.

If there's really an incompatibility between the geographiclib and
mapping toolbox versions of ecc2flat, I'll consider trying to make its
definition (and that of flat2ecc) consistent. However the documentation
for ecc2flat

https://www.mathworks.com/help/map/ref/ecc2flat.html

says that the the calling sequence f = ecc2flat(ellipsoid) "will be
removed in a future release". So perhaps I should just sit tight.

> An unrelated interesting "discovery" is that I get slightly different
> answers from MathWorks distance for the same test points if I convert
> all the points at once vs loop over them individually. Your
> geoddistance is rock steady on that front.

Yes, getting identical results from scalar and vector inputs was a
design goal.

> Note that for my simple test cases, geoddistance ranges from 74%
> (single point conversion) to 244% (10k point conversion) slower than
> distance, but that may be worth getting consistency as well as the
> accuracy improvements that you claim. Performance seems to degrade
> with number of points.

OK, thanks for the information. I don't have access to the mapping
toolbox, so I can't independently verify this.

One case where I think (hope!) that my implementation might be faster is
using geodreckon to compute multiple points on a single geodesic
(calling geodreckon with scalar lat1, lon1, azi1 and vector s12).

Also, I'd be curious how the mapping toolbox does nowadays on accuracy.
Could I ask you to check for me? Download GeodTest.dat from

http://doi.org/10.5281/zenodo.32156

Then

t=load('GeodTest.dat');
s12=t(:,7);
s12a=geoddistance(t(:,1),t(:,2),t(:,4),t(:,5));
s12b=distance(t(:,1),t(:,2),t(:,4),t(:,5),defaultellipsoid)
erra=abs(s12a-s12);
errb=abs(s12b-s12);
max(erra)
sum(~(erra<1e-3))
max(errb)
sum(~(errb<1e-3))
...

(Backwards tests are to catch NaNs.) The fact that distance at one time
might return NaNs on legal inputs disqualified it from serious use.
(For example it can't be used in algorithms which rely on the triangle
inequality.) But I'm not sure whether that's still the case.

Ray

Thanks for all your efforts on this library and for supporting the MATLAB community. Just getting started so I'm not going to rate this yet. One potential issue (not bug) to report is that there are differences between Mathworks and geographiclib implementations of the ecc2flat function. When geographiclib is on the path earlier than the mapping toolbox, then (in MATLAB 2018a at least), an error is thrown by distance.m from the mapping toolbox when a 5th input (ellipsoid) is specified. This happens regardless of whether a [SemiMajor Ecc] or referenceEllipsoid is used. I tested with wgs84Ellipsoid not a complex/prolate value.

An unrelated interesting "discovery" is that I get slightly different answers from MathWorks distance for the same test points if I convert all the points at once vs loop over them individually. Your geoddistance is rock steady on that front.

Note that for my simple test cases, geoddistance ranges from 74% (single point conversion) to 244% (10k point conversion) slower than distance, but that may be worth getting consistency as well as the accuracy improvements that you claim. Performance seems to degrade with number of points.

Thanks again!

Great Package!

help geocent_fwd says that lat and lon should be in degrees!
[X,Y,Z] = geocent_fwd(37.0328,15.065,370)
X = 4922883.48104775
Y = 1325070.28306503
Z = 3820522.46838327

To get the direction relative to one of the stations, you
should use loccart_fwd instead of geocent_fwd.

BlueEyes

Good! Some discrepancies arise between this output coming from your drift
----------------------
>> format long g
>> rad=pi/180;
>> lat=37.0328*rad;
>> lon=15.065*rad;
>> h=370.0;
>> [X,Y,Z]= geocent_fwd(lat,lon,h)
X = 6378036.70547581
Y = 29269.4077049027
Z = 71471.7403795619
and the following coming from https://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/invers3d.prl with results been att'd previously [sendspace works properly, click twice instead] and copied below (1^st station only)
First Station : CB
X = 4922883.4811 m LAT = 37 1 58.08000 North
Y = 1325070.2831 m LON = 15 3 54.00000 East
Z = 3820522.4683 m EHT = 370.0000 Meters
----------
Cheers

The sendspace link seems to require that I create an account and provide an E-mail.
I'm not prepared to do this. So please figure out another way of sending me images.

If you want the Euclidean distance between 2 points you can just convert them to
geocentric coordinates with geocent_fwd and then use Pythagoras.

BlueEyes

geoddistance enhancement proposal: mark-to-mark distance btw. two earth stations, like shown here:
https://www.sendspace.com/file/i5sg1u
Is it possible to do it and link here? Or, alternatively, to show (X,Y,Z) rectangular coords algorithm?
Thx in advance

BlueEyes

Oh yes! It works properly now, the function becomes:
---------
function n = GeodSolve59
% Check for points close with longitudes close to 180 deg apart.
n = 0;
[s12, azi1, azi2] = geoddistance(5, 0.00000000000001, 10, 180);
n = n + assertEquals(azi1, 0.000000000000035, 1.5e-14);
n = n + assertEquals(azi2, 179.99999999999996, 1.5e-14);
n = n + assertEquals(s12, 18345191.174332713, 5e-9);
end
---------

@BlueEyes You are misinterpreting the patch file. The rule is remove the lines beginning
with "-" and replace them with the lines beginning with "+" (but with the "+" removed).
(Basically, the delta in the s12 test is changed from 2.5e-9 to 5e-9.)

BlueEyes

Screenshot of the code amended
---------
https://www.sendspace.com/file/5nr3nm
---------

BlueEyes

My GNU Octave 4.0.0 package still fails, as follows:
-----
>> geographiclib_test
parse error near line 458 of file C:\Training\z_geodesia\karney\geographiclib_test.m

invalid left hand side of assignment

>>> - n = n + assertEquals(s12, 18345191.174332713, 2.5e-9);
-------
Cheers

This patch to geographiclib_test.m fixes the problem reported by BlueEyes (it relaxes
one of the tests slightly):

diff --git a/matlab/geographiclib/geographiclib_test.m b/matlab/geographiclib/geographiclib_test.m
index 2bad2f18..d301465e 100644
--- a/matlab/geographiclib/geographiclib_test.m
+++ b/matlab/geographiclib/geographiclib_test.m
@@ -457,7 +457,7 @@ function n = GeodSolve59
[s12, azi1, azi2] = geoddistance(5, 0.00000000000001, 10, 180);
n = n + assertEquals(azi1, 0.000000000000035, 1.5e-14);
n = n + assertEquals(azi2, 179.99999999999996, 1.5e-14);
- n = n + assertEquals(s12, 18345191.174332713, 2.5e-9);
+ n = n + assertEquals(s12, 18345191.174332713, 5e-9);
end

function n = GeodSolve61

BlueEyes

Just sent the reply via email. Cheers

Reply to @BlueEyes: This may be a problem with the version of Matlab you
are used or the platform on which it is being run. Please could you run
the following in your version of Matlab (or Octave) to help me diagnose
the problem:

ver
[s12, azi1, azi2] = geoddistance(5, 0.00000000000001, 10, 180);
fprintf('%.10f %.15f %.15f\n', s12, azi1, azi2)

and provide the results to me (E-mail is preferred:
charles@karney.com). Thanks.

BlueEyes

Please, be noticed of this drawback:
--------
1..
>> geographiclib_test
GeodSolve59 fail: 1
error: assert (n == 0) failed
error: called from
assert at line 92 column 11
geographiclib_test at line 56 column 3

2.. by eliminating the corresponding row, the test is successful:
% i = GeodSolve59; if i, n=n+1; fprintf('GeodSolve59 fail: %d\n', i); end
>> geographiclib_test
>>
--------
Cheers

John

Fantastic!! I have used geographiclib in Python and am very pleased it is now implemented in MATLAB.

Zhao Dejun

There's an obscure bug in geodreckon in calculating the area when an empty distance argument is given. A patch is provided in

https://sourceforge.net/p/geographiclib/news/2017/03/bug-in-matlab-function-geodreckon/

This will be included in the next release.

Bill Tandy

Exactly what I was looking for! Thank you!

Guy

The update to version 1.44 on 2015-08-14 caused a bug in geoddistance,
when it is used in Octave (not MATLAB) and is invoked with vector
arguments including nearly antipodal pairs of points on the equator.

A patch is given in
https://sourceforge.net/p/geographiclib/news/2015/08/bug--fix-in-octave-version-of-geoddistance/

This will be included in the next release.

Matthew

Chen Qiang

Great!

Updates

1.50

geodarea can now handle arbitrarily complex polygons.
Fix bug in mgrs_inv which resulted in incorrect results for UPS zones with prec = -1.
In geodreckon.m and geoddistance.m, suppress (innocuous) "warning: division by zero" messages from Octave.

1.49

Update to version 1.49

1.48.0.0

Fix BUGS in geodreckon with mixed scalar and array arguments.
Default range for longitude and azimuth is (-180d, 180d].

1.47.0.0

Improve accuracy of area calculation (fixing a flaw introduced in version 1.46).
Fix vectorization of copysignx for MATLAB (Octave was already OK).

1.46.0.0

Improve the accuracy of the solution of the inverse problem when the
longitude difference is close to 180deg.

1.45.0.0

Synchronize with GeographicLib 1.45.
tranmerc_{fwd,inv} works with mixed scalar+array args.
couple of Octave-specific fixes.
array mismatch fix for geoddistance

1.44.1.0

Fix obscure bug in geoddistance (apply to Octave only).

1.44.0.0

Synchronize with GeographicLib 1.44.

1.43.0.0

Synchronize with GeographicLib 1.43.
Fix bug in the long_unroll feature of geodreckon.
mgrs_inv now takes an optional center argument.

1.6.0.0

Reverting to zip packaging. MATLAB Central produces a garbled zip file
when submitting a mltbx file.

1.42.1.0

Version 1.42.1 Repackage as a toolbox to eliminate spurious dependency. No change in .m files.

1.4.0.0

Sync with GeographicLib 1.42. Minor changes to documentation only.

1.3.0.0

Remove bogus dependency on Robust Control Toolbox (yet again!).

1.2.0.0

Remove bogus dependency on Robust Control Toolbox (again!).

1.1.0.0

Remove bogus required product.

MATLAB Release Compatibility
Created with R2014b
Compatible with any release
Platform Compatibility
Windows macOS Linux