[R-sig-Geo] more doubts with spTransform()

[Roger Bivand]

On Tue, 14 Jul 2009, Agustin Lobo wrote:

> Inverse transform now:
>> anyellaWGS84 <- readOGR(dsn=".", layer="parcellesAnyellaWGS84")
> OGR data source with driver: ESRI Shapefile
> Source: ".", layer: "parcellesAnyellaWGS84"
> with  6  rows and  7  columns
> Feature type: wkbPoint with 2 dimensions
>> proj4string(anyellaWGS84)
> [1] " +proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +no_defs 
> +towgs84=0,0,0"
>> anyellaED50notowgs84 <- spTransform(anyellaWGS84,CRS("+proj=utm +zone=31 
> +ellps=intl +units=m +no_defs"))
>> anyellaED50towgs84 <- spTransform(anyellaWGS84,CRS("+proj=utm +zone=31 
> +ellps=intl +units=m +no_defs +towgs84=-87,-98,-121"))
>> coordinates(anyellaWGS84)[1,]
> coords.x1 coords.x2
>  421005   4683849
>> coordinates(anyellaED50notowgs84)[1,]
> coords.x1 coords.x2
> 421001.4 4683932.4
>> coordinates(anyellaED50towgs84)[1,]
> coords.x1 coords.x2
> 421097.5 4684050.9
>
> which is correct? anyellaED50notowgs84 or anyellaED50towgs84 ?
>
> I guess anyellaED50towgs84, but would like to confirm,

Good guess! In the no +towgs84 case, the ellipsoid is being used but not 
datum transformation, which leaves a contradiction between a declared 
ellipsoid and the one inherent in the default datum (WGS84). Providing the 
+towgs84 transforms to the datum expressed in the given parameters, which 
is one of the ED50 family.

The difference is anout 100m in each direction. Maybe try this on an 
easily identifiable known ground control point where you have positional 
data from several maps, like the corner of a building or a bridge pier.

Roger

>
> Thanks
>
> Agus
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no