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@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
