[R-sig-Geo] convUL to convert UTM in Lat/Long

[SchnuteJ at pac.dfo-mpo.gc.ca]

On Mon, 22 Jan 2007, Berta wrote:

> Hi R-geo users, I am trying to convert a data set in UTM (from Spain,
> zone=30) to the corresponding data set in longitude-latitude. Using
> convUL, it gives me strange results. Here I send my code, just in case
> somebody can tell me what i am missing. I suspect it has something to do
> with + proj or something that it is needed.  Thank you so much in
> advance, Berta
> 
> library(PBSmapping)
> datainUTM<-data.frame(cbind(c(559994, 559673),c(4781905,4781583), c(0, 0),
c(1,2)))
> names(datainUTM)<-c("X", "Y", "PID", "POS")
> attr(datainUTM, "projection") <-"UTM"
> attr(datainUTM, "zone") <-30
> datosLL<-convUL(datainUTM)

I'm following the thread that began with Berta's message above. Roger Bivand
correctly states that PBSmapping uses UTM coordinates measured in km, not m.
Currently, the help file doesn't mention this requirement. In the next
version, we'll fix the documentation (convUL.Rd). We may also add another
argument to "convUL" that allows a choice between "km" and "m". Thanks to
Berta and Roger for alerting us to this problem.

For the record, the following R code illustrates the identical conversion
using both PBSmapping and rgdal, based on three points mentioned in the
earlier correspondence:

# Load packages
require(PBSmapping); require(rgdal)

# Input data for PBS Mapping
datainUTM <- data.frame(cbind(
  c(559.994, 559.673, 559.835),
  c(4781.905, 4781.583, 4778.818),
  c(1,1,1), c(1,2,3)))
names(datainUTM) <- c("X", "Y", "PID", "POS")
attr(datainUTM, "projection") <- "UTM"
attr(datainUTM, "zone") <- 30

# PBS Mapping conversion
PBS.LLout <- convUL(datainUTM)

# SP conversion
SP <- SpatialPoints(cbind(1000*datainUTM$X, 1000*datainUTM$Y), 
  proj4string = CRS("+proj=utm +zone=30"))
SP.LLout <- spTransform(SP, CRS("+proj=longlat"))

# Display results
cat("UTM Input:\n"); print(datainUTM)
cat("\nPBS LL Output:\n"); print(PBS.LLout)
cat("\nSP LL Output:\n"); print(SP.LLout)

Happily (given completely independent code), both packages produce the same
results:

UTM Input:
        X        Y PID POS
1 559.994 4781.905   1   1
2 559.673 4781.583   1   2
3 559.835 4778.818   1   3

PBS LL Output:
          X        Y PID POS
1 -2.261705 43.18753   1   1
2 -2.265689 43.18466   1   2
3 -2.263995 43.15975   1   3

SP LL Output:
SpatialPoints:
     coords.x1 coords.x2
[1,] -2.261705  43.18753
[2,] -2.265689  43.18466
[3,] -2.263995  43.15975

We began developing PBS Mapping before we realized that Roger was working on
general spatial classes. In future releases, we hope to take much greater
advantage of his excellent work. PBS Mapping comes with an extensive User's
Guide (PBSmapping-UG.pdf) in the R library tree ..\library\PBSmapping. This
explains the goals of the package and includes figures that illustrate the
possibilities.

Jon

******************************************
Jon Schnute
Pacific Biological Station
3190 Hammond Bay Road
Nanaimo, BC V9T 6N7
CANADA

Tel: 250-756-7146
email: schnutej at pac.dfo-mpo.gc.ca