[R-sig-Geo] Problem with spTransform en +towgs84 ?

[Agustin Lobo]

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