[Proj] Using pj_transform API to convert WGS84 to OSGB

[Keith Sharp]

On Fri, 2006-01-20 at 08:51 -0500, Frank Warmerdam wrote:
> 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.

My celebrations were a little premature :-(

I am still not getting the same values from the command:

[kms at animal include]$ /home/kms/tmp/install/bin/cs2cs +proj=latlong +ellps=WGS84 +towgs84=0,0,0 +nodefs +to +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
-4.134063 55.705968
266000.11       647899.98 -53.12

As from my new (slightly improved) C program:

[kms at animal test]$ ./test
Input ln = -4.134063e+00  Output X = -6.234896e+06  Expected X = 266000
Input lt = 5.570597e+01  Output Y = 7.443438e+06  Expected Y = 647900

The results from the command line match the output from the conversion
routine used by Streetmap:

http://www.streetmap.co.uk/streetmap.dll?GridConvert?name=266000,647900&type=OSGrid

I am still missing something, any suggestions?

Thanks,

Keith.

Latest source:

#include <stdio.h>
#include <proj_api.h>

/*

	Test point is Glasgow Airport:
	
		Landranger: 	NS 479 660 
		WGS84 Lat:		55.705968
		WGS84 Lon:		-4.134063
		OS X:			266000
		OS Y:			647900

*/

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.999601272000000040", "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 ln = -4.134063, lt = 55.705968, at = 0;
	double x, y, z;
	int e = 0;

	x = ln * DEG_TO_RAD;
	y = lt * DEG_TO_RAD;
	z = at * DEG_TO_RAD;

	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;
	}

	if ((e = pj_transform (wgs84PJ, osgbPJ, 1, 0, &y, &x, &z)) != 0) {
		printf ("Transform failed: %s\n", pj_strerrno (e));
		return e;
	}
	printf ("Input ln = %e  Output X = %e  Expected X = 266000\n", ln, x);
	printf ("Input lt = %e  Output Y = %e  Expected Y = 647900\n", lt, y);

	return 0;
}