[Proj] How to convert a sphere to ellipsoid with correct datum?

[Jan Hartmann]

  On 09/01/10 17:37, Hermann Peifer wrote:
> On 01/09/2010 12:54, Jan Hartmann wrote:
>> The same location on the
>> world globe has *different* latlong values when computed on the basis of
>> different ellipsoids, and that is not shown by cs2cs.
> I am not so sure about different *long* values. Anyway, can't you just
> convince cs2cs to do the calculation by putting +towgs84=0,0,0 on both
> sides of the equation?
>
> $ cs2cs -v +proj=latlong +ellps=sphere +a=6371000 +towgs84=0,0,0 +to
> +proj=latlong +ellps=WGS84 +towgs84=0,0,0
>
Thanks for your reaction. That's what I thought too, but the values are 
wrong. Around 1850, the latlon values of about 1000 reference points 
were computed on the basis of a slightly smaller ellipsoid than is used 
nowadays. The modern latlon values for the same points (on the Bessel 
ellipsoid) are about 50 meters to the south. How can I use cs2cs to 
transform the old value to the new?

1) using latlon for input and output with different ellipsoids:
cs2cs -f "%8.6f" +proj=longlat +a=6376950.4 +rf=309.65 +to +proj=longlat 
+ellps=bessel
4.8838828 51.5
4.883883        51.500000 0.000000

No change

2) Adding +towgs84=0,0,0 to input and output:
cs2cs -f "%8.6f" +proj=longlat +a=6376950.4 +rf=309.65 +towgs84=0,0,0 
+to +proj=longlat +ellps=bessel +towgs84=0,0,0
4.8838828 51.5
4.883883        51.506347 -3.768025

Using PostGIS to compute the distance:
  select st_distance_spheroid(st_geomfromtext('POINT(4.883883 
51.5)',4308),st_geomfromtext('POINT(4.883883 51.506347)',4308), 
'SPHEROID["K",6376950.4,309.65]');
  st_distance_spheroid
----------------------
      706.037082633133
(1 row)

Distance far too great.

The only way I can get from the old latlon value to the new is using 
intermediate projections: from old latlon to old projected coordinates, 
from old projected coordinates to new projected coordinates using the 
old ellipsoid, and from new projected coordinates back to latlon on the 
new ellipsoid. The results look OK (within 10 meters), but I don't 
understand what I am doing here. I certainly dislike the faking of an 
old ellipsoid on a new projection.

It's important, because all pre-1900 maps have to be georeferenced this 
way, and I need to have accuracy on meter level. I still think 
transformation 1) should do the job

Jan

BTW Both latitude and longitude are different in the 1850 and modern 
tables, although the differences in latitudes are much larger.