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

[jgarcia at ija.csic.es]

Hi,
Just to add that I've also detected some problems in this sense. I've have
exported a shapefile from ArcGIS to KML, and also imported it into R, and
from there to KML. The final KML files differ in some 70m of latitude,
while the longitudes match.

# import the shapefile
> knapdale.supgeo <- readOGR(dsn="test",layer="knapdale_supgeo")
> proj4string(knapdale.supgeo)
[1] " +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000
+ellps=airy +units=m +no_defs"

# Note, the above is OSGB 1936 / British National Grid, which  equals
# CRS("+init=epsg:27700")
# Then:

knapdale.supgeoLL <- spTransform(knapdale.supgeo, CRS("+proj=longlat +
datum=WGS84"))
kmlf <- "knapdale_supgeo.kml"
writeOGR(knapdale.supgeoLL, dsn=kmlf, layer="supgeo", driver="KML")

I'm not saying that the R part is responsible for the difference, but the
mismatch is there. Some details:

ArcGIS 9.3 (Windows XP Pro)
R 2.10.1
> library(rgdal)
Loading required package: sp
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.4.2.0, released 2007/06/27
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.5.0, 22 Oct 2006
Path to PROJ.4 shared files: (autodetected)

Best regards,
Javier



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