On Thu, 2 Apr 2009, Marc Marí Dell'Olmo wrote:
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?
In datum transformation, it is hard to be precise without an accurate set of transformation parameters. These may also vary across an area of interest:
y.wgs <- c(41.4036923115558, 41.38010926480547, 41.41148169051035,
41.409534433320914)
x.wgs <- c(2.145810127258301, 2.1807003021240234, 2.2183799743652344,
2.1170568466186523)
x.utm <- c(428727134, 431587743, 434639340, 426268557)
y.utm <- c(4584118048, 4581465377, 4584855001, 4584607449)
orig <- cbind(x.utm, y.utm)/1000
library(rgdal)
out <- spTransform(SP,
CRS("+proj=utm +ellps=intl +zone=31 +towgs84=-84,-107,-120,0,0,0,0,0"))
for (i in 1:nrow(orig)) print(spDistsN1(coordinates(out)[i,,drop=FALSE],
orig[i, ]))
is:
[1] 164.4001
[1] 156.7712
[1] 167.2330
[1] 44.72934
crds <- cbind(x.wgs, y.wgs)
for (i in 2:nrow(crds)) print(spDistsN1(crds[1:(i-1),,drop=FALSE],
crds[i,], longlat=TRUE)*1000)
suggests that the points are not too far apart, under 10km, so I would not
expect differences of this scale. How were the UTM projected and WGS
geographical coordinates measured - I guess that the ED50 UTM coordinates
are from a legacy/archive source (manually digitized map?) and the WGS
geographical coordinates from field measured GPS?
Roger
spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[1,]),longlat=TRUE)[1]*1000[1] 25.66325spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[2,]),longlat=TRUE)[2]*1000[1] 15.8736spDistsN1(coordinates(SP.wgs),coordinates(SP.utm.wgs[3,]),longlat=TRUE)[3]*1000[1] 158.1777spDistsN1(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, RogerThank 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.117057coordinates(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]
-- 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
