[Proj] Using pj_transform API to convert WGS84 to OSGB
Eric Miller
EMiller at dfg.ca.gov
Fri Jan 20 12:44:04 EST 2006
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
More information about the Proj
mailing list