Skip to content

more doubts with spTransform()

4 messages · Agustin Lobo, Roger Bivand, Tomislav Hengl

Original
Agustin Lobo
# 
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,

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/e5f6f165/attachment.vcf>
Roger Bivand
#
On Tue, 14 Jul 2009, Agustin Lobo wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Tomislav Hengl
    

    

      
      #
    

  


  

      
Hi Agustin,

You should really work with some referent point for which you know both local and geographic (WGS84)
coordinates.

Here is an illustration:

http://spatial-analyst.net/wiki/index.php?title=MGI_/_Balkans_coordinate_systems#Validation_of_CRS_p
arameters 

Note that I use Google Earth imagery to validate that everything is fine.

I agree - setting up the right CRS can be a headache, but there are many 'hard' (historic,
technology connected etc) reasons for this. The good thing is that: once you got it right, you soon
forget about all these problems.


HTH,

T. Hengl

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      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:

            

      
      
      
    

        
-------------- 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>