more doubts with spTransform()
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