On Wed, 22 Sep 2010, Agustin Lobo wrote:
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)
? ? ? ? ? ? ?_
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.
I think that this is the issue - there does not seem to be a portable way of
setting the +towgs84= values in the representation used by ESRI for the
shapefile *.prj file - "DATUM["D_unknown", ..." is what comes out, and that
is maybe what QGIS tries to respect. Have you tried other drivers?
Roger
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