[R-sig-Geo] convert point and map coordinates UTM to WSG

Thu Apr 2 16:43:30 CEST 2009

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.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