Problem converting shp ED50 to WGS84 with spTransform and writeOGR
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