[R-sig-Geo] Problems with projection using spTransform (rgdal)

[Roger Bivand]

On Tue, 21 Aug 2007, Paul Hiemstra wrote:

> Dear list,
>
> When I perform the following transformation and transform back to the
> original CRS, the original coordinates are not returned:
>
> library(gstat)
> library(rgdal)
> data(meuse)
> coordinates(meuse) =~ x + y
> proj4string(meuse) = CRS("+init=EPSG:28992")

Should be: "+init=epsg:28992"

> print(bbox(meuse))
> meuse = spTransform(meuse, CRS("+proj=stere +lat_0=90 +lon_0=0.0
> +lat_ts=60.0 +a=6378.388 +b=6356.912 +x_0=0 +y_0=0"))

This is +ellps=intl, is it? I'm not sure that the +a and +b are in km, 
ought they to be in m? What about datum conversion?

This works for me:

library(rgdal)
data(meuse)
coordinates(meuse) =~ x + y
proj4string(meuse) = CRS("+init=epsg:28992")
bbox(meuse)
meuse1 = spTransform(meuse, CRS(paste("+proj=stere +lat_0=90",
   "+lon_0=0.0 +lat_ts=60.0 +a=6378388 +b=6356912 +x_0=0 +y_0=0")))
bbox(meuse1)
meuse2 <- spTransform(meuse1, CRS("+init=epsg:28992"))
bbox(meuse2)
all.equal(bbox(meuse), bbox(meuse2))

Roger


> print(bbox(meuse))
> meuse = spTransform(meuse, CRS("+init=EPSG:28992"))
> print(bbox(meuse))
>
> The same code using a transformation between EPSG:28992 (Dutch CRS) and
> EPSG:4326 (latlong) works fine.
>
> Does someone have an idea what causes this?
>
> Paul
>
> Info for my system.
> > sessionInfo()
> R version 2.5.1 (2007-06-27)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=Dutch_Netherlands.1252;LC_CTYPE=Dutch_Netherlands.1252;LC_MONETARY=Dutch_Netherlands.1252;LC_NUMERIC=C;LC_TIME=Dutch_Netherlands.1252
>
> attached base packages:
> [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"
> [7] "base"
>
> other attached packages:
>   gstat    rgdal       sp
> "0.9-39" "0.5-15" "0.9-14"
>
>

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