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

[Agustin Lobo]

I'm getting conflicting results on this issue. I will do more
and more formal testing and will report back.

Agus

2010/9/23 Roger Bivand <Roger.Bivand at nhh.no>:
> 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.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
>>>
>>
>
> --
> 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
>
>