[Proj] Verification of inverse step in projection

[Roger Oberholtzer]

On Dec 13, 2009, at 3:25 PM, José Luis García Pallero wrote:

> 2009/12/13 Roger Oberholtzer <roger at opq.se>:
>> Perhaps:
>>
>> projPJ *pjStruct;
>>
>> as pj_init_plus returns a pointer.
>>
>>
>>> pjStruct = pj_init_plus("+proj=tcc +ellps=clrk66 +lon_0=90W");
>>> if(!pjStruct)
>>>    ....
>>> if(pjStruct->inv==0)
>>>    ....
>
> But pj_init_plus() can't know if you will use the forward or the  
> inverse step.

These are two different issues:

1) You have a syntax error, which I tried to explain. That must be  
corrected no matter.

2) The function that does the actual transform is pj_transform. It  
takes to projPJ pointers. One describes the source data, and one  
describes the desired output. Swapping these two 'reverse' the  
transform.  Here is an example, with three functions. They are part of  
some larger system, but they should be clear enough to demonstrate the  
point. In my case, they  are for wgs84 <-> ed50.

int Init_wgs84_to_ed50(WGS84Info *tinfo) {

         // From GPS WGS84
         if (!((tinfo->fromProj) = pj_init_plus(
                 "+proj=latlong "
                 "+datum=WGS84 "
                 ))) return(ISERROR);

         // To ED50
         if (!((tinfo->toProj) = pj_init_plus(
                 "+proj=utm "
                 "+ellps=intl "
                 "+k=0.999600 "
                 "+x_0=500000.0 "
                 "+y_0=0.0 "
                 "+zone=30 "
                 "+units=m "
                  
"+towgs84=-125.098545,-76.000054,-156.198703,0.0,0.0,-1.129,8.30463103 "
                 "+no_defs"
                 ))) return(ISERROR);

         return(ISOK);
}

// Convert longitude/latitude to northing/easting.

int wgs84LL2NE( WGS84Info *tinfo,
                 double LAT,
                 double LONG,
                 double HEIGHT,
                 double *easting,
                 double *northing,
                 double *altitude) {

         double  lng = LONG * DEGREE_TO_RADIAN,
                 lat = LAT  * DEGREE_TO_RADIAN,
                 alt = HEIGHT;

         if (pj_transform(tinfo->fromProj,
                          tinfo->toProj,
                          1, 0,
                          &lng, &lat, &alt) != ISOK) return(ISERROR);

         *easting  = lng;
         *northing = lat;
         *altitude = alt;

         return(ISOK);
}

// Convert northing/easting to longitude/latitude.

int wgs84NE2LL( WGS84Info *tinfo,
                 double easting,
                 double northing,
                 double altitude,
                 double *LAT,
                 double *LONG,
                 double *HEIGHT) {

         double  le = easting,
                 ln = northing,
                 la = altitude;

         if (pj_transform(tinfo->toProj,
                          tinfo->fromProj,
                          1, 0,
                          &le, &ln, &la) != ISOK) return(ISERROR);

         *LONG   = le / DEGREE_TO_RADIAN;
         *LAT    = ln / DEGREE_TO_RADIAN;
         *HEIGHT = la / DEGREE_TO_RADIAN;

         return(ISOK);
}

Init_wgs84_to_ed50 - called once to initialize the transforms.

wgs84LL2NE - called with a lat/long in WGS84, returns northing/easting  
in ED50

wgs84NE2LL - called with northing/easting in ED50, returns lat/long in  
WGS84

The last two functions are used no matter what the transform is. It is  
the first function that defines the WGS84 lat/long, and the ed50  
northing/easting.

I was surprised how easy it was to get the reverse one working. Just  
reverse the first two parameters to pj_transform. I am sure there are  
places where there are problems. But for all the transforms I go  
(including OSTN02 with a grid shift), it just works.


Hope this helps. The secret is in the two calls to pj_init_plus().  
Note that both are needed even if you're are going in one direction  
only. One describes what you have, the other what you want.

>

Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Kapellgränd 7
P.O. Box 4205
SE-102 65 Stockholm, Sweden

Office: Int +46 8-615 60 20
Mobile: Int +46 70-815 1696

And remember:

It is RSofT and there is always something under construction.
It is like talking about large city with all constructions finished.
Not impossible, but very unlikely.