[Proj] Verification of inverse step in projection
Roger Oberholtzer
roger at opq.se
Sun Dec 13 12:03:09 EST 2009
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.
More information about the Proj
mailing list