[Proj] Using pj_transform API to convert WGS84 to OSGB
Frank Warmerdam
warmerdam at pobox.com
Fri Jan 20 08:51:18 EST 2006
On 1/20/06, Keith Sharp <kms at passback.co.uk> wrote:
> int
> main (int argc, char *argv[])
> {
> char *wgs84[] = { "proj=latlong", "ellps=WGS84", "towgs84=0,0,0", "nodef" };
>
> char *osgb[] = { "proj=tmerc", "lat_0=49", "lon_0=-2", "k=0.9996012717", "x_0=400000", "y_0=-100000", "ellps=airy", "towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", "units=m" };
>
> projPJ wgs84PJ;
> projPJ osgbPJ;
>
> double x = 55.863157, y = -4.432298, z = 0;
> int e = 0;
>
> if (!(wgs84PJ = pj_init (4, wgs84))) {
> printf ("Failed to initialise WGS84\n");
> return 1;
> }
>
> if (!(osgbPJ = pj_init (2, osgb))) {
> printf ("Failed to initialise OSGB\n");
> return 1;
> }
>
> printf ("Converting from WGS84 to OSGB\n");
> printf ("X: %e Y: %e Z: %e\n", x, y, z);
> if ((e = pj_transform (wgs84PJ, osgbPJ, 1, 0, &y, &x, &z)) != 0) {
> printf ("Transform failed: %s\n", pj_strerrno (e));
> return e;
> }
> printf ("X: %e Y: %e Z: %e\n", x, y, z);
> printf ("Expected result was X: 247878 and Y: 666006\n");
>
> return 0;
> }
Keith,
As mentioned, you have the x/y mixed up. But as
importantly, the pj_transform() takes geographic locations
in radians, not degrees like the user accessable interface
in cs2cs.
Multiply your values by DEG_TO_RAD to convert from degrees
to radians.
Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent
More information about the Proj
mailing list