[R-sig-Geo] Problem converting shp ED50 to WGS84 with spTransform and writeOGR

[Roger Bivand]

On Tue, 3 Dec 2013, Marc Marí Dell'Olmo wrote:

> Dear all,
>
> I'm trying to convert a map of Barcelona (in shp format) with reference
> ED50 to a kml with Google Maps reference. I have used the following syntax,
> but it creates a slightly displaced map (maybe just a few meters) and do
> not understand why!

When you understand the +datum= and +towgs84= tags, you'll understand. You 
haven't given a datum, and ED50 is not a single standard. You also did not 
quote yout sessionInfo() output, nor the startup messages from loading 
rgdal.

My guess is that you are using Windows, are using an older Windows rgdal 
binary from CRAN with PROJ.4 4.7.*, and have not read:

https://stat.ethz.ch/pipermail/r-sig-geo/2013-November/019884.html

So, guessing further, if you upgrade to the latest Windows rgdal binary 
from CRAN, which has PROJ.4 4.8.*, you will see:

> CRS("+init=epsg:23031")
CRS arguments:
  +init=epsg:23031 +proj=utm +zone=31 +ellps=intl
+towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs

which gives a three-parameter transformation from the local Spanish 
version of ED50 to WGS84 accepted by EPSG as reasonable.

If you do not give +towgs84= or a known +datum=, or they (one, other, or 
both) are not returned from +init=epsg:*, spTransform assumes that the 
datum is WGS84. The best resource for details on datums (and local 
transformation parameters) is:

http://www.asprs.org/Grids-Datums.html

with Spain July 2000, quoting the three parameters as:

\delta X = –84m ± 5m, \delta Y = –107m ± 6m, \delta Z = –120m ± 3m.

which are not quite the same as the EPSG values, so you could use those or 
values near those for possibly increased accuracy. In any case, using 
+towgs84= with the EPSG values will get you much closer than not using 
+towgs84=.

Hope this clarifies,

Roger


>
>
> setwd("C:/dir")
> barcelona <- readOGR("C:/dir", "Dist_postals_BCN")
> proj4string(barcelona) <- CRS("+proj=utm +zone=31 +ellps=intl +units=m
> +no_defs")
> # Or proj4string(barcelona) <- CRS("+init=epsg:23031")
> x1 <- spTransform(barcelona, CRS("+proj=longlat +datum=WGS84 +no_defs"))
> writeOGR(x1, "districtes.kml", "x1" ,  driver="KML")
>
>
> If you want to try, you can download the ED50 shp map here:
> https://dl.dropboxusercontent.com/u/14934021/map.zip
>
> And here de displaced kml:
> https://dl.dropboxusercontent.com/u/14934021/districtes.kml
>
> Thank you!!
>
> Marc
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
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