[Proj] Using pj_transform API to convert WGS84 to OSGB

[Frank Warmerdam]

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