You can subscribe to this list here.
2012 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}
(1) 
_{Dec}
(4) 

2013 
_{Jan}

_{Feb}

_{Mar}

_{Apr}
(1) 
_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

2014 
_{Jan}
(6) 
_{Feb}

_{Mar}

_{Apr}

_{May}
(4) 
_{Jun}

_{Jul}

_{Aug}

_{Sep}
(9) 
_{Oct}

_{Nov}

_{Dec}

From: Charles Karney <charles.karney@sr...>  20140926 00:42:42

On 09/25/2014 04:40 PM, Mike Toews wrote: > On 26 September 2014 03:56, Charles Karney <charles.karney@...> 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 
From: Mike Toews <mwtoews@gm...>  20140925 20:40:57

On 26 September 2014 03:56, Charles Karney <charles.karney@...> 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 
From: Charles Karney <charles.karney@sr...>  20140925 15:56:51

On 20140925 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 
From: Mike Toews <mwtoews@gm...>  20140925 13:02:47

On 19 September 2014 01:54, Charles Karney <charles.karney@...> wrote: > If you want different rules applied for nonsimple 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 nonsimple 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 <charles.karney@...> 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 
From: Charles Karney <charles.karney@sr...>  20140918 22:27:26

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 20140918 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 > (1f)*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 <charles.karney@...> SRI International, Princeton, NJ 085435300 Tel: +1 609 734 2312 Fax: +1 609 734 2662 
From: Mike Toews <mwtoews@gm...>  20140918 19:40:12

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 (1f)*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 
From: Charles Karney <charles.karney@sr...>  20140918 13:57:21

On 20140918 09:54, Charles Karney wrote: > However the geodesic HB only follows the equator if the longitude > difference is less than (1f)*180 deg = 179.396494 deg. Sorry, this should be: However the geodesic FH only follows the equator if the longitude difference is less than (1f)*180 deg = 179.396494 deg.  Charles Karney <charles.karney@...> SRI International, Princeton, NJ 085435300 Tel: +1 609 734 2312 Fax: +1 609 734 2662 
From: Charles Karney <charles.karney@sr...>  20140918 13:54:31

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 (1f)*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/nonsimple polygons? Definitely. Areas for nonsimple 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 counterclockwise 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 nonsimple polygons, the simple calculus definition is applied. For for the bow tie example, the triangle which is included counterclockwise 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 nonsimple 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 20140918 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 > bowtie, 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/nonsimple polygons? > Is there any further postanalysis 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 <charles.karney@...> SRI International, Princeton, NJ 085435300 Tel: +1 609 734 2312 Fax: +1 609 734 2662 
From: Mike Toews <mwtoews@gm...>  20140918 11:03:50

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 bowtie, 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/nonsimple polygons? Is there any further postanalysis 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 
From: Charles Karney <charles.karney@sr...>  20140518 17:54:10

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 <charles.karney@...> 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 >> 
From: Edward Lam <e4lam@ya...>  20140518 16:33:27

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 <charles.karney@...> 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 > 
From: Charles Karney <charles.karney@sr...>  20140517 21:04:49

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 > 
From: Edward Lam <e4lam@ya...>  20140517 19:45:15

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 
From: Charles Karney <charles.karney@sr...>  20140110 20:37:33

I have this change checked in. It will be part of 1.35. On 20140110 14:47, Ben Adler wrote: > Great. Please tell me in case this fix will *NOT* make it into 1.35.  Charles Karney <charles.karney@...> SRI International, Princeton, NJ 085435300 Tel: +1 609 734 2312 Fax: +1 609 734 2662 
From: Ben Adler <adler@in...>  20140110 19:47:38

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 
From: Charles Karney <charles.karney@sr...>  20140110 18:47:34

On 20140110 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 (googlehit at > geospatialpython) myself? > > Comparing that file with the WKT from http://spatialreference.org, the > text file seems to be missing the AXISentries. 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 
From: Ben Adler <adler@in...>  20140110 16:58:34

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 (googlehit at geospatialpython) myself? Comparing that file with the WKT from http://spatialreference.org, the text file seems to be missing the AXISentries. 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 
From: Charles Karney <charles.karney@sr...>  20140110 13:50:45

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("egm965"); 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 reuse 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 > _______________________________________________ > Geographiclibusers mailing list > Geographiclibusers@... > https://lists.sourceforge.net/lists/listinfo/geographiclibusers > 
From: Ben Adler <adler@in...>  20140110 11:39:12

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 reuse 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 
From: Charles Karney <charles.karney@sr...>  20130424 11:39:18

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 
From: Shaun Walbridge <shaun.walbridge@gm...>  20121210 20:18:18

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 <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 > > >  > 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 valueadd services > Discover what IT Professionals Know. Rescue delivers > http://p.sf.net/sfu/logmein_12329d2d > _______________________________________________ > Geographiclibusers mailing list > Geographiclibusers@... > https://lists.sourceforge.net/lists/listinfo/geographiclibusers > 
From: hugh <hdixon@bi...>  20121210 12:19:35

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 
From: hugh <hdixon@bi...>  20121208 23:52:13

On Thu, 6 Dec 2012 11:09:47 +1100 hdixon<hdixon@...> wrote: > Hi, > I am trying to build a QT app using the latest libGeographic library. > I have used libGeoGraphic successfully with QT previously, but this > time I have a linker error which I cannot resolve. > I am wanting to convert eastings / northings to lat/ long using the > GeographicLib::GeoCoords class. > > > My linker output, from qt creator, is: > ************************************************************************* > g++ Wl,rpath,/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib o > crashStat2kml main.o xmlfile.o > L/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib > L/usr/local/lib/libGeographic lphonon > L/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib > lpulsemainloopglib lpulse lglib2.0 lQtDBus L/usr/X11R6/lib64 > lQtMultimedia lQtScriptTools lQtScript lQtSvg lQtGui > lQtXmlPatterns lQtNetwork lQtXml lQtCore lpthread main.o: In > function `GeographicLib::GeoCoords::Reset(int, bool, double, > double)': /usr/local/include/GeographicLib/GeoCoords.hpp:260: > undefined reference to `GeographicLib::UTMUPS::Reverse(int, bool, > double, double, double&, double&, double&, double&, bool)' make: > Leaving directory > `/home/hugh/code/vr/crashStat2kmlbuilddesktopDesktop_Qt_4_8_1_for_GCC__Qt_SDK__Debug' /usr/local/include/GeographicLib/GeoCoords.hpp:265: > undefined reference to `GeographicLib::GeoCoords::FixHemisphere()' > collect2: error: ld returned 1 exit status make: *** [crashStat2kml] > Error 1 10:10:51: The process "/usr/bin/make" exited with code 2. > Error while building project crashStat2kml (target: Desktop) > > ************************************************************************* > > I have had trouble in the past with order of parameters in the call to > the linker, and am thinking there is something in libGeographic and QT > that conflicts or upsets. > > Does anyone have any experiance using libG with QT? andif so have > they seen this problem and can you help please! > > With thanks, > H > >  > 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 valueadd > services Discover what IT Professionals Know. Rescue delivers > http://p.sf.net/sfu/logmein_12329d2d > _______________________________________________ > Geographiclibusers mailing list > Geographiclibusers@... > https://lists.sourceforge.net/lists/listinfo/geographiclibusers Problem found in my configuration, had nothing to do with libGeographics! Sorry! H 
From: hdixon<hdixon@bi...>  20121206 00:10:00

Hi, I am trying to build a QT app using the latest libGeographic library. I have used libGeoGraphic successfully with QT previously, but this time I have a linker error which I cannot resolve. I am wanting to convert eastings / northings to lat/ long using the GeographicLib::GeoCoords class. My linker output, from qt creator, is: ************************************************************************* g++ Wl,rpath,/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib o crashStat2kml main.o xmlfile.o L/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib L/usr/local/lib/libGeographic lphonon L/home/hugh/QtSDK1.2.1/Desktop/Qt/4.8.1/gcc/lib lpulsemainloopglib lpulse lglib2.0 lQtDBus L/usr/X11R6/lib64 lQtMultimedia lQtScriptTools lQtScript lQtSvg lQtGui lQtXmlPatterns lQtNetwork lQtXml lQtCore lpthread main.o: In function `GeographicLib::GeoCoords::Reset(int, bool, double, double)': /usr/local/include/GeographicLib/GeoCoords.hpp:260: undefined reference to `GeographicLib::UTMUPS::Reverse(int, bool, double, double, double&, double&, double&, double&, bool)' make: Leaving directory `/home/hugh/code/vr/crashStat2kmlbuilddesktopDesktop_Qt_4_8_1_for_GCC__Qt_SDK__Debug' /usr/local/include/GeographicLib/GeoCoords.hpp:265: undefined reference to `GeographicLib::GeoCoords::FixHemisphere()' collect2: error: ld returned 1 exit status make: *** [crashStat2kml] Error 1 10:10:51: The process "/usr/bin/make" exited with code 2. Error while building project crashStat2kml (target: Desktop) ************************************************************************* I have had trouble in the past with order of parameters in the call to the linker, and am thinking there is something in libGeographic and QT that conflicts or upsets. Does anyone have any experiance using libG with QT? andif so have they seen this problem and can you help please! With thanks, H 