hypertidy / reproj

Reproject with the PROJ library
http://hypertidy.github.io/reproj/
3 stars 0 forks source link

document the project 2-row thing #5

Closed mdsumner closed 5 years ago

mdsumner commented 5 years ago

proj4 has a problem with a 2-row matrix, reported on 2017-01-12 and with a patch in the handling applied below.

there's an issue in the argument handling for project(). If x is a 2-row matrix then it's effectively transposed, and throws an error if the flipped longitude exceeds the prescribed latitude bounds, but not otherwise. I've attached a demonstration script, these lines cause the problem: if (d[1] == 2) { x <- xy[1, ] y <- xy[2, ] } I'd quite like to see this fixed, so I'm happy to help if need be.

(Funnily much the same as raster::extent() on a 2-row matrix, also unfixed).

#' simulate the xy argument handling of proj4::project
project_arg_handler <- 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")
      if (d[1] == 2) {
        x <- xy[1, ]
        y <- xy[2, ]
      }
      else {
        x <- xy[, 1]
        y <- xy[, 2]
      }
    }
  }
  list(x, y)
}

## 2-column 3-row matrix long-lat
project_arg_handler(cbind(c(24, -91, 10), c(0, 1, 10)))
# [[1]]
# [1]  24 -91  10
#
# [[2]]
# [1]  0  1 10

## only 2-rows, mixes x/y
project_arg_handler(cbind(c(24, -91), c(0, 1)))
# [[1]]
# [1] 24  0
#
# [[2]]
# [1] -91   1

## causes problems if longitudes out of latitude range
## this transforms the values, but is the wrong answer
proj4::project(cbind(c(24, -89), c(0, 1)), "+proj=laea")
# [,1]       [,2]
# [1,] 63453.59 -8937611.2
# [2,]     0.00   111193.5

proj4::project(cbind(c(24, -91), c(0, 1)), "+proj=laea")
##Error in proj4::project(cbind(c(24, -91), c(0, 1)), "+proj=laea") :
#  latitude or longitude exceeded limits

## the right answer without the out of bounds
proj4::project(cbind(c(24, -89, 0), c(0, 1, 0)), "+proj=laea")
# [1,]  2649210      0.0
# [2,] -8929633 155891.1
# [3,]        0      0.0
mdsumner commented 5 years ago

This was fixed: https://github.com/s-u/proj4/issues/1

mdsumner commented 5 years ago

It's not relevant to reproj which uses ptransform. But yeah, watch out.

mdsumner commented 4 years ago

And released in March 2020: https://cran.r-project.org/web/packages/proj4/NEWS