coord_sf collapses all points to 0,0

[wchase]

I'm trying to plot some points on a Mollweide projected map, but adding coord_sf to my ggplot code causes all of the points to be plotted at 0,0. Here is some code where I grab the country data, project it to Mollweide (crs=54009), then load my points data. Note that the points data has already been projected to Mollweide, so it should need no reprojection. If I plot the data without using coord_sf, the points show up as they should, but of course the projection looks wrong. But if I add coord_sf, it seems to collapse all of the points to 0,0.

library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)

#country data
countries <- ne_countries(returnclass = "sf")
world_proj <- st_transform(countries, crs = 54009)

#points, already projected to same crs
points <- data.frame(
        land = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,
                 FALSE),
           x = c(-66.89695, 143.28066, 51.63533, -76.16409, 177.77862,
                 171.00648, 101.41291, -147.60461, 117.14617, 107.37859),
           y = c(76.77662, -47.098278, -30.987178, 84.906883, -70.977776,
                 -7.352523, -78.711968, 23.036947, -51.915945, -3.458423)
)

#doesnt work
ggplot() +
  geom_sf(data = countries, fill = "grey80", col = "grey40", lwd = 0.3) +
  geom_point(data = points, aes(x, y, color = land)) +
  theme_minimal() +
  coord_sf(crs = 54009) +
  theme(axis.text = element_blank())

#works for the points, but the projection is wrong now
ggplot() +
  geom_sf(data = countries, fill = "grey80", col = "grey40", lwd = 0.3) +
  geom_point(data = points, aes(x, y, color = land)) +
  theme_minimal() +
  #coord_sf(crs = 54009) +
  theme(axis.text = element_blank())

[jlacko]

On second thought: it seems a minor tweak in datum of coords_sf does the trick

countries <- ne_countries(returnclass = "sf")
world_proj <- st_transform(countries, crs = "+proj=moll")

points <- data.frame(
  land = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,
           FALSE),
  x = c(-66.89695, 143.28066, 51.63533, -76.16409, 177.77862,
        171.00648, 101.41291, -147.60461, 117.14617, 107.37859),
  y = c(76.77662, -47.098278, -30.987178, 84.906883, -70.977776,
        -7.352523, -78.711968, 23.036947, -51.915945, -3.458423)) #

ggplot() +
  geom_sf(data = countries, fill = "grey80", col = "grey40", lwd = 0.3) +
  geom_point(data = points, aes(x, y, color = land)) +
  theme_minimal() +
  coord_sf(datum = st_crs(world_proj)) +
  theme(axis.text = element_blank())

[wchase]

Thanks for the reply! I actually found another bug, I was using the countries data in the plotting code instead of the projected data. But adding the projected data still has the same issue, and now it happens even if I remove the coord_sf. Removing coord_sf should default to whatever projection is in the first layer, which now is crs=54009. It seems like there's some issue with the projection where it is transforming the points in some incorrect way. I also tried transforming to points to an sf object and plotting with geom_sf, but I get the same result.

#country data
countries <- ne_countries(returnclass = "sf")
world_proj <- st_transform(countries, crs = 54009)

#points, already projected to same crs
points <- data.frame(
        land = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,
                 FALSE),
           x = c(-66.89695, 143.28066, 51.63533, -76.16409, 177.77862,
                 171.00648, 101.41291, -147.60461, 117.14617, 107.37859),
           y = c(76.77662, -47.098278, -30.987178, 84.906883, -70.977776,
                 -7.352523, -78.711968, 23.036947, -51.915945, -3.458423)
)

#doesnt work
ggplot() +
  geom_sf(data = world_proj, fill = "grey80", col = "grey40", lwd = 0.3) +
  geom_point(data = points, aes(x, y, color = land)) +
  theme_minimal() +
  theme(axis.text = element_blank())

[jlacko]

It looks as a bug indeed.

Out of curiosity I tried "the other" mapping package = tmap - getting the same result as you. Hard coding of the Mollweide projection did not work, and I ended up with the same map as you.

So I suspect a bug somewhere along the lines of proj.4 / gdal

points <- data.frame(
  land = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,
           FALSE),
  x = c(-66.89695, 143.28066, 51.63533, -76.16409, 177.77862,
        171.00648, 101.41291, -147.60461, 117.14617, 107.37859),
  y = c(76.77662, -47.098278, -30.987178, 84.906883, -70.977776,
        -7.352523, -78.711968, 23.036947, -51.915945, -3.458423)) %>%
  st_as_sf(coords = c("x", "y")) %>%
  st_set_crs(st_crs(world_proj))

library(tmap)

tm_shape(world_proj) + tm_polygons() +
  tm_shape(points, projection = "+proj=moll") + tm_bubbles(col = "land")

[wchase]

Very strange, thanks for the investigative help, I'll reach out to some map gurus and see if any have the answer.

[Maschette]

Hey Will,
Pretty sure it's a projection thing. I reprojected points from lon-lat to crs = 54009 and it worked fine.

library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)

countries <- ne_countries(returnclass = "sf")
world_proj <- st_transform(countries, crs = 54009)

#points, already projected to same crs
points <- data.frame(
  land = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE,
           FALSE),
  x = c(-66.89695, 143.28066, 51.63533, -76.16409, 177.77862,
        171.00648, 101.41291, -147.60461, 117.14617, 107.37859),
  y = c(76.77662, -47.098278, -30.987178, 84.906883, -70.977776,
        -7.352523, -78.711968, 23.036947, -51.915945, -3.458423)
)

#Reproject the data
library(SOmap)
points2<-SOproj(x=points$x,y=points$y, data=points$land, target="+proj=moll +lon_0=0 +x_0=0 +y_0=0+ellps=WGS84 +datum=WGS84 +units=m")
points2<-as.data.frame(points2)
names(points2)<-c( "x","y","land")

#plot
ggplot() +
  geom_sf(data = world_proj, fill = "grey80", col = "grey40", lwd = 0.3) +
  geom_point(data = points2, aes(x, y, color = land)) +
  theme_minimal() +
  theme(axis.text = element_blank())