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)
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.
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.biv...@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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
