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

[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



> -----Original Message-----
> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
> Of Agustin Lobo
> Sent: Tuesday, July 14, 2009 9:43 AM
> To: sig-geo
> Subject: [R-sig-Geo] more doubts with spTransform()
> 
> 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