Problem with spTransform en +towgs84 ?

Hi!

Might it be that spTransform is not taking the +towgs84 string into account?
I have:
proj4string(polsfl1_1000)
[1] " +proj=utm +zone=31 +ellps=intl +units=m +no_defs +towgs84=-87,-98,-121"

and do:
polsfl1_1000lonlat <- spTransform(polsfl1_1000, CRS("+proj=longlat +datum=WGS84"))
writeOGR(polsfl1_1000,  dsn="/media/TRANSCEND/FLUXPYR/FLIGHT201009/polsfl1_1000",layer="polsfl1_1000", driver="ESRI Shapefile")
writeOGR(polsfl1_1000lonlat,  dsn="/media/TRANSCEND/FLUXPYR/FLIGHT201009/polsfl1_1000lonlat.kml",layer="polsfl1_1000", driver="KML")
Then, in QGIS, with reprojection on the fly activated, the 2 vectors
are not coincident.
Instead, a polsfl1_1000lonlatv2 made from the original polsfl1_1000
shapefile in QGIS, is coincident.

As there are 2 different programs here, I'm not positive of my diagnostic. But
I think there is a problem with spTransform
Could you please provide links to screen shots and input and output data 
objects, as well as sessionInfo() output? Which OS is involved, and are 
spTransform() and QGIS using the same PROJ.4 library and metadata. Is your 
input actually:

CRS("+init=epsg:23031 +towgs84=-84,-107,-120")

with the three-parameter datum shift from:

http://www.asprs.org/resources/GRIDS/

of July 2000? Where does yours come from? The way QGIS treats the *.prj 
file is also of interest.

Roger
Thanks

Agus

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

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

Problem with spTransform en +towgs84 ?

Thread (5 messages)