On Sat, 28 Jan 2012, Agustin Lobo wrote:
Right, grid1km is the polygon object from whence cents comes from and
I copied the wrong line.
Sorry for the confusion.
The confusion is in proj4 in the ptransform() function, as:
centslonlat<- spTransform(cents,CRS("+init=epsg:4326"))
coordinates(centslonlat)[1:3,]
? ? ? ? ? ?X ? ? ? ?Y
[1,] 1.315311 42.23956
[2,] 1.327326 42.24060
[3,] 1.339342 42.24163
ptransform(data = coordinates(cents),
?src.proj = CRSargs(CRS("+init=epsg:3035")),
?dst.proj = "+proj=longlat +ellps=WGS84 +datum=WGS84"
?)[1:3,]*57.29577951308232
? ? ? ? [,1] ? ? [,2] [,3]
[1,] 1.315311 42.23956 ? ?0
[2,] 1.327326 42.24060 ? ?0
[3,] 1.339342 42.24163 ? ?0
ptransform() omits to test for geographical coordinates in the dst.proj=
argument, and consequently does not multiply by RAD_TO_DEG. The correct code
is in rgdal/src/projectit.cpp, from line 255. ptransform() is more recent
than project() in proj4, so the fact that it isn't the same (either source
or destination as longlat) got overlooked.
Hope this clarifies,
Roger