How to migrate from coord_map() to coord_sf() for orthographic projection?

[gavinsimpson]

After reporting a change in behaviour related to coord_map() and the upcoming changes to the guides system in ggplot2 v 3.5.0 tidyverse/ggplot2#5684 it seems I'll need to move the plotting code in my gratia package to use coord_sf().

The following example, which is standalone but illustrates the main points of the code used in the package, simulates data on a sphere, to which we then fit a GAM to the data using a special basis due to Duchon which {mgcv} calls a spline on the sphere via bs = "sos":

library("mgcv")
#> Loading required package: nlme
#> This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
library("ggplot2")
`sim_sos_eg_data` <- function(n = 400, seed = NULL) {
  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
    runif(1)
  }
  if (is.null(seed)) {
    RNGstate <- get(".Random.seed", envir = .GlobalEnv)
  } else {
    R.seed <- get(".Random.seed", envir = .GlobalEnv)
    set.seed(seed)
    RNGstate <- structure(seed, kind = as.list(RNGkind()))
    on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
  }
  f <- function(la, lo) { ## a test function...
    sin(lo) * cos(la - 0.3)
  }
  ## generate with uniform density on sphere...
  lo <- runif(n) * 2 * pi - pi ## longitude
  la <- runif(3 * n) * pi - pi / 2
  ind <- runif(3 * n) <= cos(la)
  la <- la[ind]
  la <- la[seq_len(n)]
  ff <- f(la, lo)
  y <- ff + rnorm(n) * 0.2 ## test data
  out <-data.frame(latitude = la * 180 / pi,
    longitude = lo * 180 / pi, y = y)
  out
}
sos_df <- sim_sos_eg_data(n = 400, seed = 0)
m_sos <- gam(y ~ s(latitude, longitude, bs = "sos", k = 60), data = sos_df)
new_df <- with(sos_df,
  expand.grid(longitude = seq(min(longitude), max(longitude), length = 15),
    latitude = seq(min(latitude), max(latitude), length = 15)))
# sm <- smooth_estimates(m_sos, n = 25)
p <- predict(m_sos, newdata = new_df)
new_df$.fitted <- p

At this point I have a data frame (tibble) containing the values I want to plot (the latitude and longitude values and the model estimated value of the response. Using coord_map() I get a nice enough plot to visualise this special smoother. If you're using the rc of ggplot 3.5.0 this will lead to some warnings and no axes get plotted but it works just fine in <3.5.0.

new_df |>
    ggplot(aes(x = longitude, y = latitude)) +
    geom_tile(aes(fill = .fitted)) +
    coord_map(projection = "orthographic",
        orientation = c(20, 0, 0)) +
    guides(fill = guide_colourbar(title = "Partial effect",
            direction = "vertical",
            barheight = grid::unit(0.25, "npc")),
        x = guide_axis(angle = 45))

This is (close to) the desired plot

but as I'm using the rc for 3.5.0 I'm not getting the axes/ticks.

I have tried variations on the theme of switching to coord_sf() and setting either the crs or default_crs arguments as implied on the coord_map() help page:

For example, I have tried:

new_df |> 
        ggplot(aes(x = longitude, y = latitude)) + 
        geom_tile(aes(fill = .fitted)) + 
        coord_sf(crs = "+proj=ortho +lat_0=20 +lon_0=-10")

new_df |> 
        ggplot(aes(x = longitude, y = latitude)) + 
        geom_tile(aes(fill = .fitted)) + 
        coord_sf(default_crs = "+proj=ortho +lat_0=20 +lon_0=-10")

new_df |> 
        ggplot(aes(x = longitude, y = latitude)) + 
        geom_tile(aes(fill = .fitted)) + 
        coord_sf(crs = sf::sf_crs(9840))

all to no avail:

so I assume that moving to coord_sf() is not sufficient and that I am missing some other steps with sf to correctly achieve an orthographic projection similar to that achieved with coord_map().

Do I need to do to something to my data frame (tibble) new_df to get coord_sf() to kick in? Or am I missing something fundamental?

Thanks in advance

[edzer]

I get the plot below with ggplot 3.4.4, do you consider that "just fine"?

[gavinsimpson]

get the plot below with ggplot 3.4.4

With coord_map()?

do you consider that "just fine"?

Yes :) It's "just fine" enough for my purposes as an automatic plot of a smooth from a GAM. If a user wants to polish this up then they are welcome to do that and gratia provides all the tools they would need to extract the data for this spline from the model.

If you are referring to the weirdness on the eastern side of the globe, that's not great but it is much less visible if you increase the number of points in the lat/long grid that we predict at; Using a 15x15 grid for the entire globe is not something someone is going to do in practice, but I kept the size low here so you weren't waiting a long time for predict() to run or ggplot() to render the plot.

[edzer]

With coord_map()?

Yes, from using your first two code chunks. Note that I'm not getting the axes/ticks either, so that doesn't appear to me as a ggplot 3.5.0 rc issue. Your issue seems to be with coord_map() and coord_sf(), which are both part of ggplot2, and with changes in ggplot2, so is this the right place to open this / discuss?

As of the weirdness on the eastern side: this is a known issue of sf: if you project a polygon and part of the polygon is "invisible" (as here: on the back-side of the globe), the whole polygon gets discarded, rather than cut back to the visible region. Computing the intersection of the visible area before projecting, as done here would be the way to go, I think a helper function in sf for the orthogonal projection would be relatively easy to realise.

[gavinsimpson]

Sorry, I think we're talking at cross purposes here.

Yes, there are no ticks in 3.4.4 and in 3.5.0, but in 3.5.0 using coord_map() and guides() raises a warning, which CRAN will ask me to fix as soon as the ggplot2 devs release this new version to CRAN. As coord_map() is superseded, the ggplot2 devs aren't going to fix the use of guides() with coord_map() so I need to find an alternative and coord_sf() is the suggested one.

Your issue seems to be with coord_map() and coord_sf(), which are both part of ggplot2, and with changes in ggplot2, so is this the right place to open this / discuss?

That's a good question: I don't know whether I am using coord_sf() incorrectly (in which case this would be the wrong place to raise this, agreed) or whether I need to get my object new_df into some other form suitable for plotting with coord_sf(). That latter point is perhaps more on topic for here?

To be honest, I don't even know if the behaviour of coord_sf() in this example (where I don't see any projection happening) is expected, or a bug with ggplot2, or my complete misunderstanding of how to plot sf objects with ggplot.

As of the weirdness on the eastern side: this is a known issue of sf...

but that image was generated using coord_map() not coord_sf(). I can't get coord_sf() to draw anything that looks like it has been projected to an orthographic projection. Which makes me wonder if I have misunderstood how this should work.

So perhaps to focus my question: do I have to convert new_df to some type of sf object before I plot it and use coord_sf()? If I do, how do I do that conversion, where what I have are the coordinates for the centres of the tiles/pixels in the grid.

[edzer]

Thanks, more clear now. Here's one way with geom_sf():

library(sf)
library(stars)
sf = st_as_stars(new_df) |> st_set_crs('OGC:CRS84') |> st_as_sf(points = FALSE)                                            
sf_o = st_transform(sf, "+proj=ortho +lat_0=20 +lon_0=0")                                                                  
sf_o = sf_o[st_is_valid(sf_o),] 
ggplot() + geom_sf(data = sf_o, aes(fill = .fitted))   

Note that your new_df contains point coordinates, and not polygons; geom_sf needs polygons if you want polygons shown on your plot. stars is used here to create a raster from points, and convert that into polygons. st_is_valid() is used to filter out polygons invisible or half visible (and hence invalid). I'll see if I can cook up a function that does the cutting first.

[edzer]

latitude values in new_df are btw odd: they have regular intervals, but the lowest value is less than half an interval size from -90; st_as_stars() assumes points are grid cell centers, which they can't be true (the bounding box then extends beyond latitude -90):

new_df$latitude |> unique() |> sort() |> diff()
# [1] 12.10202 12.10202 12.10202 12.10202 12.10202 12.10202 12.10202 12.10202 
# [9] 12.10202 12.10202 12.10202 12.10202 12.10202 12.10202
min(new_df$latitude)
# [1] -87.64025 

this explains also the hole on the North pole in all plots above. Ah, I see where this comes from; anyway it messes up the logic if you want to create a polygon coverage for the world.

[gavinsimpson]

latitude values in new_df are btw odd

Yeah; in part that was me just wanting something to illustrate the problem in ggplot in the bug report I filed with them. But it does reflect the current way in which I generate the grid of values at which I evaluate the spline on the sphere. Technically, geom_tile() also assumes that I have correctly specified the cell midpoints and width but it seemed to work "ok" despite my abusing it with data that aren't actually correctly specified.

From your other comment with a 99% solution to my problem, it seems that there isn't a simple way to just switch from coord_map() to coord_sf(), so if I need to now depend on sf and stars I should also spend a little time being less lazy and define the points on the grid properly too.

[edzer]

Here's a round Earth:

new_df <- with(sos_df,
  expand.grid(longitude = seq(-175, 175, 10),                                                                              
    latitude = seq(-85, 85, 10)))
# sm <- smooth_estimates(m_sos, n = 25)                                                                                    
p <- predict(m_sos, newdata = new_df)                                                                                      
new_df$.fitted <- p
sf = st_as_stars(new_df) |> st_set_crs('OGC:CRS84') |> st_as_sf(points = FALSE) |> st_set_agr("constant")                  
st_ortho_cut = function(x, lon_0, lat_0, radius = 9800000) {                                                               
    stopifnot(st_is_longlat(x))
    pt = st_sfc(st_point(c(lon_0, lat_0)), crs = 'OGC:CRS84')                                                              
    buf = st_buffer(pt, units::set_units(radius, m))
    ortho = paste0("+proj=ortho +lat_0=", lat_0, " +lon_0=", lon_0)                                                        
    st_transform(st_intersection(x, buf), st_crs(ortho))                                                                   
}
sf_o = st_ortho_cut(sf, lat_0 = 20, lon_0 = 0)
ggplot() + geom_sf(data = sf_o, aes(fill = .fitted))  

sf::st_make_grid() may be of help if you want to avoid dependence on stars.

[teunbrand]

Using geom_tile() generally works if you specify that non-sf layers use the '4326' i.e. long/lat CRS.
That is one potential way to avoid taking on additional dependencies
(EDIT: 'overt additional dependencies', thanks Edzer ;) You'd be swapping {mapproj} dependency for {sf} dependency).
Admittedly, this doesn't resolve the grid misalignment or hidden polygons issue.

new_df |>
  ggplot(aes(longitude, latitude, fill = .fitted)) +
  geom_tile() +
  coord_sf(crs = "+proj=ortho +lat_0=20 +lon_0=0", default_crs = 4326)

[edzer]

That is one potential way to avoid taking on additional dependencies.

The dependency is there, it's just indirect. ;-)

[gavinsimpson]

Using geom_tile() generally works if you specify that non-sf layers use the '4326' i.e. long/lat CRS.

Ahhh, @teunbrand thanks for this. So the point is that both crs and default_crs need to be used. I'll look at ?coord_map and ?coord_sf again as I didn't pick up on that meaning when I read them, and it wasn't clear to me how these arguments worked.

[gavinsimpson]

Thanks @edzer for the working sf example, that's very neat. I'll need to spend a little time figuring out how to make generating the regular grid of evaluation points work for data that only occupy part of the sphere. Longer term it would seem that moving to explicit dependencies on sf and stars would be the way to go, while using @teunbrand's solution will help me avoid the warning right now when ggplot2 3.5.0 is released and I get the email from CRAN.