[R-sig-Geo] Proj4 - latitude or longitude exceeded limits
Roger Bivand
Roger.Bivand at nhh.no
Wed Jan 18 10:59:11 CET 2017
On Wed, 18 Jan 2017, Roger Bivand wrote:
> On Tue, 17 Jan 2017, Michael Sumner wrote:
>
>> Watch out with proj4 btw, it has an obscure bug where a 2 row matrix
>> input
>> is transposed, mixing longitude and latitude in a way that triggers an
>> out
>> of bounds error for latitude, only when the mixed abs(longitude) is
>> greater
>> than 90.
>
> You mean:
>
>> library(rgdal)
>> crds <- rbind(c(90, 1), c(91, 2))
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> [,1] [,2]
> [1,] 833927.9 110682.8
> [2,] 945193.9 221604.0
>> crds <- rbind(c(90, 1), c(91, 2), c(92, 3))
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> [,1] [,2]
> [1,] 833927.9 110682.8
> [2,] 945193.9 221604.0
> [3,] 1056324.8 332866.1
>> crds <- matrix(c(91, 2), nrow=1)
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> [,1] [,2]
> [1,] 945193.9 221604
>
> and
>
>> library(proj4)
>> crds <- rbind(c(90, 1), c(91, 2))
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> Error in project(crds, "+proj=utm +zone=45 +ellps=WGS84") :
> latitude or longitude exceeded limits
>> crds <- rbind(c(90, 1), c(90, 2))
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> Error in project(crds, "+proj=utm +zone=45 +ellps=WGS84") :
>> crds <- rbind(c(90, 1), c(91, 2), c(92, 3))
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> [,1] [,2]
> [1,] 833927.9 110682.8
> [2,] 945193.9 221604.0
> [3,] 1056324.8 332866.1
>> crds <- matrix(c(91, 2), nrow=1)
>> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
> [,1] [,2]
> [1,] 945193.9 221604
>
> right?
>
> I'll take a look and if I see the bug, send a patch to Simon. Please confirm
> that this is what you see.
Lines 20-26 in proj4/R/project.R flip an input matrix with two rows, so
putting one of the >=abs(90) longitude values in y - for unknown reasons.
# Keep matrix orientation fixed (RSB/MS 170118)
# if (d[1]==2) {
# x <- xy[1,]
# y <- xy[2,]
# } else {
x <- xy[,1]
y <- xy[,2]
# }
With the offending test commented out, we're OK.
> library(proj4)
> crds <- rbind(c(90, 1), c(91, 2))
> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
[,1] [,2]
[1,] 833927.9 110682.8
[2,] 945193.9 221604.0
The workaround if proj4 is not patched would be to do:
> library(proj4)
> crds <- rbind(c(90, 1), c(91, 2))
> if (dim(crds)[1] == 2) crds <- t(crds)
> project(crds, "+proj=utm +zone=45 +ellps=WGS84")
[,1] [,2]
[1,] 833927.9 110682.8
[2,] 945193.9 221604.0
proj4::ptransform() is OK (provided an input object in degrees is
converted to radians by the user).
Could you check (edited proj4/R/project.R attached)?
Roger
>
> Best wishes,
>
> Roger
>
>>
>> Use the sf package, or rgdal, or mapproj.
>>
>> Cheers, Mike
>>
>> On Wed, Jan 18, 2017, 05:51 Roger Bivand <Roger.Bivand at nhh.no> wrote:
>>
>> > On Tue, 17 Jan 2017, aurelie.tschopp at vetsuisse.unibe.ch wrote:
>> >
>> > > Dear all,
>> > >
>> > > I am analysing GPS Data from Chad. I need to first project the
>> > > latitude
>> > > and longitudes. I use the package Proj4, the function project for
>> > > this.
>> > > I have 2 devices and with the data of the 1st device, the projection
>> > > worked perfectly. But with the data from the other device, I got this
>> > > error message: "Error in project(loc, "+proj=utm +zone=33
>> > > +ellps=WGS84
>> > > +datum=WGS84 +units=m +no_defs") : latitude or longitude exceeded
>> > > limits"
>> > >
>> > > I found that Chad was UTM 33P (Wikipedia UTM:
>> > >
>> > https: //en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
>> > > ) , but some people could use the notification 33N: 33 for the zone
>> > > ) and N
>> > > for the North hemisphere.
>> > >
>> > > My code is:
>> > > loc <- dat[,c("Longitude","Latitude")
>> > > loc_Zone33 <- project(loc, "+proj=utm +zone=33 +ellps=WGS84
>> > > +datum=WGS84
>> > > +units=m +no_defs")
>> > >
>> > > Does anybody know where my problem could be?
>> >
>> > Well, the error message is coming from inside the proj.4 library: did
>> > you
>> > look at the input object? If one device was returning decimal
>> > geographical
>> > coordinates and the other a string with degrees, minutes and seconds in
>> > one "number":
>> >
>> > library(rgdal)
>> > project(matrix(c(19.0, 15.0), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> > project(matrix(c(1900, 1500), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> >
>> > In this setting I get:
>> >
>> > > project(matrix(c(19.0, 15.0), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> > [,1] [,2]
>> > [1,] 930334.7 1662218
>> > > project(matrix(c(1900, 1500), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> > [,1] [,2]
>> > [1,] Inf Inf
>> > Warning message:
>> > In project(matrix(c(1900, 1500), nrow = 1), "+proj=utm +zone=33
>> > +ellps=WGS84") :
>> > 1 projected point(s) not finite
>> >
>> > but not your error. I suspect that the matrix you are passing contains
>> > non-geographical coordinates. If you are using the proj4 package (odd
>> > choice but there you are):
>> >
>> > > library(proj4)
>> > > project(matrix(c(19.0, 15.0), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> > [,1] [,2]
>> > [1,] 930334.7 1662218
>> > > project(matrix(c(1900, 1500), nrow=1), "+proj=utm +zone=33
>> > +ellps=WGS84")
>> > Error in project(matrix(c(1900, 1500), nrow = 1), "+proj=utm +zone=33
>> > +ellps=WGS84") :
>> > latitude or longitude exceeded limits
>> >
>> > confirming that it doesn't handle exceptions other than by a hard fail.
>> >
>> > Roger
>> >
>> >
>> > > Thank you
>> > > Aurélie
>> > >
>> > > _______________________________________________
>> > > R-sig-Geo mailing list
>> > > R-sig-Geo at r-project.org
>> > > https: //stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >
>> >
>> > --
>> > Roger Bivand
>> > Department of Economics, Norwegian School of Economics,
>> > Helleveien 30, N-5045 Bergen, Norway.
>> > voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
>> > http: //orcid.org/0000-0003-2392-6140
>> > https: //scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>> > http: //depsy.org/person/434412
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at r-project.org
>> > https: //stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>> >
>> > University of Tasmania Electronic Communications Policy (December,
>> > 2014).
>> > This email is confidential, and is for the intended recipient only.
>> > Access, disclosure, copying, distribution, or reliance on any of it by
>> > anyone outside the intended recipient organisation is prohibited and
>> > may be
>> > a criminal offence. Please delete if obtained in error and email
>> > confirmation to the sender. The views expressed in this email are not
>> > necessarily the views of the University of Tasmania, unless clearly
>> > intended otherwise.
>> > .
>> >
>>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412
-------------- next part --------------
# Copyright (c) 2007 by Simon Urbanek
#
# part of proj4 R package, license: GPL v2
project <- function(xy, proj, inverse=FALSE, degrees=TRUE, silent=FALSE, ellps.default="sphere") {
proj <- .proj2char(proj, ellps.default=ellps.default)
if (is.list(xy)) {
if (length(xy)<2) stop("input must be at least 2-dimensional")
if (length(xy)>2 && !silent) warning("more than two dimensions found, using first two")
x <- xy[[1]]
y <- xy[[2]]
} else {
d <- dim(xy)
if (is.null(d) && length(xy)==2) {
x <- xy[1]
y <- xy[2]
} else {
if (length(d) != 2 || (d[1]!=2 && d[2]!=2))
stop("input must be 2-dimensional")
# Keep matrix orientation fixed (RSB/MS 170118)
# if (d[1]==2) {
# x <- xy[1,]
# y <- xy[2,]
# } else {
x <- xy[,1]
y <- xy[,2]
# }
}
}
if (!is.numeric(x)) x <- as.numeric(x)
if (!is.numeric(y)) y <- as.numeric(y)
if (length(x) < length(y)) {
if (!silent) warning("x is shorter than y, recycling")
x <- rep(x, length.out=length(y))
}
if (length(y) < length(x)) {
if (!silent) warning("y is shorter than x, recycling")
y <- rep(y, length.out=length(x))
}
n <- length(x)
f <- 0:0
if (inverse) f <- f + 1:1
if (degrees) f <- f + 2:2
res <- .C("project", as.character(proj),
as.integer(n), x=as.double(x), y=as.double(y),
f, NAOK=TRUE, PACKAGE=.package.name)
if (is.list(xy))
list(x=res$x, y=res$y)
else
cbind(res$x, res$y)
}
More information about the R-sig-Geo
mailing list