How to plot coordinates on a spherical shape file correctly?

Question

So I'm trying to plot a bunch of coordinates on the earth and track how many coordinates are in each country. I have plotted the map and coordinates fine, but when I try to use the intersection to count how many coordinates fall within each country (polygon) it results in error. I've tried using the st_make_valid function to fix the earth shape file but it messes up the geometry. I'm new to using R so any help would be greatly appreciated.

I have used the following code to plot the earth shape file and the coordinates on top:

library(tidyverse)
library(sf)
library(rmapshaper)
library(rnaturalearth)
library(rnaturalearthdata)
library(sp)
library(raster)

###############

# Load Data

###############

# Read in data from .csv file

MeteoriteData <- read.csv("C:/Users/ChaseDickson_/Desktop/College/AERO 689/Semester Project/Meteorite Landings.csv")

# Convert these points to an SF object, specifying the X and Y

#  column names, and supplying the CRS as 4326 (which is WGS84)

MeteoriteData.sf <- st_as_sf(MeteoriteData, coords=c('long', 'lat'), crs=4326)

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

MeteoriteMap <- ggplot(data = world) +
  geom_sf() +
  geom_sf(data = MeteoriteData.sf, size = 0.5, shape = 23, fill = "darkred") +
  theme_bw()

MeteoriteMap

Which gives the following plot

However, when getting the intersection of the code I used this

intersection <- st_intersection(x = world, y = MeteoriteData.sf)

But it gave the error

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
  Loop 96 is not valid: Edge 743 crosses edge 998

To fix this I changed the world sf by adding st_make_valid like so:

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

Now this allows the intersection function to work as such:

intersection <- st_intersection(x = world, y = MeteoriteData.sf)

int_result <- intersection %>%
  group_by(sovereignt) %>%
  count()

And the output is recorded shown below

However, this messes the countries (polygons) up in the plot and will give inaccurate data as the new earth shape file is shown below:

Any help figuring out how to maintain the first plot, but still get the intersection function and count to work after adding st_make_valid would be greatly appreciated!

Spectroscopy

Answer 1

The {rnaturalearth} package has had a long and productive history, but - kind of like the similar {maps} package - it belongs to a different, less demanding era. You should consider doing a Marie Kondo on it: thank it for its service, and let it go.

So instead of trying to repair its failings look for a different instance of world dataset, which is a very common and thus standardized use case.

Consider this piece of code, and note that it is not a single piece of wrong geometry, but 6 (out of 241). Correcting them one by one would be a fruitless task.

library(sf)

rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")|>
  st_is_valid() |>
  table()

# FALSE  TRUE 
#     6   235 

My preferred source of the world country data is the {giscoR} package, which interfaces GISCO spatial dataset, ultimately maintained by Eurostat.

It is very handy, known to be valid and actively maintained.

giscoR::gisco_get_countries(resolution = "20") |>
  st_is_valid() |>
  table()

# TRUE 
#  257 

The rest of your code - the intersection and plotting part - should work just fine once you get rid of invalid geometries.

Change

what a beautiful answer lol! – Vaden

Answer 2

Just a note regarding st_intersection() (and shamelessly building on top of Jindra's answer.. ) -
this happens to be an extremely expensive method to retrieve hit count, especially if resulting geometries are just a by-product that will be disregarded (and it's expensive for this too, check st_join())

For counting you can save a lot (like ~250ms vs ~16s(!) for 100 random points on my machine) by opting for st_intersects():

library(sf)
library(ggplot2)
countries <-  giscoR::gisco_get_countries(resolution = "20")

set.seed(1)
random_points <- data.frame(x = runif(100,-180,180), y = runif(100,-90,90)) |> 
  st_as_sf(coords = c("x","y"), crs = "WGS84")

#> Measure st_intersection():
system.time({
  countries_isect <- st_intersection(countries, random_points)
})
#>    user  system elapsed 
#>   15.97    0.47   16.50

#> Measure st_intersects():
system.time({
  countries$hits <- lengths(st_intersects(countries, random_points))
})
#>    user  system elapsed 
#>    0.22    0.00    0.22

ggplot(countries) +
  geom_sf(data = countries) +
  geom_sf(data = random_points) +
  theme_void()

Results, i.e. non-zero lengths(st_intersects(countries, random_points)):

countries[countries$hits > 0,c("NAME_ENGL","hits")] |> st_drop_geometry()
#>              NAME_ENGL hits
#> 4           Antarctica    7
#> 13           Australia    1
#> 30              Brazil    4
#> 31               China    3
#> 48           Greenland    1
#> 55              Canada    3
#> 116         Kazakhstan    1
#> 117               Laos    1
#> 125           Cambodia    1
#> 141         Mauritania    1
#> 155               Oman    1
#> 162           Paraguay    1
#> 184               Mali    1
#> 185 Russian Federation    5
#> 221      United States    3
#> 223          Venezuela    1
#> 249           Thailand    1

^{Created on 2022-11-23 with reprex v2.0.2}

Cleavable