Thank you very much for your answer. Now it works! Marc
2013/12/3 Roger Bivand <[email protected]> > 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 >> [email protected] >> 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: [email protected] > [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
