1\name{project}
2\alias{project}
3
4\title{Projection of coordinate matrices}
5\description{
6  Interface to the PROJ.4 library of projection functions for geographical position data, no datum transformation possible. Use \code{spTransform()} for extended support.
7}
8\usage{
9project(xy, proj, inv = FALSE, use_ob_tran=FALSE, legacy=TRUE,
10 allowNAs_if_not_legacy=FALSE, coordOp = NULL, verbose = FALSE,
11 use_aoi=TRUE)
12}
13
14\arguments{
15  \item{xy}{ 2-column matrix of coordinates }
16  \item{proj}{ character string of projection arguments; the arguments must be entered exactly as in the PROJ.4 documentation, in particular there cannot be any white space in +<arg>=<value> strings, and successive such strings can only be separated by blanks. }
17  \item{inv}{ default FALSE, if TRUE inverse projection to geographical coordinates }
18  \item{use_ob_tran}{default FALSE, if TRUE and \dQuote{+proj=ob_tran}, use General Oblique Transformation with internalised from/to projection reversal; the user oblique transforms forward rather than inverse.}
19  \item{legacy}{default TRUE, if FALSE, use transform C functions (enforced internally for Windows 32-bit platforms)}
20  \item{allowNAs_if_not_legacy}{used if legacy is FALSE, default FALSE; introduced to handle use of NAs as object separators in \pkg{oce}}
21  \item{coordOp}{default NULL, for PROJ >= 6 used to pass through a pre-defined coordinate operation}
22  \item{verbose}{default FALSE, for PROJ >=6 used to show the coordinate operation used}
23  \item{use_aoi}{With PROJ >= 6, use the area of interest defined as the range of \code{xy} in limiting the search for candidate coordinate operations; set FALSE if \code{use_ob_tran} is TRUE}
24}
25\details{
26  Full details of projection arguments available from website below, and examples in file "epsg" in the data directory installed with PROJ.4.
27
28Note that from PROJ.4 4.9.3, the definition of UTM is changed from TMERC to ETMERC; see example.
29}
30\value{
31  A two column matrix with projected coordinates.
32}
33\references{\url{https://proj.org/}}
34\author{Barry Rowlingson, Roger Bivand \email{Roger.Bivand@nhh.no}}
35
36\note{ The locations of Hawaii and Alaska in the data source are (putting it mildly) arbitrary, please avoid airlines using these positions.}
37
38\seealso{ \code{\link[sp]{CRS-class}}, \code{\link{spTransform-methods}} }
39
40\examples{
41data(state)
42res <- project(cbind(state.center$x, state.center$y),
43 "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84")
44res1 <- project(res, "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84",
45 inv=TRUE)
46summary(res1 - cbind(state.center$x, state.center$y))
47plot(cbind(state.center$x, state.center$y), asp=1, type="n")
48text(cbind(state.center$x, state.center$y), state.abb)
49plot(res, asp=1, type="n")
50text(res, state.abb)
51broke_proj <- FALSE
52pv <- .Call("PROJ4VersionInfo", PACKAGE="rgdal")[[2]]
53# https://github.com/OSGeo/PROJ/issues/1525
54if (pv >= 600 && pv < 620) broke_proj <- TRUE
55if (!broke_proj) {
56crds <- matrix(data=c(9.05, 48.52), ncol=2)
57a <- project(crds, paste("+proj=ob_tran +o_proj=longlat",
58 "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"),
59 use_ob_tran=TRUE)
60a
61#should be (-5.917698, -1.87195)
62project(a, paste("+proj=ob_tran +o_proj=longlat",
63 "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"),
64 inv=TRUE, use_ob_tran=TRUE)
65#added after posting by Martin Ivanov
66}
67#
68getPROJ4VersionInfo()
69# Test for UTM == TMERC (<= 4.9.2) or UTM == ETMERC (> 4.9.2)
70nhh <- matrix(c(5.304234, 60.422311), ncol=2)
71nhh_utm_32N_P4 <- project(nhh, "+init=epsg:3044")
72nhh_tmerc_P4 <- project(nhh, paste("+proj=tmerc +k=0.9996 +lon_0=9",
73 "+x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
74nhh_etmerc_P4 <- project(nhh, paste("+proj=etmerc +k=0.9996 +lon_0=9",
75 "+x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
76all.equal(nhh_utm_32N_P4, nhh_tmerc_P4, tolerance=1e-9, scale=1)
77# UTM == TMERC: PROJ4 <=4.9.2
78all.equal(nhh_utm_32N_P4, nhh_etmerc_P4, tolerance=1e-9, scale=1)
79# UTM == ETMERC: PROJ4 > 4.9.2
80unis <- matrix(c(15.653453, 78.222504), ncol=2)
81unis_utm_33N_P4 <- project(unis, "+init=epsg:3045")
82unis_tmerc_P4 <- project(unis, paste("+proj=tmerc +k=0.9996 +lon_0=15",
83 "+x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
84unis_etmerc_P4 <- project(unis, paste("+proj=etmerc +k=0.9996 +lon_0=15",
85 "+x_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
86all.equal(unis_utm_33N_P4, unis_tmerc_P4, tolerance=1e-9, scale=1)
87# UTM == TMERC: PROJ4 <=4.9.2
88all.equal(unis_utm_33N_P4, unis_etmerc_P4, tolerance=1e-9, scale=1)
89# UTM == ETMERC: PROJ4 > 4.9.2
90#pv <- attr(getPROJ4VersionInfo(), "short")
91#if (pv < 500) {
92# valgrind leakages in some cases for PROJ >= 5; many non-projection proj values added
93# available projections and their inverses if provided
94# For >=4.9.3 returns non-finite points rather than needing crash protection
95projs <- as.character(projInfo()$name)
96res <- logical(length(projs))
97names(res) <- projs
98msgs <- character(length(projs))
99names(msgs) <- projs
100owarn <- options("warn")$warn
101options(warn=2L)
102for (i in seq(along=res)) {
103  iprs <- paste("+proj=", projs[i], sep="")
104  xy <- try(project(cbind(0, 0), iprs, legacy=TRUE, use_aoi=FALSE), silent=TRUE)
105  if (inherits(xy, "try-error")) {
106    res[i] <- NA
107    msgs[i] <- paste("fwd:", strsplit(xy, "\n")[[1]][2])
108  } else if(any(abs(xy) > 1e+08)) {
109    res[i] <- NA
110    msgs[i] <- paste("fwd: huge value")
111  } else {
112    out <- try(project(xy, iprs, inv=TRUE, legacy=TRUE, use_aoi=FALSE), silent=TRUE)
113    if (inherits(out, "try-error")) {
114      res[i] <- NA
115      msgs[i] <- paste("inv:", strsplit(out, "\n")[[1]][2])
116    } else {
117      res[i] <- isTRUE(all.equal(cbind(0,0), out))
118    }
119  }
120}
121options(warn=owarn)
122df <- data.frame(res=unname(res), msgs=unname(msgs), row.names=names(res))
123# projection and inverse projection failures
124# fwd: missing parameters
125# inv: mostly inverse not defined
126df[is.na(df$res),]
127# inverse not equal to input
128# (see http://lists.maptools.org/pipermail/proj/2011-November/006015.html)
129df[!is.na(df$res) & !df$res,]
130# inverse equal to input
131row.names(df[!is.na(df$res) & df$res,])
132#}
133# oce data representation with NAs
134ll <- structure(c(12.1823368669203, 11.9149630062421, 12.3186076188739,
13512.6207597184845, 12.9955172054652, 12.6316117692658, 12.4680041846297,
13612.4366882666609, NA, NA, -5.78993051516384, -5.03798674888479,
137-4.60623015708619, -4.43802336997614, -4.78110320396188, -4.99127125409291,
138-5.24836150474498, -5.68430388755925, NA, NA), .Dim = c(10L,
1392L), .Dimnames = list(NULL, c("longitude", "latitude")))
140try(xy0 <- project(ll, "+proj=moll", legacy=TRUE))
141if (!PROJis6ormore()) { # legacy=TRUE PROJ >= 6
142try(xy1 <- project(ll, "+proj=moll", legacy=FALSE, allowNAs_if_not_legacy=FALSE))
143try(xy2 <- project(ll, "+proj=moll", legacy=FALSE, allowNAs_if_not_legacy=TRUE))
144if (exists("xy0")) all.equal(xy0, xy2)
145}
146if (!exists("xy0")) xy0 <- structure(c(1217100.8468177, 1191302.229156,
1471232143.28841193, 1262546.27733232, 1299648.82357849, 1263011.18154638,
1481246343.17808186, 1242654.33986052, NA, NA, -715428.207551599,
149-622613.577983058, -569301.605757784, -548528.530156422, -590895.949857199,
150-616845.926397351, -648585.161643274, -702393.1160979, NA, NA),
151.Dim = c(10L, 2L), .Dimnames = list(NULL, c("longitude", "latitude")))
152try(ll0 <- project(xy0, "+proj=moll", inv=TRUE, legacy=TRUE))
153if (!PROJis6ormore()) { # legacy=TRUE PROJ >= 6
154try(ll1 <- project(xy0, "+proj=moll", inv=TRUE, legacy=FALSE, allowNAs_if_not_legacy=FALSE))
155try(ll2 <- project(xy0, "+proj=moll", inv=TRUE, legacy=FALSE, allowNAs_if_not_legacy=TRUE))
156if (exists("ll0")) all.equal(ll0, ll2)
157}
158if (exists("ll0")) all.equal(ll0, ll)
159}
160\keyword{spatial}
161