[Proj] Using pj_transform API to convert WGS84 to OSGB
Keith Sharp
kms at passback.co.uk
Fri Jan 20 12:11:50 EST 2006
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;
}
More information about the Proj
mailing list