[Proj] Using pj_transform API to convert WGS84 to OSGB

[Eric Miller]

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.

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

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. 

Let this compile time idiom count your array elements for you:

size_t els = sizeof(array)/sizeof(array[0]);

Obviously, it only works when the definition of the array is in scope.
 
>>> kms at passback.co.uk 1/20/2006 9:11 AM >>>
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;
}

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