convert point and map coordinates UTM to WSG

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.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 <Roger.Bivand at nhh.no>
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 <p.hiemstra at geo.uu.nl>

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
R-sig-Geo at stat.math.ethz.ch
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: Roger.Bivand at nhh.no

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: Roger.Bivand at nhh.no

convert point and map coordinates UTM to WSG

Thread (4 messages)