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

[Roger Bivand]

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