[Proj] Using pj_transform API to convert WGS84 to OSGB

[Roger Oberholtzer]

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