Problem with spTransform en +towgs84 ?

Roger,

I apologize for not including version details etc, here they are:
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.2, released 2010/04/23
Path to GDAL shared files: /usr/share/gdal17
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)

version

              _
platform       i486-pc-linux-gnu
arch           i486
os             linux-gnu
system         i486, linux-gnu
status
major          2
minor          11.1
year           2010
month          05
day            31
svn rev        52157
language       R
version.string R version 2.11.1 (2010-05-31)

Also, I've put SPDF polsfl1_1000 here:
https://sites.google.com/site/openfiles2/home/polsfl1_1000.rda?attredirects=0&d=1

Finally, the towgs84 parameters come from the last paragraph in the pdf 
for The Kingdom of Spain (July) listed in the page you mention. (BTW, I 
recall having seen another list with slightly different parameters but 
cannot find it anymore) I had to put the same 3 parameters for the 
definition of the projection in QGIS (which qgis does not do by default, 
had to define a custom projection to get the datum shift working). The 
towgs parameters are not in the proj file of the shapefile set, at least 
not in the ones created by writeOGR.

Hope this helps to clarify the problem. I'm still not certain that
spTransform is doing anything wrong.

Thanks for your help,

Agus


2010/9/21 Roger Bivand <Roger.Bivand at nhh.no>:

On Tue, 21 Sep 2010, Agustin Lobo wrote:

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

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