Plot SpatialPoints (LAT/LON) on MULTIPOLYGON sf data.frame using R

Question

I want to plot points on a multipolygone map (shps from naturalearthdata.com). Ideally using LAT-LON coordinates. However my points are not displayed, end up somewhere or slip out of the plot after zooming into the map (modifying the extent(), respectively)

I tried rasterizing the world_wgs84 map and then using the points() function to plot the coordinates. This works but gives me no option to plot borders, change country color fillings later. So I am looking for a way to plot both objects onto each other while keeping the full information of the MULTIPOLYGONE data.frame.

require(sf)

#set path
path_<-"R/Shapefile Natural Earth Data/ne_50m_admin_0_countries/"
shp_<-"ne_50m_admin_0_countries.shp"

#import
world_wgs84 <- st_read(paste0(path_,shp_))

plot(world_wgs84[,8])


# point coordinates
city_<-data.frame("LON"=0, "LAT"=0)
coordinates(city_)<-~LON+LAT
crs(city_)<-crs(world_wgs84)


# plotting
plot(world_wgs84[,4], col=1)
points(city_, col="red", pch=20, cex=2)

This gives me just the map, no points. (Point should appear around the equator, at the same latitude as Greenwich)

This is what my world_wgs84 sf object looks like:

> world_wgs84
Simple feature collection with 242 features and 161 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS:  WGS 84

and these are my SpatialPoints:

> city_
class       : SpatialPoints 
features    : 1 
extent      : 10, 10, 53, 53  (xmin, xmax, ymin, ymax)
crs         : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs

any ideas how to plot these to object onto each other?

Answer 1

This is how to plot points using lon/lat coordinates on top of natural earth shapes:

library(rnaturalearth)
library(tidyverse)

points <- tibble(
  x = c(-100, -90), # lon
  y = c(0, 0), # lat
  value = c("A", "B")
)

world <- ne_countries(scale = "medium", returnclass = "sf")

world %>%
  ggplot() +
  geom_sf(aes(fill = continent)) +
  geom_point(data = points, mapping = aes(x, y, color = value))

^{Created on 2022-02-18 by the reprex package (v2.0.0)}

Answer 2

This gives me just the map, no points. (Point should appear around the equator, at the same latitude as Greenwich)

They are, however for some reason plot doesn't care about coordinates. The easiest solution with plot() would be to create an empty plot with defined area and then add other object to this plot. I do prefer plot() over ggplot() because of simplicity.

library(sf)
library(sp)

world_wgs84 <- st_read("/home/sapi/projekty/test/ne/ne_50m_admin_0_countries.shp")

Let's create an empty plot:

plot(x = 0, y = 0, cex = 0, 
     xlim = c(-180, 180), ylim = c(-90,90), 
     axes = TRUE, xlab = "", ylab = "")

And add our objects with add = TRUE:

plot(world_wgs84[,8], add = TRUE)

city_<-data.frame("LON"=0, "LAT"=0)
coordinates(city_) <- c("LON", "LAT")
plot(city_, add = TRUE, pch = 20, cex = 2, col = "red")

^{Created on 2022-02-20 by the reprex package (v2.0.1)}