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

[Roger Bivand]

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