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