[Proj] Using pj_transform API to convert WGS84 to OSGB

[Keith Sharp]

On Fri, 2006-01-20 at 09:44 -0800, Eric Miller wrote:
> Don't multiply height values with DEG_TO_RAD.  Heights are in meters
> relative to the ellipsoid surface.  But, since you start with zero,
> that doesn't really impact anything.  If you don't have height values,
> pj_transform supports passing NULL for the Z-value pointer.

Ok, switched to using NULL instead of z as my test case doesn't include
height.

> Your k-factor appears to be slightly different between command line
> and the program.

That was me tinkering, I've made them both the same again.

> But, I think your biggest problem is in your pj_init call.  The first
> argument is supposed to be the number of elements in the array.  Your
> argument for osgb is only 2.  And none of those first two says don't
> use default values. 

Yep, this was a bug introduced by me randomly changing things with
little or no understanding in an attempt to make it work.  Even changing
the pj_init call of OSGB to use 9 as the element count it still gives
different answers:

[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

[kms at animal test]$ ./test 
Input ln = -4.134063e+00  Output X = -6.373681e+06  Expected X = 2.660000e+05
Input lt = 5.570597e+01  Output Y = 8.231288e+06  Expected Y = 6.479000e+05

I must be missing something fundamental, I've been looking at the cs2cs
source and I cannot see what it's main () and process() functions are
doing that I am not.

Keith.

Slightly fixed up 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.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 ln = -4.134063, lt = 55.705968;
	double x, y;
	int e = 0;
	
	double east = 266000;
	double nrth = 647900;

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

	if (!(wgs84PJ = pj_init (4, wgs84))) {
		printf ("Failed to initialise WGS84\n");
		return 1;
	}

	if (!(osgbPJ = pj_init (9, osgb))) {
		printf ("Failed to initialise OSGB\n");
		return 1;
	}

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

	return 0;
}