geographiclib-users — GeographicLib discussion

[Geographiclib-users] Regarding geographiclib cve scanning

[Pawan B. <Paw...@kp...>]

Hi Team,

We want to scan the CVE for this package.
Could you please provide the CVE_PRODUCT name for geographiclib?

Thanks & Regards,
Pawan
​KPIT Technologies,Pune
This message contains information that may be privileged or confidential and is the property of the KPIT Technologies Ltd. It is intended only for the person to whom it is addressed. If you are not the intended recipient, you are not authorized to read, print, retain copy, disseminate, distribute, or use this message or any part thereof. If you receive this message in error, please notify the sender immediately and delete all copies of this message. KPIT Technologies Ltd. does not accept any liability for virus infected mails.

Re: [Geographiclib-users] Geographiclib and 3D

[Charles K. <cha...@sr...>]

If you start with latitude/longitude/height-above-ellipsoid, then you
can use LocalCartesian::Foward to convert to a normal 3d coordinates.
Use standard 3d geometry to move to the second point according to
heading, pitch, and distance.  Then use LocalCartesian::Reverse to get
back to lat/long/height.


________________________________________
From: Theuns Heydenrych <the...@gm...>
Sent: Tuesday, March 28, 2023 7:39 AM
To: geo...@li...
Subject: [EXTERNAL] [Geographiclib-users] Geographiclib and 3D

Hi
I would like to use Geographiclib to work out a forward location by providing a start location, heading, pitch and a distance.
The start location is in 3D space lat, lon and height, and I would like to work out the next point in 3D space providing the above mentioned parameters.

Is this possible by using Geographiclib?

Regards
Theuns Heydenrych

[Geographiclib-users] Geographiclib and 3D

[Theuns H. <the...@gm...>]

Hi
I would like to use Geographiclib to work out a forward location by
providing a start location, heading, pitch and a distance.
The start location is in 3D space lat, lon and height, and I would like to
work out the next point in 3D space providing the above mentioned
parameters.

Is this possible by using Geographiclib?

Regards
Theuns Heydenrych

[Geographiclib-users] Excessive object allocation in Geodesic

[Andrea A. <and...@gm...>]

Hi,
I am investigating the slowness of a software using Geodesic lib to perform
calculations.
The code is calling Geodesic.inverse, which leads to the following report
by Java Mission Control:

------------------------
The most allocated type is likely 'double[]', most commonly allocated by:
Geodesic.C1f(double, double[]) (34,8 %)
Frequently allocated types are good places to start when trying to reduce
garbage collections. Look at where the most common types are being
allocated to see if many instances are created along the same call path.
Try to reduce the number of instances created by invoking the most commonly
taken paths less.
------------------------

The stack trace looks like the following:

void net.sf.geographiclib.Geodesic.C1f(double, double[]) 16856
Geodesic$LengthsV net.sf.geographiclib.Geodesic.Lengths(double, double,
double, double, double, double, double, double, double, double, int,
double[], double[]) 16852
Geodesic$Lambda12V net.sf.geographiclib.Geodesic.Lambda12(double, double,
double, double, double, double, double, double, double, double, boolean,
double[], double[], double[]) 13132
Geodesic$InverseData net.sf.geographiclib.Geodesic.InverseInt(double,
double, double, double, int) 13132
GeodesicData net.sf.geographiclib.Geodesic.Inverse(double, double, double,
double, int) 13132
GeodesicData net.sf.geographiclib.Geodesic.Inverse(double, double, double,
double) 13132

Let's have a look at the source code of Geodesic.C1f:

// The coefficients C1[l] in the Fourier expansion of B1
  protected static void C1f(double eps, double c[]) {
    final double coeff[] = {
      // C1[1]/eps^1, polynomial in eps2 of order 2
      -1, 6, -16, 32,
      // C1[2]/eps^2, polynomial in eps2 of order 2
      -9, 64, -128, 2048,
      // C1[3]/eps^3, polynomial in eps2 of order 1
      9, -16, 768,
      // C1[4]/eps^4, polynomial in eps2 of order 1
      3, -5, 512,
      // C1[5]/eps^5, polynomial in eps2 of order 0
      -7, 1280,
      // C1[6]/eps^6, polynomial in eps2 of order 0
      -7, 2048,
    };
    double
      eps2 = GeoMath.sq(eps),
      d = eps;
    int o = 0;
    for (int l = 1; l <= nC1_; ++l) { // l is index of C1p[l]
      int m = (nC1_ - l) / 2;         // order of polynomial in eps^2
      c[l] = d * GeoMath.polyval(m, coeff, o, eps2) / coeff[o + m + 1];
      o += m + 2;
      d *= eps;
    }
  }


Now... the only allocation I see is the coeff array... maybe the original
author thought, since it was a final, that it would
be allocated just once? Instead it seems it's allocated for every call (and
it's really mutable, only the object reference itself
is final, arrays in java cannot be made constant, as far as I know).
My guess is that moving these coefficients (and many others used in the
other methods) out in a "static final" field may solve the issue.
Does that make sense?

Btw, I looked for a bug tracker, but could not find one, so I've reported
the issue here.

Cheers
Andrea

Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

On 09/25/2014 04:40 PM, Mike Toews wrote:
> On 26 September 2014 03:56, Charles Karney <cha...@sr...> wrote:
>> You've got the MPFR requirements wrong...
>>
>> I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
>> to MPFR C++ which is a single header file (mpreal.h) you can stick into
>> the include directory.
>>
>
> OK thanks, I see the distinction between MPFR and MPFR C++; I had only
> installed the former. The later is not yet available in Debian stable,
> but is in testing (Debian Jessie). Regardless, it's a simple header
> file that can easily be downloaded and placed in /usr/local/include,
> and the project compiles normally.
>
> -Mike
>

The way to set the MPFR precision in the utility programs is as follows:

   echo '0 0;0 90;90 0' |
   tr ';' '\n' |
   GEOGRAPHICLIB_DIGITS=1000 tools/Planimeter -p 100

which means calculate with 1000 bits of precision and print the results
with ~100 decimal digits of precision.

   --Charles

Re: [Geographiclib-users] Find S12 using planimeter

[Mike T. <mw...@gm...>]

On 26 September 2014 03:56, Charles Karney <cha...@sr...> wrote:
> You've got the MPFR requirements wrong...
>
> I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
> to MPFR C++ which is a single header file (mpreal.h) you can stick into
> the include directory.
>

OK thanks, I see the distinction between MPFR and MPFR C++; I had only
installed the former. The later is not yet available in Debian stable,
but is in testing (Debian Jessie). Regardless, it's a simple header
file that can easily be downloaded and placed in /usr/local/include,
and the project compiles normally.

-Mike

Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

On 2014-09-25 09:01, Mike Toews wrote:
> I'm unable to build a higher precision version. On Debian Wheezy, I
> can't compile with either GEOGRAPHICLIB_PRECISION 4 or 5, since the
> supported releases of Boost 1.49 or MPFR 3.1.0 don't meet the
> requirements for v 1.37 of GeographicLib, which requires either Boost
> 1.54 or MPFR 3.5.9 (in CMakeLists.txt). So I'll test this out some
> other time.
>

You've got the MPFR requirements wrong...

I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
to MPFR C++ which is a single header file (mpreal.h) you can stick into
the include directory.  The MPFR C++ download link is

   http://www.holoborodko.com/pavel/mpfr/#download

I will fix the documentation in CMakeLists.txt and in doxygen to clarify
this.

By the way, at some point as you increase the precision you need to
switch to using Planimeter with the -E option.  The geodesics are then
computed with elliptic integrals.  The area calculation uses a 30th
order series (I never figured out how to express the are integral in
terms of elliptic integrals); so eventually you will run out of accuracy
on areas (but it's probably good for at least 1000 digits for WGS84).

   --Charles

Re: [Geographiclib-users] Find S12 using planimeter

[Mike T. <mw...@gm...>]

On 19 September 2014 01:54, Charles Karney <cha...@sr...> wrote:
> If you want different rules applied for non-simple polygons, then you
> must decompose such polygons into the union of simple polygons.  E.g.,
> change the bow tie quadrilateral
>
>     1 -1
>     -1 -1
>     1 1
>     -1 1
>
> to the hexagon
>
>     1 -1
>     -1 -1
>     0 0
>     -1 1
>     1 1
>     0 0

But of course finding the intersection point (i.e. 0 0) to make the
valid geometry is something altogether different! And yes, I see the
basic approach using a Gnomic projection. What I was looking to see if
computing an area on a non-simple polygon provides a similar answer as
other algorithms, and that it doesn't catch an unexpected exception
(e.g. divide by zero).

On 19 September 2014 10:27, Charles Karney <cha...@sr...> wrote:
> The errors you get will very much depend on the type of polygon.  0.1m^2
> per edge is a worst case; it appears that the errors cancel to some
> extent.
>
> Version 1.37 allows you to compile GeographicLib with Boost or MPFR
> allowing you to use either quad or arbitrary precision.  So perhaps you
> can gauge the errors in your case yourself.  (I can possibly get you a
> binary version of GeographicLib + MPFR for Windows.)  If you are on a
> Linux machine you can also compile using long doubles and get an extra
> 11 bits of precision.
>
> As a practical example...  I have a dataset of the boundaries of Poland
> which has
>
>     67801 edges
>     mean edge 53m
>     min edge 0.06m
>     max edge 10km
>
> Planimeter gives the area as 312679715911.98608 m^2.  The true area is
> 312679715911.9856346137 m^2.  So the error in this case is less than
> 0.001 m^2.

These errors are small enough to me. I might expect the largest errors
to be for simplified polygons with few vertices and many sharp angles.

I'm unable to build a higher precision version. On Debian Wheezy, I
can't compile with either GEOGRAPHICLIB_PRECISION 4 or 5, since the
supported releases of Boost 1.49 or MPFR 3.1.0 don't meet the
requirements for v 1.37 of GeographicLib, which requires either Boost
1.54 or MPFR 3.5.9 (in CMakeLists.txt). So I'll test this out some
other time.

-Mike

Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

The errors you get will very much depend on the type of polygon.  0.1m^2
per edge is a worst case; it appears that the errors cancel to some
extent.

Version 1.37 allows you to compile GeographicLib with Boost or MPFR
allowing you to use either quad or arbitrary precision.  So perhaps you
can gauge the errors in your case yourself.  (I can possibly get you a
binary version of GeographicLib + MPFR for Windows.)  If you are on a
Linux machine you can also compile using long doubles and get an extra
11 bits of precision.

As a practical example...  I have a dataset of the boundaries of Poland
which has

     67801 edges
     mean edge 53m
     min edge 0.06m
     max edge 10km

Planimeter gives the area as 312679715911.98608 m^2.  The true area is
312679715911.9856346137 m^2.  So the error in this case is less than
0.001 m^2.

   --Charles

On 2014-09-18 15:39, Mike Toews wrote:
> Great, thanks, I think that clears up my questions with my assumptions
> on the quadrilateral for S12. And yes, the absolute error of residual
> with the technique described previously increases approaching the
> (1-f)*180 deg edges.
>
> The planimeter approach is what I'm interested in. Out of curiosity,
> are there any test datasets that can be used to check the accuracy of
> polygon area calculations? Also, how does the error of polygon area
> calculation accumulate? For instance, I see the expected error of S12
> is 0.1 m^2. Does the error relate to the number of vertices in the
> polygon in combination to the length of each edge? I just looking for
> a general idea, nothing exhaustive.
>
> Thanks in advance,
> -Mike
>

-- 
Charles Karney <cha...@sr...>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

Re: [Geographiclib-users] Find S12 using planimeter

[Mike T. <mw...@gm...>]

Great, thanks, I think that clears up my questions with my assumptions
on the quadrilateral for S12. And yes, the absolute error of residual
with the technique described previously increases approaching the
(1-f)*180 deg edges.

The planimeter approach is what I'm interested in. Out of curiosity,
are there any test datasets that can be used to check the accuracy of
polygon area calculations? Also, how does the error of polygon area
calculation accumulate? For instance, I see the expected error of S12
is 0.1 m^2. Does the error relate to the number of vertices in the
polygon in combination to the length of each edge? I just looking for
a general idea, nothing exhaustive.

Thanks in advance,
-Mike

Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

On 2014-09-18 09:54, Charles Karney wrote:
> However the geodesic HB only follows the equator if the longitude
> difference is less than (1-f)*180 deg = 179.396494 deg.

Sorry, this should be:

However the geodesic FH only follows the equator if the longitude
difference is less than (1-f)*180 deg = 179.396494 deg.

-- 
Charles Karney <cha...@sr...>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

You raise four issues:

(A) The cause of the discrepancy between GeodSolve's area and that
computed by Planimeter is as follows:

GeodSolve -f reports the area of the quadrilateral

     AF along a meridian
     FH along the equator
     HB along a meridian
     BA along a geodesic

In contrast, Planimeter reports the area of the quadrilateral

     AF along the shortest geodesic
     FH along the shortest geodesic
     HB along the shortest geodesic
     BA along the shortest geodesic

The geodesics AF and HB will be along meridians (for oblate ellipsoids).
However the geodesic HB only follows the equator if the longitude
difference is less than (1-f)*180 deg = 179.396494 deg.  In your
example, it picks a shorter geodesic through the northern hemisphere.
In fact, there are two equally short geodesics and it might be picked
the other one.

So, here's what the man page for Planimeter has to say about this:

        The edges of the polygon are given by the shortest geodesic
        between consecutive vertices.  In certain cases, there may be two
        or many such shortest geodesics, and in that case, the polygon is
        not uniquely specified by its vertices.  This only happens with
        very long edges (for the WGS84 ellipsoid, any edge shorter than
        19970 km is uniquely specified by its end points).  In such
        cases, insert an additional vertex near the middle of the long
        edge to define the boundary of the polygon.

(B) As far as getting the same results with the C library as with
GeodSolve,

Call either geod_geninverse or geod_gendirect.  As you can see from

   http://geographiclib.sourceforge.net/html/C/geodesic_8h.html

These take an optional pointer argument pS12 which returns the same area
as reported by GeodSolve.

(C) Do you need to worry about simple/non-simple polygons?  Definitely.
Areas for non-simple polygons can be defined in multiple ways and none
of them is the "right" way.

For simple polygons (not self intersecting), the area reported is
positive for counter-clockwise traversal.  If this area exceeds half the
area of the ellipsoid, the negative of the "rest" of the area of the
ellipsoid is returned.  (The -s and -r flags to Planimeter govern this
behavior.)  The poles of the ellipsoid are taken care of for you.

For non-simple polygons, the simple calculus definition is applied.  For
for the bow tie example, the triangle which is included
counter-clockwise counts as a positive area and the other one as a
negative area and the sum will be close to zero for a symmetric bow tie.
Similarly the area a polygon which encircles a point twice will be two
times the area of the same sized polygon which encircle it one.

If you want different rules applied for non-simple polygons, then you
must decompose such polygons into the union of simple polygons.  E.g.,
change the bow tie quadrilateral

     1 -1
     -1 -1
     1 1
     -1 1

to the hexagon

     1 -1
     -1 -1
     0 0
     -1 1
     1 1
     0 0

(D) Do you needs to worry about the poles?  No, not if you're using
Planimeter.  You can supply the boundary of Antarctica to Planimeter and
get the correct area.

I hope this answers your questions.

   --Charles

On 2014-09-18 07:03, Mike Toews wrote:
> Hi,
>
> I'm trying to use the C library (v1.32) to find polygon areas, but I'm
> having difficulty getting consistent results with the GeodTest.dat
> data set. From what I understand, to find the area under the geodesic,
> S12, the AFHB quadrilateral can be built with corners (lat1,lon1),
> (0,lon1), (0,lon2), and (lat2,lon2).
>
> Example 1 is line 1 of GeodTest.dat
>
> $ sed -n 1p GeodTest.dat
> 36.530042355041 0 176.125875162171 -48.164270779097768864
> 5.762344694676510456 175.334308316285410561 9398502.0434687
> 84.663858149358862201 6333544.7732452481809 -559418252332.321555
>
> or input for planimeter:
>
> 36.530042355041 0
> 0 0
> 0 5.762344694676510456
> -48.164270779097768864 5.762344694676510456
>
> With output:
> 4 19421014.749204 -559418252332.3
>
> This is great, since the areas are within the expected precision.
> However, do note that this quadrilateral is not simple, since it has a
> bow-tie, which is understood not be used by this function. Similar
> success with absolute residuals <= 0.1 m^2 is found from 55.3498% of
> the 500k lines in the dataset, even with simple quadrilaterals like
> line 1270.
>
> Example 2 is line 381 of GeodTest.dat, which forms a simple quadrilateral
>
> $ sed -n 381p GeodTest.dat
> 42.496576764546 0 .165548738367 61.200095175201199603
> 179.666196014678657282 179.746893753529690115 8509305.8221784
> 76.480111135775561028 6211557.6845057436865 127219347433402.487683
>
> or input for planimeter:
>
> 42.496576764546 0
> 0 0
> 0 179.666196014678657282
> 61.200095175201199603 179.666196014678657282
>
> With output:
> 4 39997539.52627776 47480061825561.484
>
> This isn't too great, since the absolute residual of area calculations
> is 79739285607841. or a ratio of about 0.373
>
> Now I know that GeographicLib (v1.37) can get the correct S12 for this
> example to GeodTest.dat, e.g.:
>
> $ ./GeodSolve -i -f
> 42.496576764546 0 61.200095175201199603 179.666196014678657282
> 42.49657676 0.00000000 0.16554874 61.20009518 179.66619601
> 179.74689375 8509305.822 76.48011114 6211557.685 0.2368660383
> 0.2380528208 127219347433402
>
> However, with the C++ tool Planimeter from this same version (v1.37),
> I get similar results to the C library. The only difference between
> the two is the precision of output, which is really just a formatting
> difference. So I suspect I'm doing my area under the geodesic
> calculation with the planimeter utilities incorrectly.
>
> My question is how can I achieve the same answer to GeodSolve with the
> C library? Also, do I need to worry about simple/non-simple polygons?
> Is there any further post-analysis if the polygon encircles a pole (in
> regards to last paragraph of section 6 of "Algorithms for geodesics"
> paper), or does the library handle these situations seamlessly?
>
> Many thanks,
>
> -Mike


-- 
Charles Karney <cha...@sr...>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

[Geographiclib-users] Find S12 using planimeter

[Mike T. <mw...@gm...>]

Hi,

I'm trying to use the C library (v1.32) to find polygon areas, but I'm
having difficulty getting consistent results with the GeodTest.dat
data set. From what I understand, to find the area under the geodesic,
S12, the AFHB quadrilateral can be built with corners (lat1,lon1),
(0,lon1), (0,lon2), and (lat2,lon2).

Example 1 is line 1 of GeodTest.dat

$ sed -n 1p GeodTest.dat
36.530042355041 0 176.125875162171 -48.164270779097768864
5.762344694676510456 175.334308316285410561 9398502.0434687
84.663858149358862201 6333544.7732452481809 -559418252332.321555

or input for planimeter:

36.530042355041 0
0 0
0 5.762344694676510456
-48.164270779097768864 5.762344694676510456

With output:
4 19421014.749204 -559418252332.3

This is great, since the areas are within the expected precision.
However, do note that this quadrilateral is not simple, since it has a
bow-tie, which is understood not be used by this function. Similar
success with absolute residuals <= 0.1 m^2 is found from 55.3498% of
the 500k lines in the dataset, even with simple quadrilaterals like
line 1270.

Example 2 is line 381 of GeodTest.dat, which forms a simple quadrilateral

$ sed -n 381p GeodTest.dat
42.496576764546 0 .165548738367 61.200095175201199603
179.666196014678657282 179.746893753529690115 8509305.8221784
76.480111135775561028 6211557.6845057436865 127219347433402.487683

or input for planimeter:

42.496576764546 0
0 0
0 179.666196014678657282
61.200095175201199603 179.666196014678657282

With output:
4 39997539.52627776 47480061825561.484

This isn't too great, since the absolute residual of area calculations
is 79739285607841. or a ratio of about 0.373

Now I know that GeographicLib (v1.37) can get the correct S12 for this
example to GeodTest.dat, e.g.:

$ ./GeodSolve -i -f
42.496576764546 0 61.200095175201199603 179.666196014678657282
42.49657676 0.00000000 0.16554874 61.20009518 179.66619601
179.74689375 8509305.822 76.48011114 6211557.685 0.2368660383
0.2380528208 127219347433402

However, with the C++ tool Planimeter from this same version (v1.37),
I get similar results to the C library. The only difference between
the two is the precision of output, which is really just a formatting
difference. So I suspect I'm doing my area under the geodesic
calculation with the planimeter utilities incorrectly.

My question is how can I achieve the same answer to GeodSolve with the
C library? Also, do I need to worry about simple/non-simple polygons?
Is there any further post-analysis if the polygon encircles a pole (in
regards to last paragraph of section 6 of "Algorithms for geodesics"
paper), or does the library handle these situations seamlessly?

Many thanks,

-Mike

Re: [Geographiclib-users] Comparing distances between geographic coordinates

[Charles K. <cha...@sr...>]

Most people are poor at knowing where the bottlenecks in their code will
be.  So my recommendation would be to use the geodesic distance and
only worry about finding a substitute if that proves too slow.

On 05/18/2014 12:33 PM, Edward Lam wrote:
> Hi Charles,
>
> Yes, I want the metric to correlate to real distances. Some of the hardware I'm looking at can be very underpowered (think 400 MHz clock speeds) which may (or may not) matter. I figured I'd ask anyway if there some proxy distance measure I could use that was faster to compute.
>
> Thanks for your suggestion regardless.
>
> -Edward
>
>> On May 17, 2014, at 5:04 PM, Charles Karney <cha...@sr...> wrote:
>>
>> Lots of distance metric will satisfy the triangle inequality
>>
>> * geodesic distance on ellipsoid
>> * Euclidean distance on ellipsoid
>> * surface distance on sphere
>> * Euclidean distance on sphere
>> * a Manhattan distance |dlat| + |dlon|
>>
>> etc.  However, I suspect you want a stronger condition, namely
>>
>>   d(A,B) < d(A,C)
>>
>> implies that B really is closer to A than C.  For this I recommend just
>> using the true geodesic distance (e.g., from GeographicLib).  You can do
>> a million such distance calculations in a couple of seconds.
>>
>> On 05/17/2014 03:45 PM, Edward Lam wrote:
>>> Hi,
>>>
>>> I trying to write an app that performs route planning through a graph of
>>> points given by their GPS latitude/longitude coordinates. So to do this,
>>> I need to compare the relative distances between these points such that
>>> the triangle inequality holds. What is the best/fastest way to do so?
>>>
>>> There's a great deal of description on the web on how to compute
>>> distances between two GPS coordinates ranging from approximate ones
>>> based on the haversine formula via an idealized sphere to more accurate
>>> ellipsoidal ones like the one in GeographicLib.
>>>
>>> For this application, I don't need real distances between the points,
>>> just some metric so that the triangle inequality holds. I'm tempted to
>>> use the haversine formula since it is cheaper than the alternatives but
>>> it has distortions depending on the chosen radius that may affect the
>>> triangle inequality? Or am I over thinking this? Is there some
>>> projection I can use that is both cheap and accurate?
>>>
>>> Thanks,
>>> -Edward
>>

Re: [Geographiclib-users] Comparing distances between geographic coordinates

[Edward L. <e4...@ya...>]

Hi Charles,

Yes, I want the metric to correlate to real distances. Some of the hardware I'm looking at can be very underpowered (think 400 MHz clock speeds) which may (or may not) matter. I figured I'd ask anyway if there some proxy distance measure I could use that was faster to compute.

Thanks for your suggestion regardless.

-Edward

> On May 17, 2014, at 5:04 PM, Charles Karney <cha...@sr...> wrote:
> 
> Lots of distance metric will satisfy the triangle inequality
> 
> * geodesic distance on ellipsoid
> * Euclidean distance on ellipsoid
> * surface distance on sphere
> * Euclidean distance on sphere
> * a Manhattan distance |dlat| + |dlon|
> 
> etc.  However, I suspect you want a stronger condition, namely
> 
>  d(A,B) < d(A,C)
> 
> implies that B really is closer to A than C.  For this I recommend just
> using the true geodesic distance (e.g., from GeographicLib).  You can do
> a million such distance calculations in a couple of seconds.
> 
> On 05/17/2014 03:45 PM, Edward Lam wrote:
>> Hi,
>> 
>> I trying to write an app that performs route planning through a graph of
>> points given by their GPS latitude/longitude coordinates. So to do this,
>> I need to compare the relative distances between these points such that
>> the triangle inequality holds. What is the best/fastest way to do so?
>> 
>> There's a great deal of description on the web on how to compute
>> distances between two GPS coordinates ranging from approximate ones
>> based on the haversine formula via an idealized sphere to more accurate
>> ellipsoidal ones like the one in GeographicLib.
>> 
>> For this application, I don't need real distances between the points,
>> just some metric so that the triangle inequality holds. I'm tempted to
>> use the haversine formula since it is cheaper than the alternatives but
>> it has distortions depending on the chosen radius that may affect the
>> triangle inequality? Or am I over thinking this? Is there some
>> projection I can use that is both cheap and accurate?
>> 
>> Thanks,
>> -Edward
> 

Re: [Geographiclib-users] Comparing distances between geographic coordinates

[Charles K. <cha...@sr...>]

Lots of distance metric will satisfy the triangle inequality

* geodesic distance on ellipsoid
* Euclidean distance on ellipsoid
* surface distance on sphere
* Euclidean distance on sphere
* a Manhattan distance |dlat| + |dlon|

etc.  However, I suspect you want a stronger condition, namely

   d(A,B) < d(A,C)

implies that B really is closer to A than C.  For this I recommend just
using the true geodesic distance (e.g., from GeographicLib).  You can do
a million such distance calculations in a couple of seconds.

On 05/17/2014 03:45 PM, Edward Lam wrote:
> Hi,
>
> I trying to write an app that performs route planning through a graph of
> points given by their GPS latitude/longitude coordinates. So to do this,
> I need to compare the relative distances between these points such that
> the triangle inequality holds. What is the best/fastest way to do so?
>
> There's a great deal of description on the web on how to compute
> distances between two GPS coordinates ranging from approximate ones
> based on the haversine formula via an idealized sphere to more accurate
> ellipsoidal ones like the one in GeographicLib.
>
> For this application, I don't need real distances between the points,
> just some metric so that the triangle inequality holds. I'm tempted to
> use the haversine formula since it is cheaper than the alternatives but
> it has distortions depending on the chosen radius that may affect the
> triangle inequality? Or am I over thinking this? Is there some
> projection I can use that is both cheap and accurate?
>
> Thanks,
> -Edward
>

[Geographiclib-users] Comparing distances between geographic coordinates

[Edward L. <e4...@ya...>]

Hi,

I trying to write an app that performs route planning through a graph of points given by their GPS latitude/longitude coordinates. So to do this, I need to compare the relative distances between these points such that the triangle inequality holds. What is the best/fastest way to do so?

There's a great deal of description on the web on how to compute distances between two GPS coordinates ranging from approximate ones based on the haversine formula via an idealized sphere to more accurate ellipsoidal ones like the one in GeographicLib.

For this application, I don't need real distances between the points, just some metric so that the triangle inequality holds. I'm tempted to use the haversine formula since it is cheaper than the alternatives but it has distortions depending on the chosen radius that may affect the triangle inequality? Or am I over thinking this? Is there some projection I can use that is both cheap and accurate?

Thanks,
-Edward 

Re: [Geographiclib-users] Projecting ECEF into UTM

[Charles K. <cha...@sr...>]

I have this change checked in.  It will be part of 1.35.

On 2014-01-10 14:47, Ben Adler wrote:
> Great. Please tell me in case this fix will *NOT* make it into 1.35.

-- 
Charles Karney <cha...@sr...>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

Re: [Geographiclib-users] Projecting ECEF into UTM

[Ben A. <ad...@in...>]

Hello,

 > Please give an example for code that your compiler rejects.  I believe
> [...]

I tried and cannot, I must have confused that with some other 
parameters. Sorry.

>> Why does GeographicLib::UTMUPS::EncodeEPSG(11, true) return 32610?
>> Shouldn't it be 32611 (http://spatialreference.org/ref/epsg/32611/)?
>
> That's a careless bug.  Please apply the attached patch.  (Only the part
> affecting EncodeEPSG is essential.) Thanks for catching this.

Great. Please tell me in case this fix will *NOT* make it into 1.35.

Charles, that was great support! Thank you.

Ben

Re: [Geographiclib-users] Projecting ECEF into UTM

[Charles K. <cha...@sr...>]

On 2014-01-10 11:58, Ben Adler wrote:
> Charles,
>
> brilliant! Thank you for these quick and helpful answers! I tried
> implementing them, here's my problems:
>
> Why do the GeographicLib::UTMUPS::Forward/Reverse method not use const
> for input parameters? That breaks a lot of my code. Of course I could
> make relevant members mutable or use const_cast, but that seems weird.

Please give an example for code that your compiler rejects.  I believe
that

     static void Forward(real lat, real lon, ...) { ... }

has the same signature as

     static void Forward(const real lat, const real lon, ...) { ... }

(I.e., the const says how lat/lon can be used inside the function but
nothing about what happens outside the function.)

> Now that I have the EPSG code, I need to get the WKT. Grepping for wkt
> in the GeographicLib sources didn't yield any results, so I guess I need
> to search and extract the relevant line from EPSG2WKT.TXT (google-hit at
> geospatialpython) myself?
>
> Comparing that file with the WKT from http://spatialreference.org, the
> text file seems to be missing the AXIS-entries. Also, of course, the
> height information is missing. Maybe it's out of GeographicLib's scope,
> but I wonder how exactly the WKT for e.g. UTM 11N with height above the
> WGS84 ellipsoid in meters should look like.

Probably the AXIS entries are optional?  These EPSG codes deal only with
the horizontal axes.  There are separate (less well established)
conventions for the vertical datum.

> Why does GeographicLib::UTMUPS::EncodeEPSG(11, true) return 32610?
> Shouldn't it be 32611 (http://spatialreference.org/ref/epsg/32611/)?

That's a careless bug.  Please apply the attached patch.  (Only the part
affecting EncodeEPSG is essential.) Thanks for catching this.

   --Charles

Re: [Geographiclib-users] Projecting ECEF into UTM

[Ben A. <ad...@in...>]

Charles,

brilliant! Thank you for these quick and helpful answers! I tried 
implementing them, here's my problems:

Why do the GeographicLib::UTMUPS::Forward/Reverse method not use const 
for input parameters? That breaks a lot of my code. Of course I could 
make relevant members mutable or use const_cast, but that seems weird.

Now that I have the EPSG code, I need to get the WKT. Grepping for wkt 
in the GeographicLib sources didn't yield any results, so I guess I need 
to search and extract the relevant line from EPSG2WKT.TXT (google-hit at 
geospatialpython) myself?

Comparing that file with the WKT from http://spatialreference.org, the 
text file seems to be missing the AXIS-entries. Also, of course, the 
height information is missing. Maybe it's out of GeographicLib's scope, 
but I wonder how exactly the WKT for e.g. UTM 11N with height above the 
WGS84 ellipsoid in meters should look like.

Why does GeographicLib::UTMUPS::EncodeEPSG(11, true) return 32610? 
Shouldn't it be 32611 (http://spatialreference.org/ref/epsg/32611/)?

Thanks again!
Ben

Re: [Geographiclib-users] Projecting ECEF into UTM

[Charles K. <cha...@sr...>]

Ben,

* You need to convert geocentric coordinates to geographic coordinates
   using the Geocentric class (Geocentric::WGS84 is a static instance
   for the WGS84 ellipsoid).

* UTMUPS::Forward converts from geographic coordinates to UTM (or UPS
   near the poles).  This knows the rules for zone assignments (including
   the Norway and Svalbard exceptions and the UTM/UPS boundary).  The
   "setzone" argument to this function lets you override the standard
   assignment.

* UTMUPS::EncodeEPSG converts a zone+hemisphere to an EPSG code.

* The TransverseMercator class just does the actual projection.  UTMUPS
   calls this class and additionally takes care of the zone and the false
   origin.  Typically you won't need to invoke the TransverseMercator
   class directly

* The TransverseMercatorExact class is a variant of TransverseMercator
   which allows the projection to cover the entire globe.  You don't need
   to use this if you're just converting to UTM.

* The geodetic height returned by Geocentric::WGS84.Reverse is the
   height above the ellipsoid.  This might be what you want.  (There's no
   such thing as "UTM height".)  If you need the height about the geoid,
   use the Geoid class to determine the height of the geoid above the
   ellipsoid.  You will need to install the data for the geoid, see

     http://geographiclib.sourceforge.net/html/geoid.html

* "Latitude", without further qualification, always mean geographic
   (geodetic) latitude in GeographicLib.

* My guess is the GeographicLib is a better fit to your needs than
   proj.4.  The strength of proj.4 is in dealing with a wide range of
   projections.

Here's the outline of what you want to do.  See also the example code
for the Geocentric and UTMUPS classes.

    #include <GeographicLib/Geocentric.hpp>
    #include <GeographicLib/UTMUPS.hpp>
    #include <GeographicLib/Geoid.hpp>

    using namespace GeographicLib;

    double X = ..., Y = ..., Z = ...;
    double lat, lon, hae;
    Geocentric::WGS84.Reverse(X, Y, Z, lat, lon, h);
    int zone;
    bool northp;
    double x, y;
    UTMUPS::Forward(lat, lon, zone, northp, x, y);
    int epsg = UTMUPS::EncodeEPSG(zone, northp);
    Geoid egm96("egm96-5");
    hmsl = hae + Geoid::ELLIPSOIDTOGEOID * egm96(lat, lon);

On 01/10/2014 06:38 AM, Ben Adler wrote:
> Hello GeographicLib,
>
> I'm rather new to geodesy and trying to project ECEF/geocentric
> coordinates into UTMUPS/WGS84 with the height given relative to the
> WGS84 ellipsoid. And I have a few questions:
>
> 1) For a given set of ECEF coordinates, I first need to get the
> corresponding UTM zone's central meridian and description. This might
> sound trivial at first, but then there's Norway and all the other
> exceptions close to the north pole. And there's UPS for the poles.
>
>    a) How can I do this with GeographicLib? Which classes/methods should
> I use?
>
>    b) How can I get the WKT or EPSG code for the chosen UTM/UPS zone?
>
> 2) Looking at GeographicLib::TransverseMercatorExact, it seems I must
> input lat/lon in order to retrieve the desired UTM coordinates. This
> means I must first convert ECEF into geodetic, which is also handled by
> GeographicLib. But,
>
>    a) do I need to input geodetic or geocentric latitude?
>
>    b) The docs say "The latitude of origin is taken to be the equator".
> While true for northern zones, Does that mean I need to introduce a
> false northing (of -10M meters) in southern zones when using GeographicLib?
>
>    c) can I just re-use the geodetic WGS84 height as UTM height without
> conversion? I think so, but I have erred often in the past :)
>
> I'm not sure whether GeographicLib is suited to my needs, or whether I
> should use another lib for this, e.g. PROJ.4. Any hints on this?
>
> thanks!
> ben
> _______________________________________________
> Geographiclib-users mailing list
> Geo...@li...
> https://lists.sourceforge.net/lists/listinfo/geographiclib-users
>

[Geographiclib-users] Projecting ECEF into UTM

[Ben A. <ad...@in...>]

Hello GeographicLib,

I'm rather new to geodesy and trying to project ECEF/geocentric 
coordinates into UTMUPS/WGS84 with the height given relative to the 
WGS84 ellipsoid. And I have a few questions:

1) For a given set of ECEF coordinates, I first need to get the 
corresponding UTM zone's central meridian and description. This might 
sound trivial at first, but then there's Norway and all the other 
exceptions close to the north pole. And there's UPS for the poles.

  a) How can I do this with GeographicLib? Which classes/methods should 
I use?

  b) How can I get the WKT or EPSG code for the chosen UTM/UPS zone?

2) Looking at GeographicLib::TransverseMercatorExact, it seems I must 
input lat/lon in order to retrieve the desired UTM coordinates. This 
means I must first convert ECEF into geodetic, which is also handled by 
GeographicLib. But,

  a) do I need to input geodetic or geocentric latitude?

  b) The docs say "The latitude of origin is taken to be the equator". 
While true for northern zones, Does that mean I need to introduce a 
false northing (of -10M meters) in southern zones when using GeographicLib?

  c) can I just re-use the geodetic WGS84 height as UTM height without 
conversion? I think so, but I have erred often in the past :)

I'm not sure whether GeographicLib is suited to my needs, or whether I 
should use another lib for this, e.g. PROJ.4. Any hints on this?

thanks!
ben

Re: [Geographiclib-users] defining a projection (modified)

[Charles K. <cha...@sr...>]

Hugh,

You can perform this projection (but without the false easting and false
northing offsets) with

   TransverseMercatorProj -l 145 -e 6378160 1/298.25 -k 1

E.g.,

echo -37 145.5 | TransverseMercatorProj -l 145 -e 6378160 1/298.25 -k 1
->
44506.1549112743 -4096642.0780917532 -0.3009124466146239 1.0000243910016182

echo 44506.1549112743 -4096642.0780917532 |
TransverseMercatorProj -l 145 -e 6378160 1/298.25 -k 1 -r
->
-37.000000000000007 145.500000000000000 -0.3009124466146240 
1.0000243910016182

The 3rd and 4th numbers in the output are the meridian convergence and
scale.  Applying the false easting and false northing offsets is
straightforward.  Use "TransverseMercatorProj --help" for full
documentation on this utility.

You can also carry out the projection in C++ code using the
TransverseMercator class.

   --Charles

On Mon, Dec 10, 2012 at 4:19 AM, hugh <hdixon@...> wrote:
 > Hi,
 > I am wanting to convert eastings and northings into lats and longs.  I
 > am a civil engineer, not a surveyor, but have dabbled with projections
 > in the past so know a little.
 >
 > I have been told my E / N are in a pseudo AMG (UTM) projection.  I
 > have been given the following parameters:
 >
 > PSEUDO AMG 1966
 >
 > Projection: Transverse_Mercator
 > False_Easting: 500000.000000
 > False_Northing: 10000000.000000
 > Central_Meridian: 145.000000
 > Scale_Factor: 1.000000
 > Latitude_Of_Origin: 0.000000
 > Linear Unit: Meter (1.000000)
 >
 > Geographic Coordinate System: GCS_Australian_1966
 > Angular Unit: Degree (0.017453292519943299)
 > Prime Meridian: Greenwich (0.000000000000000000)
 > Datum: D_Australian_1966
 >   Spheroid: Australian
 >     Semimajor Axis: 6378160.000000000000000000
 >     Semiminor Axis: 6356774.719195305400000000
 >     Inverse Flattening: 298.250000000000000000
 >
 >
 > To me this looks like the standard UTM, except for the scale factor.
 > Is this correct?
 >
 > Any someone suggest an open source library I can use to convert these
 > co ordinates into lats and longs.
 >
 > Thanks
 > Hugh

Re: [Geographiclib-users] defining a projection (modified)

[Shaun W. <sha...@gm...>]

Proj.4 is widely used for this, and might be a good starting point.
Transverse mercator is '+proj=tmerc' in Proj.4 parlance:
  http://trac.osgeo.org/proj/

On Mon, Dec 10, 2012 at 4:19 AM, hugh <hd...@bi...> wrote:

> Hi,
> I am wanting to convert eastings and northings into lats and longs.
> I am a civil engineer, not a surveyor, but have dabbled with
> projections in the past so know a little.
>
> I have been told my E / N are in a pseudo AMG (UTM) projection.  I have
> been given the following parameters:
>
>
> PSEUDO AMG 1966
>
> Projection: Transverse_Mercator
> False_Easting: 500000.000000
> False_Northing: 10000000.000000
> Central_Meridian: 145.000000
> Scale_Factor: 1.000000
> Latitude_Of_Origin: 0.000000
> Linear Unit: Meter (1.000000)
>
> Geographic Coordinate System: GCS_Australian_1966
> Angular Unit: Degree (0.017453292519943299)
> Prime Meridian: Greenwich (0.000000000000000000)
> Datum: D_Australian_1966
>   Spheroid: Australian
>     Semimajor Axis: 6378160.000000000000000000
>     Semiminor Axis: 6356774.719195305400000000
>     Inverse Flattening: 298.250000000000000000
>
>
> To me this looks like the standard UTM, except for the scale factor.
> Is this correct?
>
> Any someone suggest an open source library I can use to convert these
> co ordinates into lats and longs.
>
> Thanks
> Hugh
>
>
> ------------------------------------------------------------------------------
> LogMeIn Rescue: Anywhere, Anytime Remote support for IT. Free Trial
> Remotely access PCs and mobile devices and provide instant support
> Improve your efficiency, and focus on delivering more value-add services
> Discover what IT Professionals Know. Rescue delivers
> http://p.sf.net/sfu/logmein_12329d2d
> _______________________________________________
> Geographiclib-users mailing list
> Geo...@li...
> https://lists.sourceforge.net/lists/listinfo/geographiclib-users
>