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

Jan Hartmann j.l.h.hartmann at uva.nl
Wed Sep 1 13:06:32 EST 2010

```  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.
```