Instantly share code, notes, and snippets.

# fzenoni/ortho_view.R

Last active November 2, 2023 11:06
Show Gist options
• Save fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1 to your computer and use it in GitHub Desktop.
Preserve polygons in orthographic view
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # Download earth data first # https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip library(sf) library(lwgeom) library(dplyr) library(ggplot2) library(mapview) # Read the data mini_world <- read_sf('data/ne_110m_land/ne_110m_land.shp') # Define the orthographic projection # Choose lat_0 with -90 <= lat_0 <= 90 and lon_0 with -180 <= lon_0 <= 180 lat <- 45 lon <- 2 ortho <- paste0('+proj=ortho +lat_0=', lat, ' +lon_0=', lon, ' +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs') # Define the polygon that will help you finding the "blade" # to split what lies within and without your projection circle <- st_point(x = c(0,0)) %>% st_buffer(dist = 6371000) %>% st_sfc(crs = ortho) # Project this polygon in lat-lon circle_longlat <- circle %>% st_transform(crs = 4326) # circle_longlat cannot be used as it is # You must decompose it into a string with ordered longitudes # Then complete the polygon definition to cover the hemisphere if(lat != 0) { circle_longlat <- st_boundary(circle_longlat) circle_coords <- st_coordinates(circle_longlat)[, c(1,2)] circle_coords <- circle_coords[order(circle_coords[, 1]),] circle_coords <- circle_coords[!duplicated(circle_coords),] # Rebuild line circle_longlat <- st_linestring(circle_coords) %>% st_sfc(crs = 4326) if(lat > 0) { rectangle <- list(rbind(circle_coords, c(X = 180, circle_coords[nrow(circle_coords), 'Y']), c(X = 180, Y = 90), c(X = -180, Y = 90), c(X = -180, circle_coords[1, 'Y']), circle_coords[1, c('X','Y')])) %>% st_polygon() %>% st_sfc(crs = 4326) } else { rectangle <- list(rbind(circle_coords, c(X = 180, circle_coords[nrow(circle_coords), 'Y']), c(X = 180, Y = -90), c(X = -180, Y = -90), c(X = -180, circle_coords[1, 'Y']), circle_coords[1, c('X','Y')])) %>% st_polygon() %>% st_sfc(crs = 4326) } circle_longlat <- st_union(st_make_valid(circle_longlat), st_make_valid(rectangle)) } # This visualization shows the visible emisphere in red ggplot() + geom_sf(data = mini_world) + geom_sf(data = circle_longlat, color = 'red', fill = 'red', alpha = 0.3) # A small negative buffer is necessary to avoid polygons still disappearing in a few pathological cases # I should not change the shapes too much visible <- st_intersection(st_make_valid(mini_world), st_buffer(circle_longlat, -0.09)) %>% st_transform(crs = ortho) # DISCLAIMER: This section is the outcome of trial-and-error and I don't claim it is the best approach # Resulting polygons are often broken and they need to be fixed # Get reason why they're broken broken_reason <- st_is_valid(visible, reason = TRUE) # First fix NA's by decomposing them # Remove them from visible for now na_visible <- visible[is.na(broken_reason),] visible <- visible[!is.na(broken_reason),] # Open and close polygons na_visible <- st_cast(na_visible, 'MULTILINESTRING') %>% st_cast('LINESTRING', do_split=TRUE) na_visible <- na_visible %>% mutate(npts = npts(geometry, by_feature = TRUE)) # Exclude polygons with less than 4 points na_visible <- na_visible %>% filter(npts >=4) %>% select(-npts) %>% st_cast('POLYGON') # Fix other broken polygons broken <- which(!st_is_valid(visible)) for(land in broken) { result = tryCatch({ # visible[land,] <- st_buffer(visible[land,], 0) # Sometimes useful sometimes not visible[land,] <- st_make_valid(visible[land,]) %>% st_collection_extract() }, error = function(e) { visible[land,] <<- st_buffer(visible[land,], 0) }) } # Bind together the two tables visible <- rbind(visible, na_visible) # Final plot ggplot() + geom_sf(data = circle, #fill = 'aliceblue') + # if you like the color fill = NA) + geom_sf(data=st_collection_extract(visible)) + coord_sf(crs = ortho)

Hi @fzenoni! Wow! I cannot thank you enough for sharing this! I won't get a chance to dig into it until next week, but I'll be sure to drop any questions and results here as soon as I do. Very much appreciated!

### csaybar commented Aug 27, 2019

Thanks for sharing @fzenoni, it is simply brilliant!

Your plots in the sf issue look great! So I tried to reproduce them.

The first warning appeared after line 58:
`although coordinates are longitude/latitude, st_union assumes that they are planar`

The following plot looked as expected. But i wasn't able to create the object `visible`. Lines 67-68 resulted in the following error

``````dist is assumed to be in decimal degrees (arc_degrees).
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) :
Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point -132.71000788443121 54.040009315423447 at -132.71000788443121 54.040009315423447.
In st_buffer.sfc(circle_longlat, -0.09) :
st_buffer does not correctly buffer longitude/latitude data
``````

I'm very grateful for any advice how to get this working.

Session info
``````> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 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=de_DE.UTF-8       LC_NUMERIC=C
[3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8
[5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8
[7] LC_PAPER=de_DE.UTF-8       LC_NAME=C
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] mapview_2.7.0 ggplot2_3.2.1 dplyr_0.8.3   lwgeom_0.1-7  sf_0.8-0

loaded via a namespace (and not attached):
[1] Rcpp_1.0.3         later_1.0.0        pillar_1.4.2       compiler_3.6.1
[5] base64enc_0.1-3    class_7.3-15       tools_3.6.1        digest_0.6.23
[9] viridisLite_0.3.0  satellite_1.0.2    lifecycle_0.1.0    tibble_2.1.3
[13] gtable_0.3.0       lattice_0.20-38    png_0.1-7          pkgconfig_2.0.3
[17] rlang_0.4.2        shiny_1.4.0        DBI_1.0.0          crosstalk_1.0.0
[21] fastmap_1.0.1      e1071_1.7-3        raster_3.0-7       withr_2.1.2
[25] htmlwidgets_1.5.1  webshot_0.5.2      stats4_3.6.1       classInt_0.4-2
[29] leaflet_2.0.3      grid_3.6.1         tidyselect_0.2.5   glue_1.3.1
[33] R6_2.4.1           sp_1.3-2           farver_2.0.1       purrr_0.3.3
[37] magrittr_1.5       codetools_0.2-16   promises_1.1.0     scales_1.1.0
[41] htmltools_0.4.0    units_0.6-5        assertthat_0.2.1   xtable_1.8-4
[45] mime_0.7           colorspace_1.4-1   httpuv_1.5.2       KernSmooth_2.23-16
[49] lazyeval_0.2.2     munsell_0.5.0      crayon_1.3.4
``````
Geo Libraries
``````> library(sf)
Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
> library(lwgeom)
Linking to liblwgeom 2.5.0dev r16016, GEOS 3.8.0, PROJ 6.2.1
``````

### fzenoni commented Dec 9, 2019

Hi @rnuske, thank you for testing the script.

Don't worry about the warning, it is a standard print when you manipulate or intersect geometries in lat/lon coordinates. To avoid it, in principle one should:

1. project to some planar CRS
2. do the operation
3. project back to lat/lon.

Unless you are looking after precise results (e.g. compute the centroid of a region with great accuracy), usually it is not worth it.

Concerning the error, I was able to reproduce it.
I have no explanation about why it used to work before, but I have now replaced line 67 with

``````visible <- st_intersection(st_make_valid(mini_world), st_buffer(circle_longlat, -0.09)) %>%
``````

That should do it.

### rnuske commented Dec 9, 2019

With the new line 67 it also works for me.

### rnuske commented Dec 9, 2019

just a short follow up question: Do you happen to know how to color the oceans?

I am very sorry I missed your message.
To display the color you may uncomment https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1#file-ortho_view-r-L108 and comment https://gist.github.com/fzenoni/ef23faf6d1ada5e4a91c9ef23b0ba2c1#file-ortho_view-r-L109. That will completely hide the graticule, though.

To preserve (and control) the graticule, you need to use the `ggplot2`'s `theme()` function.

This is inspired by the code I showed in the original Github issue.

``````ggplot() +
geom_sf(data = circle,
fill = 'aliceblue') +
geom_sf(data=st_collection_extract(visible)) +
coord_sf(crs = ortho) +
theme(
panel.grid.major = element_line(
color = 'red', linetype = 2, size = 0.5),
panel.ontop = TRUE,
panel.background = element_rect(fill = NA)
)
``````

You might also want to play with the alpha level, so that the graticule goes below the land masses.
It's a quirky solution because one can't perfectly control the color of the sea anymore, nor the color of the graticule.

``````# Final plot
ggplot() +
geom_sf(data = circle,
fill = 'aliceblue', alpha=0.5)+
geom_sf(data=st_collection_extract(visible), fill='darkgreen', color = 'black') +
coord_sf(crs = ortho) +
theme(
panel.grid.major = element_line(
color = 'red', linetype = 2, size = 0.5)
)
``````

### rnuske commented Dec 17, 2019

Thanks a lot for your suggestions @fzenoni. They helped a lot to get what I wanted. I should definitely learn more ggplot!

Did I understand correctly that the graticule can only be behind everything else (in the background) or in front of everything (in the foreground) but not between the layers?

### fzenoni commented Dec 17, 2019

I am glad it worked out for you!

Did I understand correctly that the graticule can only be behind everything else (in the background) or in front of everything (in the foreground) but not between the layers?

After a lot of trial and error, I believe this is correct. But it would be better to ask this question to the `ggplot2` maintainers.

### jlguilherme commented Mar 11, 2021

This is brilliant! Thanks for sharing!

### rasmerla commented Oct 11, 2021

This is totally brilliant! Thanks a lot! Got this link trough a question a Stack exchange when I tried to solve it on my own.

I haven't really wrapped my head around the effects of using s2 geometry, but I did not get this script to work unless I set sf_use_s2(FALSE).

### fzenoni commented Oct 12, 2021

Thank you @rasmerla, I really appreciate that.

The script was written before `sf` 1.0, and since then my professional activities shifted a bit, so I haven't had the opportunity yet to adapt it to the new `s2` dependency. This vignette gives me the feeling that things got easier though.

### rCarto commented Nov 2, 2023

Here is an example of ortho view using s2:
https://gist.github.com/rCarto/02dff288244691a53b3c7b55ee44ef86
It's probably not bulletproof though.