[R-sig-Geo] problem with lcca inverse projection

[Daniel Kelley]

Hello. I am having difficulties with the inverse version of the lcca projection, so I made a test code that shows the problem. I’m enclosing a script after my address information, along with the output. Basically, as you can see by the “BAD” flags at the ends of the lines, my tests with lcca are failing, whereas those with lcc are passing. I’m trying two spots on the globe to project and then unproject, and two lat_0 values.

Q: Am I doing something wrong, or is there a problem with the lcca projection? 

Thanks in advance, for any advice.

Dan.

Dan Kelley, Oceanography Department, Dalhousie University






R version 3.4.1 (2017-06-30) -- "Single Candle"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(rgdal)
> LON <- -135
> for (p in c("lcca", "lcc ")) {
+     for (lon0 in c(-180, -170)) {
+         for (lat0 in c(-90, -80)) {
+             for (LAT in c(-60, 60)) {
+                 proj <- sprintf("+proj=%s +lon_0=%.0f +lat_0=%.0f",
+                                 p, lon0, lat0)
+                 xy <- project(cbind(LON, LAT), proj=proj)
+                 ll <- project(xy, proj=proj, inv=TRUE)
+                 dlon <- abs(ll[1,1] - LON)
+                 dlat <- abs(ll[1,2] - LAT)
+                 ok <- if (dlon < 0.0001 && dlat < 0.0001) " OK" else "BAD"
+                 cat(sprintf("\"%s\" : %4.0fE %3.0fN -> %9.0fm %9.0fm -> %8.3fE %8.3fN %s\n",
+                             proj, LON, LAT, xy[1,1], xy[1,2], ll[1,1], ll[1,2], ok))
+             }
+         }
+         cat("\n")
+     }
+ }
"+proj=lcca +lon_0=-180 +lat_0=-90" : -135E -60N ->   2475298m   2475298m ->   45.000E -120.000N BAD
"+proj=lcca +lon_0=-180 +lat_0=-90" : -135E  60N ->  25074305m  25074305m ->   45.000E -240.000N BAD
"+proj=lcca +lon_0=-180 +lat_0=-80" : -135E -60N ->   2378510m   1307652m ->   42.223E -117.855N BAD
"+proj=lcca +lon_0=-180 +lat_0=-80" : -135E  60N ->  22316311m  21727017m ->   42.223E -224.927N BAD

"+proj=lcca +lon_0=-170 +lat_0=-90" : -135E -60N ->   2007862m   2867523m ->   45.000E -120.000N BAD
"+proj=lcca +lon_0=-170 +lat_0=-90" : -135E  60N ->  20339263m  29047477m ->   45.000E -240.000N BAD
"+proj=lcca +lon_0=-170 +lat_0=-80" : -135E -60N ->   1926825m   1678569m ->   42.223E -117.855N BAD
"+proj=lcca +lon_0=-170 +lat_0=-80" : -135E  60N ->  18078384m  25207137m ->   42.223E -224.927N BAD

"+proj=lcc  +lon_0=-180 +lat_0=-90" : -135E -60N ->  13525767m -25044282m -> -135.000E  -60.000N  OK
"+proj=lcc  +lon_0=-180 +lat_0=-90" : -135E  60N ->   2588932m  -4793661m -> -135.000E   60.000N  OK
"+proj=lcc  +lon_0=-180 +lat_0=-80" : -135E -60N ->  13525767m  32572859m -> -135.000E  -60.000N  OK
"+proj=lcc  +lon_0=-180 +lat_0=-80" : -135E  60N ->   2588932m  52823480m -> -135.000E   60.000N  OK

"+proj=lcc  +lon_0=-170 +lat_0=-90" : -135E -60N ->  10693582m -26378205m -> -135.000E  -60.000N  OK
"+proj=lcc  +lon_0=-170 +lat_0=-90" : -135E  60N ->   2046831m  -5048984m -> -135.000E   60.000N  OK
"+proj=lcc  +lon_0=-170 +lat_0=-80" : -135E -60N ->  10693582m  31238936m -> -135.000E  -60.000N  OK
"+proj=lcc  +lon_0=-170 +lat_0=-80" : -135E  60N ->   2046831m  52568157m -> -135.000E   60.000N  OK