[Proj] Using pj_transform API to convert WGS84 to OSGB
Eric Miller
EMiller at dfg.ca.gov
Fri Jan 20 10:42:13 PST 2006
Your "y" and "x" arguments are in the wrong order in the pj_transform call. I didn't notice it before.
>>> kms at passback.co.uk 1/20/2006 10:26 AM >>>
On Fri, 2006-01-20 at 09:44 -0800, Eric Miller wrote:
> 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;
}
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
More information about the Proj
mailing list