Charles Karney | 7 Dec 2011 15:10

Earth gravity (and magnetic) models in GeographicLib

I've added earth gravity models to GeographicLib.  This allows you to
compute quantities like the gravitational field, the deflections of the
vertical, and the geoid height directly from the EGM2008 (or EGM96 or
EGM84) model coefficients.  For details see

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

My development of this capability in GeographicLib relied, in part, on
the code for EGM_Geoid_CalculatorClass by Mathieu Peyréga.

GeographicLib is able to reproduce accurately the gridded data for
EGM2008 published by the NGA.  Single point calculations for EGM2008 are
rather slow (78 ms for the geoid height).  However the computations can
be sped up by a factor of 800 by computing several points on a circle of
latitude and by using, for example, OpenMP, for parallel processing.
Thus, on my 8-processor desktop machine I can compute the global 1' grid
of geoid heights in 40 mins (i.e., about 10 us per point).

P.S. Magnetic models are also supported, see

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

--

-- 
Charles Karney <charles.karney <at> sri.com>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662
_______________________________________________
Proj mailing list
(Continue reading)

Cyril Bondue | 7 Dec 2011 21:58
Picon
Favicon

Proj4js, cs2cs and pj_transform

Hello fellow coders!

I'm writing to ask for your help. I need to convert data from NTF (Paris) / Lambert Centre France (epsg 27562) to WGS84 (epsg 4326). To test Proj4 I used the web interface Proj4js, and i could do the conversion without much troubles. Here are the two projPj I used, and working :
 - source : +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs
 - dest : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
Good! But I need to do the conversion in a C++ program, so I downloaded the source, downloaded the NAD for France (who were already in the archive BTW), compiled with MinGW (i'm doing a Windows software), compilation runs without problems, great!

Now, troubles ^^
Here is my code :

projPJ source = pj_init_plus("+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs");
projPJ dest = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs");

double x = 535003.285368d;
double y = 200358.358412d;
pj_transform(source, dest, 1, 1, &x, &y, NULL);

pj_transform returns 0 (OK), but x and y are wrong. I have x=0.0259181 and y=0.816814 while the correct results are 1.485 and 46.8 (rounded up), results that Proj4Js confirms.
I also tried the same transformation on cs2cs :
cs2cs.exe +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356 515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
535003.285368 200358.358412
Results are different and also wrong (or i can't read them right) x=1d29'6"E, y=46d48'N 44.005

I have taken care of settings PROJ_LIB with the nad folder. So i must do something wrong, i've tried many things but i don't understand what happens, so here I am, asking you for some assistance. Can you see what's going on?

Thank you in advance!

Cyril
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
nicolas david | 7 Dec 2011 22:57
Picon

Re: Proj4js, cs2cs and pj_transform

Hi

pj_transform use lat/long coordinate in radian. So you have to convert them to degree with  "RAD_TO_DEG "  defined at the value  57.29577951308232 in proj_api.h

Nicolas

2011/12/7 Cyril Bondue <c.bondue <at> hotmail.com>
Hello fellow coders!

I'm writing to ask for your help. I need to convert data from NTF (Paris) / Lambert Centre France (epsg 27562) to WGS84 (epsg 4326). To test Proj4 I used the web interface Proj4js, and i could do the conversion without much troubles. Here are the two projPj I used, and working :
 - source : +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs
 - dest : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
Good! But I need to do the conversion in a C++ program, so I downloaded the source, downloaded the NAD for France (who were already in the archive BTW), compiled with MinGW (i'm doing a Windows software), compilation runs without problems, great!

Now, troubles ^^
Here is my code :

projPJ source = pj_init_plus("+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs");
projPJ dest = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs");

double x = 535003.285368d;
double y = 200358.358412d;
pj_transform(source, dest, 1, 1, &x, &y, NULL);

pj_transform returns 0 (OK), but x and y are wrong. I have x=0.0259181 and y=0.816814 while the correct results are 1.485 and 46.8 (rounded up), results that Proj4Js confirms.
I also tried the same transformation on cs2cs :
cs2cs.exe +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356 515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
535003.285368 200358.358412
Results are different and also wrong (or i can't read them right) x=1d29'6"E, y=46d48'N 44.005

I have taken care of settings PROJ_LIB with the nad folder. So i must do something wrong, i've tried many things but i don't understand what happens, so here I am, asking you for some assistance. Can you see what's going on?

Thank you in advance!

Cyril

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Jean-Claude Repetto | 8 Dec 2011 09:35
Picon
Favicon

Re: Proj4js, cs2cs and pj_transform

Le 07/12/2011 21:58, Cyril Bondue a écrit :

> Correct results are 1.485 and 46.8 (rounded
> up), results that Proj4Js confirms.

> I also tried the same transformation on cs2cs :
> cs2cs.exe +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
> +x_0=600000 +y_0=200000 +a=6378249.2 +b=6356 515
> +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs +to
> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
> 535003.285368 200358.358412
> Results are different and also wrong (or i can't read them right)
> x=1d29'6"E, y=46d48'N 44.005"

No, these results are correct.

Jean-Claude

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

Roger Oberholtzer | 8 Dec 2011 12:37
Picon
Gravatar

Projection advice

I have a system that starts with WGS84 lat/long, projects them so as to
obtain northings/eastings, does some fiddling with the locations, and
then converts them back in to WGS84 lat/long. The conversion from/to
lat/long is done by proj. The fiddling with the northings/eastings is
done based on some INS that provides relative movement information in
meters. Only the WGS84 lat/longs are kept.

Currently, I use a projection for the local region. The thought was that
this would best 'match' the INS movements, which are of course always
local. I would like to generalize this solution by perhaps not using a
local projection. Instead use a generic projection everywhere. It seems
reasonable to do this as (1) the projection will not be kept, and (2)
the movements applied to the projections are relative.

My question is; is there a 'best' projection to use for such a
procedure?

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Noel Zinn (cc | 8 Dec 2011 13:41
Picon

Re: Projection advice

Roger,

Presumably you want a conformal projection so that the scale distortion is 
uniform in all directions.  Does the INS fiddling in Northings/Easting 
account for scale distortion (point scale factor) at all, or does it presume 
there is no scale distortion?  If the latter, then the larger the area your 
generic projection covers, the greater the scale distortion your INS 
fiddling will ignore.  That may not be a good change from your current, 
local-projection procedure.

Noel

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn <at> hydrometronics.com (email)
http://www.hydrometronics.com (website)

-----Original Message----- 
From: Roger Oberholtzer
Sent: Thursday, December 08, 2011 5:37 AM
To: proj <at> lists.maptools.org
Subject: [Proj] Projection advice

I have a system that starts with WGS84 lat/long, projects them so as to
obtain northings/eastings, does some fiddling with the locations, and
then converts them back in to WGS84 lat/long. The conversion from/to
lat/long is done by proj. The fiddling with the northings/eastings is
done based on some INS that provides relative movement information in
meters. Only the WGS84 lat/longs are kept.

Currently, I use a projection for the local region. The thought was that
this would best 'match' the INS movements, which are of course always
local. I would like to generalize this solution by perhaps not using a
local projection. Instead use a generic projection everywhere. It seems
reasonable to do this as (1) the projection will not be kept, and (2)
the movements applied to the projections are relative.

My question is; is there a 'best' projection to use for such a
procedure?

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj 

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Roger Oberholtzer | 8 Dec 2011 16:10
Picon
Gravatar

Re: Projection advice

On Thu, 2011-12-08 at 06:41 -0600, Noel Zinn (cc) wrote:
> Roger,
> 
> Presumably you want a conformal projection so that the scale distortion is 
> uniform in all directions.  

Yes.

> Does the INS fiddling in Northings/Easting 
> account for scale distortion (point scale factor) at all, or does it presume 
> there is no scale distortion?  If the latter, then the larger the area your 
> generic projection covers, the greater the scale distortion your INS 
> fiddling will ignore.  That may not be a good change from your current, 
> local-projection procedure.

The INS is comprised of 3D gyroscopes, multi-axis accelerometers, a
wheel pulse transducer, and in some cases inclinometers. From these is
obtained very accurate relative movement information. Drift is minimal
over the distances it is expected to be used. These devices are scaled
and offset against a well-defined surface in the general area (within
the same country...) as they are used. Gravitational and surface curve
effects are minimal over the short distances they are used. At least
that is our working model.

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Noel Zinn (cc | 8 Dec 2011 19:28
Picon

Re: Projection advice

GPS, too?  Anyway, I take your response as an affirmative (i.e. scale 
accounted for).  Nevertheless, I'd advise sticking with local (and 
disposable) projections, e.g. an oblique stereographic centered in your 
(small) work area with (probably) a scale of unity at the center.  The 
"best" wide-area projection may not be good enough for some (or many) work 
areas.  -Noel

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn <at> hydrometronics.com (email)
http://www.hydrometronics.com (website)

-----Original Message----- 
From: Roger Oberholtzer
Sent: Thursday, December 08, 2011 9:10 AM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Projection advice

On Thu, 2011-12-08 at 06:41 -0600, Noel Zinn (cc) wrote:
> Roger,
>
> Presumably you want a conformal projection so that the scale distortion is
> uniform in all directions.

Yes.

> Does the INS fiddling in Northings/Easting
> account for scale distortion (point scale factor) at all, or does it 
> presume
> there is no scale distortion?  If the latter, then the larger the area 
> your
> generic projection covers, the greater the scale distortion your INS
> fiddling will ignore.  That may not be a good change from your current,
> local-projection procedure.

The INS is comprised of 3D gyroscopes, multi-axis accelerometers, a
wheel pulse transducer, and in some cases inclinometers. From these is
obtained very accurate relative movement information. Drift is minimal
over the distances it is expected to be used. These devices are scaled
and offset against a well-defined surface in the general area (within
the same country...) as they are used. Gravitational and surface curve
effects are minimal over the short distances they are used. At least
that is our working model.

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se 

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Roger Oberholtzer | 8 Dec 2011 21:33
Picon
Gravatar

Re: Projection advice

On Thu, 2011-12-08 at 12:28 -0600, Noel Zinn (cc) wrote:
> GPS, too?  Anyway, I take your response as an affirmative (i.e. scale 
> accounted for).  Nevertheless, I'd advise sticking with local (and 
> disposable) projections, e.g. an oblique stereographic centered in your 
> (small) work area with (probably) a scale of unity at the center.  The 
> "best" wide-area projection may not be good enough for some (or many) work 
> areas.  -Noel

GPS is always the tricky part. When is it really good? A Kahlman filter
is the best practical solution. For that we use systems like those from
Oxford Technical Solutions. In the case I am discussing here, we are
messing around a bit. 

But how to make this projection choice automatic? That is what I am
trying to accomplish.

-- 

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Noel Zinn (cc | 8 Dec 2011 22:28
Picon

Re: Projection advice

> But how to make this projection choice automatic?

GPS would solve that, of course.  Once you have a position in your area, you 
define that as the central point of an oblique stereographic, for example. 
If no GPS, then perhaps an atlas!  BTW, thanks for the OxTS lead.  -Noel

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn <at> hydrometronics.com (email)
http://www.hydrometronics.com (website)

-----Original Message----- 
From: Roger Oberholtzer
Sent: Thursday, December 08, 2011 2:33 PM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Projection advice

On Thu, 2011-12-08 at 12:28 -0600, Noel Zinn (cc) wrote:
> GPS, too?  Anyway, I take your response as an affirmative (i.e. scale
> accounted for).  Nevertheless, I'd advise sticking with local (and
> disposable) projections, e.g. an oblique stereographic centered in your
> (small) work area with (probably) a scale of unity at the center.  The
> "best" wide-area projection may not be good enough for some (or many) work
> areas.  -Noel

GPS is always the tricky part. When is it really good? A Kahlman filter
is the best practical solution. For that we use systems like those from
Oxford Technical Solutions. In the case I am discussing here, we are
messing around a bit.

But how to make this projection choice automatic? That is what I am
trying to accomplish.

-- 

Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer <at> ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj 

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

Gmane