[R-sig-Geo] more doubts with spTransform()

[Agustin Lobo]

Confirmed using Global Mapper, one of the
few packages performing correct on-the-fly reprojection.
(The demo is free of charge, although you cannot save your work
it's very useful to check things out. And I'm
using it on linux through wine). In Global Mapper,
both shapes exported from anyellaWGS84 and
anyellaED50towgs84 are coincident.
This is not the case in QGIS, which seems to ignore the towgs84 field.
(Although I have to try writing that field on my own in QGIS).

Just to note that this is not a paranoid obsession.
People using non-WGS84 datums and saving to
kml to visualize (or let others visualize) spatial
features in google maps will get significant errors,
as google maps is on WGS84. In this particular case,
field plots on grassland were in the forest according
to the previous kml & google maps. Perhaps
an explicit mention of the kml problem would be
worth in the spTransform help page.

Thanks a lot for your help.

Agus



Roger Bivand wrote:
> On Tue, 14 Jul 2009, Agustin Lobo wrote:
> 
>> Inverse transform now:
>>> anyellaWGS84 <- readOGR(dsn=".", layer="parcellesAnyellaWGS84")
>> OGR data source with driver: ESRI Shapefile
>> Source: ".", layer: "parcellesAnyellaWGS84"
>> with  6  rows and  7  columns
>> Feature type: wkbPoint with 2 dimensions
>>> proj4string(anyellaWGS84)
>> [1] " +proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +no_defs 
>> +towgs84=0,0,0"
>>> anyellaED50notowgs84 <- spTransform(anyellaWGS84,CRS("+proj=utm +zone=31 
>> +ellps=intl +units=m +no_defs"))
>>> anyellaED50towgs84 <- spTransform(anyellaWGS84,CRS("+proj=utm +zone=31 
>> +ellps=intl +units=m +no_defs +towgs84=-87,-98,-121"))
>>> coordinates(anyellaWGS84)[1,]
>> coords.x1 coords.x2
>>  421005   4683849
>>> coordinates(anyellaED50notowgs84)[1,]
>> coords.x1 coords.x2
>> 421001.4 4683932.4
>>> coordinates(anyellaED50towgs84)[1,]
>> coords.x1 coords.x2
>> 421097.5 4684050.9
>>
>> which is correct? anyellaED50notowgs84 or anyellaED50towgs84 ?
>>
>> I guess anyellaED50towgs84, but would like to confirm,
> 
> Good guess! In the no +towgs84 case, the ellipsoid is being used but not 
> datum transformation, which leaves a contradiction between a declared 
> ellipsoid and the one inherent in the default datum (WGS84). Providing 
> the +towgs84 transforms to the datum expressed in the given parameters, 
> which is one of the ED50 family.
> 
> The difference is anout 100m in each direction. Maybe try this on an 
> easily identifiable known ground control point where you have positional 
> data from several maps, like the corner of a building or a bridge pier.
> 
> Roger
> 
>>
>> Thanks
>>
>> Agus
>>
> 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: alobolistas.vcf
Type: text/x-vcard
Size: 251 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20090714/1128e598/attachment.vcf>