Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
dankelley opened this issue Jul 15, 2021 · 25 comments
Open

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

dankelley opened this issue Jul 15, 2021 · 25 comments
Assignees
Labels
development mapping pinned Prevent automatic stale designation
Projects

Comments

@dankelley
Copy link
Owner

dankelley commented Jul 15, 2021

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 dankelley self-assigned this Jul 15, 2021
@github-actions github-actions bot added this to Needs triage in Issues Jul 15, 2021
@dankelley
Copy link
Owner Author

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).

Screen Shot 2021-07-15 at 10 15 55 AM

@dankelley
Copy link
Owner Author

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.)

Screen Shot 2021-07-15 at 10 28 33 AM

@dankelley
Copy link
Owner Author

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?

Screen Shot 2021-07-15 at 11 18 45 AM

@dankelley
Copy link
Owner Author

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)

Screen Shot 2021-07-15 at 11 31 33 AM

@dankelley
Copy link
Owner Author

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
Copy link
Owner Author

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
Copy link
Collaborator

clayton33 commented Jul 15, 2021

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
Copy link
Owner Author

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
Copy link
Collaborator

updated my comment above, which now includes sessionInfo()

@dankelley
Copy link
Owner Author

@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
Copy link
Collaborator

richardsc commented Jul 16, 2021

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
Copy link
Owner Author

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
Copy link
Owner Author

dankelley commented Jul 16, 2021

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

@j-harbin
Copy link
Contributor

j-harbin commented Jul 16, 2021

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
Copy link
Contributor

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

@dankelley dankelley moved this from Needs triage to High priority in Issues Jul 16, 2021
@richardsc
Copy link
Collaborator

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
Copy link

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

@richardsc
Copy link
Collaborator

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
Copy link

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
Copy link
Collaborator

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
Copy link

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
Copy link
Owner Author

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
Copy link
Collaborator

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.

@dankelley dankelley moved this from High priority to Low priority in Issues Aug 7, 2021
@stale
Copy link

stale bot commented Apr 17, 2022

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.

@stale stale bot added the stale label Apr 17, 2022
@dankelley dankelley added the pinned Prevent automatic stale designation label Apr 17, 2022
@stale stale bot removed the stale label Apr 17, 2022
@dankelley
Copy link
Owner Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
development mapping pinned Prevent automatic stale designation
Projects
Issues
Low priority
Development

No branches or pull requests

5 participants