[Proj] ED50 from C API

[Roger Oberholtzer]

On Thu, 2007-03-29 at 09:05 -0700, Eric Miller wrote:
> I don't know about the "esri" file, but you have two different ellipsoids and no information to convert between the different models.  Don't you need a datum or towwg84 parameter for the ED50 definition?  I think the ED50 towgs84 parameters are missing because there is more than one candidate.  See, http://crs.bkg.bund.de/crs-eu/
> 
> There isn't really much difference between using cs2cs and pj_transform.  The init strings are the same.  Same for proj and pj_fwd/pj_inv.

I have since tried a few ED50 definitions to see which is closest to the
target locations:

"+towgs84=-87,-98,-121 "     // "EUROPEAN 1950, Mean"

"+towgs84=-87,-96,-120 "     // "EUROPEAN 1950, Western Europe"

"+towgs84=-84,-107,-120 "       // "EUROPEAN 1950, Portugal &
Spain" (best)

"+towgs84=-86,-98,-119 "     // "EUROPEAN 1979"

The best is "+towgs84=-84,-107,-120 ", which gets me to within 5 meters.
Still too far, but better.

It still remains a mystery to me how anyone actually uses the esri file
when it is quite often incomplete. To the point it does not work at all.
Especially for a conversion as 'big' as one used for all roads in Spain.

Perhaps proj needs its own file that contains complete working
definitions? proj is a great library. But without proper definitions is
is less useful than it could be. I understand the complexity of the
thing since so many things can be converted into so many things,
resulting in a zillion combinations. But there must be something
possible. Or is the solution staring me in the face and I just don't see
it? Do all users of proj make all the definitions locally (read
"reinvent the wheel")?

> >>> On 3/29/2007 at 3:40 AM, Roger Oberholtzer <roger at opq.se> wrote:
> > Let me ask this another way:
> > 
> > proj has definitions for many conversions in the 'esri' file that that
> > it provides. These definitions appear to be used as changes from some
> > default initial settings. I say that as many say very little. The rest
> > must be specified elsewhere. When I want to add a new conversion via the
> > C API, I expect I could start by looking for an existing definition in
> > the 'esri' file. As the list can probably tell, each time I find I need
> > to add a new conversion, I never get it right and then ask the list.
> > Each time the corrections seem (to me) to be different. So, I am missing
> > some important part of the puzzle. Does no one else on the list use the
> > C API very much?
> > 
> > On Wed, 2007-03-28 at 16:04 +0200, Roger Oberholtzer wrote:
> >> I would like to convert WGS84 lat/longs from a GPS receiver into
> >> northings and eastings in ED50, fixed in UTM zone 30N. I have tried the
> >> following, but the results are always 100-150 meters off from what I
> >> expect. What have I missed? I based this in the ESRI file in proj 4.4.9
> >> 
> >> Initialization (error checking removed here for clarity):
> >> 
> >>         // From GPS WGS84
> >>         fromProj = pj_init_plus( "+proj=latlong "
> >>                                                 "+ellps=WGS84 "
> >>                                                 "+towgs84=0,0,0 "
> >>                                                 "+no_defs"):
> >> 
> >>         // To ED50
> >>         toProj   = pj_init_plus( "+proj=utm "
> >>                                                 "+ellps=intl "
> >>                                                 "+zone=30 "
> >>                                                 "+units=m "
> >>                                                 "+no_defs"):
> >> 
> >> The call for each location:
> >> 
> >>         lng = LONG * DEGREE_TO_RADIAN;
> >>         lat = LAT  * DEGREE_TO_RADIAN;
> >>         alt = 0; // We do not give this a height
> >> 
> >>         pj_transform(fromProj, toProj, 1, 0, &lng, &lat, &alt);
> >> 
> >> For example: 
> >> 
> >> Input WGS84: lat=37.57508059 lng=-0.93032997 alt=0.00000000
> >> 
> >> Obtains: Easting=682771.459 Northing=4160841.260
> >> 
> >> But I am expecting: Easting=682876.090 Northing=4160892.020 
> >> 
> >> There is a difference in location of 116.293 meters.
> >> 
> >> The error is in the initialization step. I suspect it is something that
> >> is not in the ESRI file that I need to add when doing this from the C
> >> API. But what?
> >> 
> > -- 
> > Roger Oberholtzer
> > 
> > OPQ Systems / Ramböll RST
> > 
> > Ramböll Sverige AB
> > Kapellgränd 7
> > P.O. Box 4205
> > SE-102 65 Stockholm, Sweden
> > 
> > Tel: Int +46 8-615 60 20
> > Fax: Int +46 8-31 42 23
> > 
> > _______________________________________________
> > Proj mailing list
> > Proj at lists.maptools.org 
> > http://lists.maptools.org/mailman/listinfo/proj
> 
> 
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
-- 
Roger Oberholtzer

OPQ Systems / Ramböll RST

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

Tel: Int +46 8-615 60 20
Fax: Int +46 8-31 42 23