Mollweide projection ignores axis breaks in geom_sf()

[fzenoni]

When using a Mollweide projection (crs = 54009), ggplot2 ignores the axis breaks given by scale_x|y_continuous(). See reprex and output.

library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(ggplot2)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

grid <- st_sf(st_make_grid(what = 'polygons'))

ggplot(grid) +
  geom_sf() +
  scale_y_continuous(breaks = c(-60, -30, 0, 30, 60)) +
  coord_sf(crs = st_crs(54009)) + # Mollweide
  #coord_sf(crs = st_crs(54030)) + # Robinson
  theme(
    panel.background = element_blank(),
    panel.ontop = TRUE,
    panel.grid.major = element_line(color = 'red',
                                    size = 1))

Created on 2020-04-11 by the reprex package (v0.3.0)

Comparison with a Robinson projection.

[clauswilke]

I think this is somehow related to the bounding box that is being sent to sf::st_graticule(). Pinging @edzer so he can provide input. What follows is my investigation so far. It's incomplete, but I want to record what I have figured out up to this point.

Update: The info in this paragraph is not correct. Ignore. First, I'm pretty certain that in the current version of coord_sf() breaks need to be specified in projected coordinates, not in long/lat. However, this is not properly documented in sf::st_graticule() and so I'm somewhat confused. I'm also confused why it seems to work for the Robinson projection. @edzer, any input on this?

Second, the graticule calculation seems to depend on the bounding box provided, and if the bounding box is too large then st_graticule() seems to ignore the long and lat parameters.

library(sf)
#> Warning: package 'sf' was built under R version 3.6.2
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(ggplot2)

grid <- st_sf(st_make_grid(what = 'polygons'))

p <- ggplot(grid) + geom_sf() +
  scale_x_continuous(breaks = 1e6*c(-10, -5, 0, 5, 10)) +
  coord_sf(crs = st_crs(54009)) # Mollweide

b <- ggplot_build(p)
xbreaks <- b$layout$get_scales(1)$x$breaks
xbreaks
#> [1] -1e+07 -5e+06  0e+00  5e+06  1e+07
ybreaks <- b$layout$get_scales(1)$y$breaks
ybreaks
#> list()
#> attr(,"class")
#> [1] "waiver"

# this is the bounding box used by ggplot2 when it calculates the graticule
bbox <- c(-19844105, -9922053, 19844105, 9922053)

# this is what the call to st_graticule() produces. Too many
# graticule lines of type "E"
sf::st_graticule(
  bbox,
  crs = st_crs(54009),
  lat = NULL,
  lon = xbreaks,
  datum = st_crs(4326),
  ndiscr = 100
)
#> Simple feature collection with 28 features and 10 fields
#> Attribute-geometry relationship: 2 constant, 0 aggregate, 0 identity, 8 NA's
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: -18040100 ymin: -9020048 xmax: 18040100 ymax: 9020048
#> CRS:            EPSG:54009
#> First 10 features:
#>    degree type degree_label                       geometry x_start y_start
#> 1    -180    E         <NA> LINESTRING (-1.104637e-09 -...    -180  -89.91
#> 2    -160    E         <NA> LINESTRING (-9.818998e-10 -...    -160  -89.91
#> 3    -140    E         <NA> LINESTRING (-8.591623e-10 -...    -140  -89.91
#> 4    -120    E         <NA> LINESTRING (-7.364248e-10 -...    -120  -89.91
#> 5    -100    E         <NA> LINESTRING (-6.136874e-10 -...    -100  -89.91
#> 6     -80    E         <NA> LINESTRING (-4.909499e-10 -...     -80  -89.91
#> 7     -60    E         <NA> LINESTRING (-3.682124e-10 -...     -60  -89.91
#> 8     -40    E         <NA> LINESTRING (-2.454749e-10 -...     -40  -89.91
#> 9     -20    E         <NA> LINESTRING (-1.227375e-10 -...     -20  -89.91
#> 10      0    E         <NA> LINESTRING (0 -9020048, 0 -...       0  -89.91
#>    x_end y_end angle_start angle_end plot12
#> 1   -180 89.91          90        90   TRUE
#> 2   -160 89.91          90        90   TRUE
#> 3   -140 89.91          90        90   TRUE
#> 4   -120 89.91          90        90   TRUE
#> 5   -100 89.91          90        90   TRUE
#> 6    -80 89.91          90        90   TRUE
#> 7    -60 89.91          90        90   TRUE
#> 8    -40 89.91          90        90   TRUE
#> 9    -20 89.91          90        90   TRUE
#> 10     0 89.91          90        90   TRUE

# this is the graticule stored by ggplot2, which looks the same
b$layout$panel_params[[1]]$graticule
#> Simple feature collection with 28 features and 10 fields
#> Attribute-geometry relationship: 2 constant, 0 aggregate, 0 identity, 8 NA's
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 0.04545455 ymin: 0.04545455 xmax: 0.9545455 ymax: 0.9545455
#> CRS:            NA
#> First 10 features:
#>    degree type degree_label                       geometry   x_start   y_start
#> 1    -180    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999955 0.4999955
#> 2    -160    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999960 0.4999955
#> 3    -140    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999965 0.4999955
#> 4    -120    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999970 0.4999955
#> 5    -100    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999975 0.4999955
#> 6     -80    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999980 0.4999955
#> 7     -60    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999985 0.4999955
#> 8     -40    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999990 0.4999955
#> 9     -20    E         <NA> LINESTRING (0.5 0.04545455,... 0.4999995 0.4999955
#> 10      0    E         <NA> LINESTRING (0.5 0.04545455,... 0.5000000 0.4999955
#>        x_end     y_end angle_start angle_end plot12
#> 1  0.4999955 0.5000045          90        90   TRUE
#> 2  0.4999960 0.5000045          90        90   TRUE
#> 3  0.4999965 0.5000045          90        90   TRUE
#> 4  0.4999970 0.5000045          90        90   TRUE
#> 5  0.4999975 0.5000045          90        90   TRUE
#> 6  0.4999980 0.5000045          90        90   TRUE
#> 7  0.4999985 0.5000045          90        90   TRUE
#> 8  0.4999990 0.5000045          90        90   TRUE
#> 9  0.4999995 0.5000045          90        90   TRUE
#> 10 0.5000000 0.5000045          90        90   TRUE

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

[clauswilke]

Sorry, I got confused about how st_graticule() works. It does indeed take limits in longitude and latitude, not in projected coordinates. The problem seems to be that for the Mollweide projection this doesn't work when the bounding box contains the entire globe. If we limit the bounding box to a smaller area, things work as expected:

# approx. 22x22 degree box around 0, 0
bbox <- c(-1100000, -1350000, 1100000, 1350000)
sf::st_graticule(
  bbox,
  crs = sf::st_crs(54009), # Mollweide
  lat = c(-5, -3, 0, 3, 5),
  lon = c(-5, -3, 0, 3, 5),
  datum = sf::st_crs(4326),
  ndiscr = 100
)
#> Simple feature collection with 10 features and 10 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: -1100000 ymin: -1350000 xmax: 1100000 ymax: 1350000
#> CRS:            EPSG:54009
#>    degree type degree_label                       geometry    x_start
#> 1      -5    E "5"*degree*W LINESTRING (-495469.5 -1350...  -495469.5
#> 2      -3    E "3"*degree*W LINESTRING (-297281.7 -1350...  -297281.7
#> 3       0    E   "0"*degree LINESTRING (0 -1350000, 0 -...        0.0
#> 4       3    E "3"*degree*E LINESTRING (297281.7 -13500...   297281.7
#> 5       5    E "5"*degree*E LINESTRING (495469.5 -13500...   495469.5
#> 6      -5    N "5"*degree*S LINESTRING (-1100000 -61792... -1100000.0
#> 7      -3    N "3"*degree*S LINESTRING (-1100000 -37086... -1100000.0
#> 8       0    N   "0"*degree LINESTRING (-1100000 0, -10... -1100000.0
#> 9       3    N "3"*degree*N LINESTRING (-1100000 370869... -1100000.0
#> 10      5    N "5"*degree*N LINESTRING (-1100000 617923... -1100000.0
#>       y_start     x_end     y_end angle_start angle_end plot12
#> 1  -1350000.0 -495469.5 1350000.0    90.47687  89.52313   TRUE
#> 2  -1350000.0 -297281.7 1350000.0    90.28612  89.71388   TRUE
#> 3  -1350000.0       0.0 1350000.0    90.00000  90.00000   TRUE
#> 4  -1350000.0  297281.7 1350000.0    89.71388  90.28612   TRUE
#> 5  -1350000.0  495469.5 1350000.0    89.52313  90.47687   TRUE
#> 6   -617923.6 1100000.0 -617923.6     0.00000   0.00000   TRUE
#> 7   -370869.7 1100000.0 -370869.7     0.00000   0.00000   TRUE
#> 8         0.0 1100000.0       0.0     0.00000   0.00000   TRUE
#> 9    370869.7 1100000.0  370869.7     0.00000   0.00000   TRUE
#> 10   617923.6 1100000.0  617923.6     0.00000   0.00000   TRUE

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

Similarly:

library(sf)
#> Warning: package 'sf' was built under R version 3.6.2
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(ggplot2)

grid <- st_sf(st_make_grid(what = 'polygons'))

ggplot(grid) +
  geom_sf() +
  scale_x_continuous(breaks = c(-5, -3, 0, 3, 5)) +
  scale_y_continuous(breaks = c(-5, -3, 0, 3, 5)) +
  coord_sf(
    crs = st_crs(54009), # Mollweide
    xlim = c(-1100000, 1100000),
    ylim = c(-1350000, 1350000)
  ) + 
  theme(
    panel.background = element_blank(),
    panel.ontop = TRUE,
    panel.grid.major = element_line(color = 'red', size = 1)
  )

^{Created on 2020-04-11 by the reprex package (v0.3.0)}

So it's not clear to me that there's anything to be done in ggplot2. It's an sf problem. At a minimum, st_graticule() should issue a warning if it ignores the values of the lat and long arguments.

[clauswilke]

I have filed an issue with the sf project: r-spatial/sf#1361

[thomasp85]

closing here