Dear Roger, thank you for your help. There is something I don't quite understand. If I compare distances between the converted coordinates (SP.utm.wgs) with the coordinates (SP.wgs) that I use to compare with (in pairs), the distances are more variant than I expected. These go from approximately 15 meters to 192. I thought that the error would be more or less constant. Is it normal?
> spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[1,]),longlat=TRUE)[1]*1000 [1] 25.66325 > spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[2,]),longlat=TRUE)[2]*1000 [1] 15.8736 > spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[3,]),longlat=TRUE)[3]*1000 [1] 158.1777 > spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[4,]),longlat=TRUE)[4]*1000 [1] 194.8797 Thank you, Marc 2009/4/1 Roger Bivand <[email protected]> > > On Wed, 1 Apr 2009, Marc Marí Dell'Olmo wrote: > >> Dear all, >> >> I have tried to transform coordinates from UTM Zone-31 ED-50 Spain Portugal >> (EUR-D code in this web >> http://earth-info.nga.mil/GandG/coordsys/onlinedatum/CountryEuropeTable.html# >> EURD) to coordinates WSG84 ( google maps coordinates) with the >> "ptransform" of the proj4 library and I can not do it right (first syntax). >> I have tried also using the "spTransform" of the rgdal library (second >> syntax) and I have not been successful. >> >> I attach the syntax used, I got 4 points with coordinates ED ESP-50 ( "utm" >> in the syntax) and their respective coordinates WSG84 (wgs.gold "in the >> syntax). Can anyone tell me where is the error? > > So far several. One is reversing x and y values, another saying "latlong" > when you mean "longlat", a third is the units of your input UTM (I think mm, > not m), here, I'm going from longlat to UTM: > > library(rgdal) > > y.wgs <- c(41.4036923115558, 41.38010926480547, 41.41148169051035, > 41.409534433320914) > x.wgs <- c(2.145810127258301, 2.1807003021240234, 2.2183799743652344, > 2.1170568466186523) > SP <- SpatialPoints(cbind(x.wgs, y.wgs), > proj4string=CRS("+proj=longlat +datum=WGS84")) > spTransform(SP, > CRS("+proj=utm +ellps=intl +zone=31 +towgs84=-84,-107,-120,0,0,0,0,0")) > > The remaining differences are from the inadequacy of the +towgs84 string; > Grids & Datums (http://www.asprs.org/resources/grids/) - like your source - > gives +/- ranges for the transformation parameters, so you'll need to fish > unless you can get a better authority (and there are few better than Grids & > Datums). > > Hope this helps, > > Roger > > > >> >> Thank you very much, >> >> Marc >> >>> #First sintax >>> library(proj4) >>> x.utm <- c(428727134,431587743,434639340,426268557) >>> y.utm <- c(4584118048,4581465377,4584855001,4584607449) >>> x.wgs <- >> >> c(41.4036923115558,41.38010926480547,41.41148169051035,41.409534433320914) >>> >>> y.wgs <- >> >> c(2.145810127258301,2.1807003021240234,2.2183799743652344,2.1170568466186523) >> >>> >>> utm <- cbind(x.utm, y.utm) >>> wgs.new <- ptransform(utm, "+proj=utm +ellps=intl +zone=31 >> >> +towgs84=-84,-107,-120,0,0,0,0,0 +no_defs", >> + "+proj=latlong +datum=WGS84") >>> >>> wgs.new >> >> [,1] [,2] [,3] >> [1,] -2.236353 1.570775 39.63333 >> [2,] -2.236353 1.570775 39.63333 >> [3,] -2.236353 1.570775 39.63333 >> [4,] -2.236353 1.570775 39.63333 >>> >>> (wgs.gold <- cbind(x.wgs, y.wgs)) >> >> x.wgs y.wgs >> [1,] 41.40369 2.145810 >> [2,] 41.38011 2.180700 >> [3,] 41.41148 2.218380 >> [4,] 41.40953 2.117057 >>> >>> >>> >>> #Second sintax >>> library(rgdal) >>> (data.utm <- data.frame(x.utm,y.utm,x.wgs,y.wgs)) >> >> x.utm y.utm x.wgs y.wgs >> 1 428727134 4584118048 41.40369 2.145810 >> 2 431587743 4581465377 41.38011 2.180700 >> 3 434639340 4584855001 41.41148 2.218380 >> 4 426268557 4584607449 41.40953 2.117057 >>> >>> coordinates(data.utm) <- c("x.utm", "y.utm") >>> proj4string(data.utm) <- CRS("+proj=utm +ellps=intl +zone=31 >> >> +towgs84=-84,-107,-120,0,0,0,0,0 +no_defs") >>> >>> (wgs.new2 <- spTransform(data.utm, CRS("+proj=latlong +datum=WGS84"))) >> >> coordinates x.wgs y.wgs >> 1 (-128.134, 89.9988) 41.40369 2.145810 >> 2 (-128.134, 89.9988) 41.38011 2.180700 >> 3 (-128.134, 89.9988) 41.41148 2.218380 >> 4 (-128.134, 89.9988) 41.40953 2.117057 >> >> >> >> 2009/3/23 Paul Hiemstra <[email protected]> >> >>> Hi, >>> >>> To transform coordinates you can use the spTransform function from the >>> rgdal pacakge. It uses proj4 to do the transformations. Proj uses a string >>> to describe a certain projection, the so called proj4string. If you can >>> determine the proj4string for both your projections you can transform >>> between them. The first string (WGS84) looks like this: >>> >>> "+proj=longlat +datum=WGS84" >>> >>> The second (UTM) something like: >>> >>> "+proj=utm +zone=yourzone +datum=yourdatum" >>> >>> See the documentation of spTransform for more information on how to use the >>> function. >>> >>> cheers and hth, >>> Paul >>> >>> >>> Marc Mar? Dell'Olmo wrote: >>> >>>> Dear R-users, >>>> >>>> Can anyone tell me how to convert point coordinates from WSG84 (google >>>> maps coordinates) format to UTM HUSO-31 ED-50. And if I have a map >>>> (shp format) with UTM HUSO-31 ED-50 coordinates, how can I convert to >>>> WSG84 (google maps coordinates) coordinates. Should I install some >>>> external software? >>>> >>>> Can anyone tell me the instructions? >>>> >>>> Thank you, >>>> >>>> Marc >>>> >>>> _______________________________________________ >>>> R-sig-Geo mailing list >>>> [email protected] >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>>> >>> >>> >>> -- >>> Drs. Paul Hiemstra >>> Department of Physical Geography >>> Faculty of Geosciences >>> University of Utrecht >>> Heidelberglaan 2 >>> P.O. Box 80.115 >>> 3508 TC Utrecht >>> Phone: +3130 274 3113 Mon-Tue >>> Phone: +3130 253 5773 Wed-Fri >>> http://intamap.geo.uu.nl/~paul <http://intamap.geo.uu.nl/%7Epaul> >>> >>> >> >> [[alternative HTML version deleted]] >> >> > > -- > Roger Bivand > Economic Geography Section, Department of Economics, Norwegian School of > Economics and Business Administration, Helleveien 30, N-5045 Bergen, > Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 > e-mail: [email protected] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
