EOW (development: edge-of-world) #1848

dankelley · 2021-07-15T13:15:48Z

This issue thread is just for me to take notes, and so I won't bother to explain either my intentions or the work I'm doing. If my ideas end up going somewhere, I'll explain them separately.

No file mentioned here is necessarily up-to-date on GH. Any comment may supersede previous comments. And I'll feel free to go back and edit comments as I see fit.

dankelley · 2021-07-15T13:18:56Z

Hm sandbox/dk/eow/eow02.R with a lat_0 gives odd results. These seem independent of the number of steps in the bisection search. Next, I'll try using a grid pattern, as a way to check on whether the bisector has got caught on some sort of weird echo (not the right word, but, after all, this issue is just for my own notes).

dankelley · 2021-07-15T13:31:40Z

Hm, The weird outline is not a result of bisection, as the grid test shows. I think map2lonlat() must be doing something odd.

Note the interesting thing at the bottom-centre of the plot, which the grid reveals whereas the bisection didn't. (The bisection gives us much higher precision of the boundary along a line. This grid is pretty slow because it's doing more computations.)

dankelley · 2021-07-15T14:21:40Z

Focussing in on one spot, and doing direct calls to the sf function (to remove possibility that the hiccough is caused by an oce wrapper), I see that these x,y points are just not invertible in this projection.

Maybe the plot is (somehow) using a different projection?
Maybe the inverse computation is just wrong?

dankelley · 2021-07-15T14:34:22Z

A test case: -75.612 N and -85.352 E (using google maps to find the spot. I like their new interface, but it won't let you centre on either pole ... sheesh)

dankelley · 2021-07-15T16:26:08Z

I made a test that is independent of oce, as shown in the reprex pasted below. Red dots in the diagram indicate that the projected points can be reverse-projected. By contrast, the blue dots can be forward-projected, but the reverse fails to converge.

The points display a lon-lat grid. The South pole is visible near the bottom, flanked by problematic (blue) points. Problematic points can also be seen at around 11 and 13 o'clock on the diagram.

I have asked @richardsc, @clayton33 and @j-harbin to try to reproduce this, so that I can learn whether it's a problem that is particular to my machine. (I am using a beta version of macOS.). If they also have problems, I'll raise an issue on the sf github page.

# Demonstrate problem with +proj=ortho forward/backward
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.0, PROJ 7.2.1
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"
proj0 <- "+proj=longlat +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

^{Created on 2021-07-15 by the reprex package (v2.0.0)}

dankelley · 2021-07-15T16:51:47Z

Hm, rgdal can do this. (I've not checked yet on whether it reliably gives NAs outside the earth ... that's what I need.)

library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.5-23, (SVN revision 1121)
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
#> Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/gdal
#> GDAL binary built with GEOS: TRUE 
#> Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
#> Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
#> Linking to sp version:1.4-5
#> To mute warnings of possible GDAL/OSR exportToProj4() degradation,
#> use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- project(lonlat, proj)
#> Warning in project(lonlat, proj): 5318 projected point(s) not finite
LONLAT <- project(xy, proj, inv=TRUE)
#> Warning in project(bb, proj, inv = inv, use_aoi = FALSE): 2 projected point(s)
#> not finite
#> Warning in project(xy, proj, inv = TRUE): 5318 projected point(s) not finite
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

^{Created on 2021-07-15 by the reprex package (v2.0.0)}

clayton33 · 2021-07-15T17:48:11Z

Results from my test :

# Demonstrate problem with +proj=ortho forward/backward
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
packageVersion("sf")
#> [1] '0.9.5'
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"
proj0 <- "+proj=longlat +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(c(proj_from_crs(from), proj_from_crs(to)),
#> as.matrix(pts), : one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18363)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.1252 
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] sf_0.9-5
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.5         rstudioapi_0.11    knitr_1.30         magrittr_2.0.1    
#>  [5] units_0.6-7        tidyselect_1.1.0   R6_2.4.1           rlang_0.4.11      
#>  [9] dplyr_1.0.0        stringr_1.4.0      styler_1.4.1       highr_0.8         
#> [13] tools_4.0.2        grid_4.0.2         xfun_0.18          KernSmooth_2.23-17
#> [17] e1071_1.7-3        DBI_1.1.0          withr_2.3.0        class_7.3-17      
#> [21] htmltools_0.5.0    ellipsis_0.3.1     yaml_2.2.1         digest_0.6.26     
#> [25] tibble_3.0.3       lifecycle_0.2.0    crayon_1.3.4       purrr_0.3.4       
#> [29] vctrs_0.3.2        fs_1.5.0           glue_1.4.1         evaluate_0.14     
#> [33] rmarkdown_2.4      reprex_2.0.0       stringi_1.5.3      compiler_4.0.2    
#> [37] pillar_1.4.6       generics_0.0.2     backports_1.1.10   classInt_0.4-3    
#> [41] pkgconfig_2.0.3

^{Created on 2021-07-15 by the reprex package (v2.0.0)}

dankelley · 2021-07-15T17:50:03Z

Any chance of adding a sessionInfo() note? (PS. I submitted an issue to r-spatial/sf and perhaps someone there will want to know this info.)

clayton33 · 2021-07-15T17:56:16Z

updated my comment above, which now includes sessionInfo()

dankelley · 2021-07-15T20:02:23Z

@clayton33 can you do sf::sf_extSoftVersion() and post the results at r-spatial/sf#1729? That might help in the analysis of how things are different. Many thanks. Dan.

richardsc · 2021-07-16T01:11:33Z

Here's my test. I'm actually not sure if it worked or not. Edited to add the results of sf::sf_extSoftVersion() at the bottom:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
packageVersion("sf")
#> [1] '0.9.8'
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"
proj0 <- "+proj=longlat +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] sf_0.9-8
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.6         compiler_4.1.0     pillar_1.6.1       highr_0.9         
#>  [5] class_7.3-19       tools_4.1.0        digest_0.6.27      tibble_3.1.2      
#>  [9] evaluate_0.14      lifecycle_1.0.0    pkgconfig_2.0.3    rlang_0.4.11      
#> [13] reprex_2.0.0       DBI_1.1.1          cli_2.5.0          rstudioapi_0.13   
#> [17] yaml_2.2.1         xfun_0.23          e1071_1.7-6        withr_2.4.2       
#> [21] stringr_1.4.0      dplyr_1.0.6        knitr_1.33         generics_0.1.0    
#> [25] fs_1.5.0           vctrs_0.3.8        tidyselect_1.1.1   classInt_0.4-3    
#> [29] grid_4.1.0         glue_1.4.2         R6_2.5.0           fansi_0.4.2       
#> [33] rmarkdown_2.8      purrr_0.3.4        magrittr_2.0.1     ps_1.6.0          
#> [37] htmltools_0.5.1.1  ellipsis_0.3.2     units_0.7-1        assertthat_0.2.1  
#> [41] KernSmooth_2.23-20 utf8_1.2.1         stringi_1.6.2      proxy_0.4-25      
#> [45] crayon_1.4.1

^{Created on 2021-07-15 by the reprex package (v2.0.0)}

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H
       "3.8.0"        "3.0.4"        "7.0.0"         "true"         "true"

dankelley · 2021-07-16T08:50:50Z

Wow, @richardsc is getting very different results. I will be summarizing all the results in a text file, for transmission to the "sf" team. Thanks!

dankelley · 2021-07-16T09:10:47Z

I've put notes in the sandbox/dk/eow/README.md file.

j-harbin · 2021-07-16T10:04:44Z

Results from my test:

# Demonstrate problem with +proj=ortho forward/backward
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
#> Linking to GEOS 3.9.1, GDAL 3.3.0, PROJ 7.2.1
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"
proj0 <- "+proj=longlat +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.4
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] sf_1.0-0
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.7         compiler_4.1.0     pillar_1.6.1       highr_0.9         
#>  [5] class_7.3-19       tools_4.1.0        digest_0.6.27      tibble_3.1.2      
#>  [9] evaluate_0.14      lifecycle_1.0.0    pkgconfig_2.0.3    rlang_0.4.11      
#> [13] reprex_2.0.0       DBI_1.1.1          cli_3.0.0          rstudioapi_0.13   
#> [17] yaml_2.2.1         xfun_0.24          e1071_1.7-7        withr_2.4.2       
#> [21] stringr_1.4.0      dplyr_1.0.7        knitr_1.33         generics_0.1.0    
#> [25] fs_1.5.0           vctrs_0.3.8        tidyselect_1.1.1   classInt_0.4-3    
#> [29] grid_4.1.0         glue_1.4.2         R6_2.5.0           fansi_0.5.0       
#> [33] rmarkdown_2.9      purrr_0.3.4        magrittr_2.0.1     htmltools_0.5.1.1 
#> [37] ellipsis_0.3.2     units_0.7-2        KernSmooth_2.23-20 utf8_1.2.1        
#> [41] stringi_1.6.2      proxy_0.4-26       crayon_1.4.1

^{Created on 2021-07-16 by the reprex package (v2.0.0)}

sf::sf_extSoftVersion()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
#>        "3.8.1"        "3.2.1"        "7.2.1"         "true"         "true" 
#>           PROJ 
#>        "7.2.1"

^{Created on 2021-07-16 by the reprex package (v2.0.0)}

j-harbin · 2021-07-16T10:35:09Z

I updated my comment to include the output from sf::sf_extSoftVersion()

richardsc · 2021-07-16T18:00:28Z

For completeness I'm going to update my Ubuntu sf (and s2 and wk) packages and re-run the test. It's just a long process because they all have to compile (no binaries for Linux).

paleolimbot · 2021-07-16T18:13:28Z

You should set up RStudio package manager public instance! (Pre compiled binaries for common Linux distributions...)!

richardsc · 2021-07-16T18:27:30Z

Maybe a good idea, except that I can't run RStudio inside the WSL Ubuntu ...

s2 takes about 100 years to compile on my machine. I'm going to harass the package author for doing such a terrible job.

paleolimbot · 2021-07-16T18:54:55Z

No need for RStudio, just set the option(repos = "https://packagemanager.rstudio.com/cran/__linux__/bionic/latest"). See https://packagemanager.rstudio.com/client/#/repos/1/overview

It might be worth checking the output of

sf::sf_proj_pipelines(
  sf::st_crs("+proj=ortho +lat_0=-20"),
  sf::st_crs("+proj=longlat")
)
#> Candidate coordinate operations found:  1 
#> Strict containment:     FALSE 
#> Axis order auth compl:  FALSE 
#> Source:  PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]],
#>             ID["EPSG",6326]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8901]]],
#>     CONVERSION["unknown",
#>         METHOD["Orthographic",
#>             ID["EPSG",9840]],
#>         PARAMETER["Latitude of natural origin",-20,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["False easting",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1,
#>                 ID["EPSG",9001]]]] 
#> Target:  GEOGCRS["unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433],
#>         ID["EPSG",8901]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433,
#>                 ID["EPSG",9122]]]] 
#> Best instantiable operation has accuracy: 0 m
#> Description: Inverse of unknown + axis order change (2D) 
#> Definition:  +proj=pipeline +step +inv +proj=ortho +lat_0=-20 +lon_0=0
#>              +x_0=0 +y_0=0 +ellps=WGS84 +step +proj=unitconvert
#>              +xy_in=rad +xy_out=deg

...to see if a different operation is getting used.

(FWIW, I get the problematic output described by Dan originally on MacOS)

richardsc · 2021-07-16T18:58:39Z

After sf/s2/wk upgrade everything looks fine now:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
packageVersion("sf")
#> [1] '1.0.1'
#> [1] '0.9.8'
proj <- "+proj=ortho +lat_0=-20 +datum=WGS84 +no_defs"
proj0 <- "+proj=longlat +datum=WGS84 +no_defs"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)
ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

^{Created on 2021-07-16 by the reprex package (v2.0.0)}

> sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sf_1.0-1     reprex_2.0.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         pillar_1.6.1       compiler_4.1.0     highr_0.9
 [5] class_7.3-19       tools_4.1.0        digest_0.6.27      evaluate_0.14
 [9] lifecycle_1.0.0    tibble_3.1.2       pkgconfig_2.0.3    rlang_0.4.11
[13] DBI_1.1.1          cli_2.5.0          rstudioapi_0.13    yaml_2.2.1
[17] xfun_0.23          e1071_1.7-6        withr_2.4.2        dplyr_1.0.6
[21] knitr_1.33         generics_0.1.0     fs_1.5.0           vctrs_0.3.8
[25] classInt_0.4-3     grid_4.1.0         tidyselect_1.1.1   glue_1.4.2
[29] R6_2.5.0           processx_3.5.2     fansi_0.4.2        rmarkdown_2.8
[33] callr_3.7.0        purrr_0.3.4        clipr_0.7.1        magrittr_2.0.1
[37] htmltools_0.5.1.1  ps_1.6.0           ellipsis_0.3.2     units_0.7-1
[41] assertthat_0.2.1   utf8_1.2.1         KernSmooth_2.23-20 proxy_0.4-25
[45] crayon_1.4.1

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H
       "3.8.0"        "3.0.4"        "7.0.0"         "true"         "true"
          PROJ
       "7.0.0"

paleolimbot · 2021-07-16T19:02:28Z

The problem is almost 100% the "new" ellipsoidal formulation of the orthographic projection. You can force the spherical version by using +R=6378137 +f=0. See https://proj.org/operations/projections/ortho.html#id1 . This changed in 7.2.0 (PROJ), which is why Clark isn't seeing it (he's on the dinosaur ubuntu bionic).

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1

proj <- "+proj=ortho +lat_0=-20 +R=6378137 +f=0"
proj0 <- "OGC:CRS84"

grid <- expand.grid(lon=seq(-180,180,2.5), lat=seq(-90,90,2.5))
lonlat <- cbind(grid$lon, grid$lat)
xy <- sf_project(proj0, proj, lonlat, keep=TRUE)
#> Warning in CPL_proj_direct(from_to, as.matrix(pts), keep, warn,
#> authority_compliant): one or more projected point(s) not finite
LONLAT <- sf_project(proj, proj0, xy, keep=TRUE)

ok <- is.finite(LONLAT[,1]) & is.finite(LONLAT[,2])

par(mar=c(3,3,1,1), mgp=c(2,0.7,0))
plot(xy[,1]/1e3, xy[,2]/1e3, asp=1,
     xlab="Easting Coordinate [km]", ylab="Northing Coordinate [km]",
     pch=20, cex=0.3, col=ifelse(ok, "red", "blue"))
grid()

^{Created on 2021-07-16 by the reprex package (v0.3.0)}

dankelley · 2021-07-16T19:12:19Z

Re @paleolimbot's "After sf/s2/wk upgrade everything looks fine now:" comment, what are the steps to do these things? Upgrade from CRAN binary/source, or from github/other source?

I'm concerned for students trying to do this sort of thing in September, and I'd love to have clear instructions for them, to avoid the problem. This projection is great for teaching, because I can show things about say the North Atlantic, with the Pacific being hidden, and then I can pop over to the Pacific and they can see that the ideas carry over -- that reasoning and a search for underlying principles is useful in physical oceanography.

richardsc · 2021-07-16T19:39:32Z

Actually that was my comment about the sf/s2/wk upgrade and all it was was an install.packages(‘sf’) call.

But Dewey’s other point about my ubuntu proj libraries being out of date is probably more relevant. I’ll look later tonight.

stale · 2022-04-17T04:58:18Z

An automated message from the stale bot -- this issue has been marked as stale because no comments have been made in 30 days. The goal of this action is to remind reporter and developers alike that progress seems to have stalled. Please note that the issue will be automatically closed 30 days from now, unless a comment is made, or the pinned label is applied. Thank you.

dankelley · 2022-04-17T09:55:15Z

Note to @richardsc -- I added the 'pinned' label to stop the stalebot from closing this in a month. I think something is whacky with the stalebot because this morning I got something like 6 notifications at once about stale things. So I pinned them all, and also made it 100 days until the first message from the bot.

dankelley added the development label Jul 15, 2021

dankelley self-assigned this Jul 15, 2021

dankelley mentioned this issue Jul 15, 2021

problem with inverse projection of "ortho" r-spatial/sf#1729

Closed

dankelley mentioned this issue Jul 16, 2021

should oce use rgdal::spTransform()? #1849

Closed

dankelley added the mapping label Jul 18, 2021

stale bot added the stale label Apr 17, 2022

dankelley added the pinned Prevent automatic stale designation label Apr 17, 2022

stale bot removed the stale label Apr 17, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

EOW (development: edge-of-world) #1848

EOW (development: edge-of-world) #1848

dankelley commented Jul 15, 2021 •

edited

Loading

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

clayton33 commented Jul 15, 2021 •

edited

Loading

dankelley commented Jul 15, 2021

clayton33 commented Jul 15, 2021

dankelley commented Jul 15, 2021

richardsc commented Jul 16, 2021 •

edited

Loading

dankelley commented Jul 16, 2021

dankelley commented Jul 16, 2021 •

edited

Loading

j-harbin commented Jul 16, 2021 •

edited

Loading

j-harbin commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

dankelley commented Jul 16, 2021

richardsc commented Jul 16, 2021

stale bot commented Apr 17, 2022

dankelley commented Apr 17, 2022

EOW (development: edge-of-world) #1848

EOW (development: edge-of-world) #1848

Comments

dankelley commented Jul 15, 2021 • edited Loading

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

dankelley commented Jul 15, 2021

clayton33 commented Jul 15, 2021 • edited Loading

dankelley commented Jul 15, 2021

clayton33 commented Jul 15, 2021

dankelley commented Jul 15, 2021

richardsc commented Jul 16, 2021 • edited Loading

dankelley commented Jul 16, 2021

dankelley commented Jul 16, 2021 • edited Loading

j-harbin commented Jul 16, 2021 • edited Loading

j-harbin commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

richardsc commented Jul 16, 2021

paleolimbot commented Jul 16, 2021

dankelley commented Jul 16, 2021

richardsc commented Jul 16, 2021

stale bot commented Apr 17, 2022

dankelley commented Apr 17, 2022

dankelley commented Jul 15, 2021 •

edited

Loading

clayton33 commented Jul 15, 2021 •

edited

Loading

richardsc commented Jul 16, 2021 •

edited

Loading

dankelley commented Jul 16, 2021 •

edited

Loading

j-harbin commented Jul 16, 2021 •

edited

Loading