[Proj] Using pj_transform API to convert WGS84 to OSGB
Roger Oberholtzer
roger at opq.se
Sun Jan 22 05:56:51 EST 2006
On Fri, 2006-01-20 at 18:26 +0000, Keith Sharp wrote:
> 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.
Here is what I use. The _plus functions make this easier. They take one
string. So there are NO commas between my strings. I just like to make
the call easier to read.
if (!fromProj = pj_init_plus(
"+proj=latlong "
"+ellps=WGS84 "
"+towgs84=0,0,0 "
"+nodefs"))) return(ISERROR);
if (!(toProj = pj_init_plus(
"+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 "
"+no_defs"))) return(ISERROR);
Then, setting LONG and LAT to the location in dd.mmmmm, I use it with:
lng = LONG * DEGREE_TO_RADIAN;
lat = LAT * DEGREE_TO_RADIAN;
alt = 0; // We do not give this a height
if (pj_transform(&fromProj, &toProj, 1, 0, &lng, &lat, &alt)
!= ISOK) return(ISERROR);
*localx = lng;
*localy = lat;
*localz = alt; // Whatever...
--
Roger Oberholtzer
More information about the Proj
mailing list