How to move Crimea from Russia to Ukraine in maps with R

Natural Earth’s de facto on-the-ground policy conflicts with de jure boundaries. Use {sf} and R to relocate parts of country shapes.
r
tidyverse
ggplot
ojs
observable plot
gis
maps
Author
Published

Thursday, February 13, 2025

Doi

The Natural Earth Project

The Natural Earth Project provides high quality public domain geographic data with all sorts of incredible detail, at three resolutions: high (1:10m), medium (1:50m), and low (1:110m). I use their data all the time in my own work and research, and the {rnaturalearth} package makes it really easy to get their data into R for immediate mapping. I mean, look at this!

library(tidyverse)
library(sf)
library(rnaturalearth)

# Set some colors
ukr_blue <- "#0057b7"  # Blue from the Ukrainian flag
ukr_yellow <- "#ffdd00"  # Yellow from the Ukrainian flag
rus_red <- "#d62718"  # Red from the Russian flag

clr_ocean <- "#d9f0ff"
clr_land <- "#facba6"

# CARTOColors Prism (https://carto.com/carto-colors/)
carto_prism = c(
  "#5F4690", "#1D6996", "#38A6A5", "#0F8554", "#73AF48", "#EDAD08", 
  "#E17C05", "#CC503E", "#94346E", "#6F4070", "#994E95", "#666666"
)

ne_countries(scale = 110) |> 
  filter(admin != "Antarctica") |> 
  ggplot() + 
  geom_sf(aes(fill = continent), color = "white", linewidth = 0.1) +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(crs = "+proj=robin") +
  theme_void()

Natural Earth’s de facto policy

Maps are intensely political things. There are dozens of disputes over maps and land and territories (e.g., Palestine, Western Sahara, Northern Cyprus, Taiwan, Kashmir, etc.), and many UN member states don’t recognize other UN member states (Israel isn’t recognized by many Arab states; Pakistan doesn’t recognize Armenia).

The Natural Earth Project’s official policy for disputed territories is to reflect on-the-ground de facto control over each piece of land1:

Natural Earth Vector draws boundaries of sovereign states according to defacto status. We show who actually controls the situation on the ground. For instance, we show China and Taiwan as two separate states. But we show Palestine as part of Israel.

Though they claim that this de facto policy “is rigorous and self consistent”, it gets them in trouble a lot. For instance, there are nearly two dozen issues on GitHub about Crimea, which is illegally occupied by Russia but de jure part of Ukraine. There are huge debates over the ethics of the de facto policy.

Treating the Natural Earth de facto policy as a de facto policy

I’m not weighing in on that policy here! I don’t super like it—it makes it really hard to map Palestine, for instance—but it is what it is. In this post I’m treating the de facto policy as the de facto situation of the data.

Natural Earth de jure points of view

Natural Earth’s solution for disputed territories is to offer different options to reflect country-specific de jure points of view. They offer pre-built high resolution shapefiles for 31 different points of views, so it’s possible to download data that reflect de jure boundaries for a bunch of different countries. Their other shapefiles all have columns like fclass_us, fclass_ua, and so on for doing… something?… with the point of view. I can’t figure out how these columns work beyond localization stuff (i.e. changing place names based on the point of view). The documentation doesn’t say much about how to actually use these different points of view, and pre-built medium and low resolution maps don’t exist yet.

For example, the US doesn’t de jure-ily recognize the Russian occupation of Crimea, so if we download the pre-built high resolution (10m) version of the world from the US point of view, we can see Crimea as part of Ukraine (we have to download this manuallyrnaturalearth::ne_countries() doesn’t support POV files):

world_10_us <- read_sf("ne_10m_admin_0_countries_usa/ne_10m_admin_0_countries_usa.shp")

world_10_us |> 
  filter(ADMIN == "Ukraine") |> 
  ggplot() +
  geom_sf(fill = ukr_blue) +
  theme_void()

Unfortunately, the pre-built point-of-view datasets only exist for the 10m high resolution data. If we want to show medium or low resolution maps, we’re stuck with the de facto version of the map, which means Crimea will be shown as part of Russia. Here’s the low resolution version of Ukraine, with Crimea in Russia:

world <- ne_countries(scale = 110, type = "map_units")

ukraine <- world |> filter(admin == "Ukraine")
russia <- world |> filter(admin == "Russia")

ukraine_bbox <- ukraine |> 
  st_buffer(dist = 100000) |>  # Add 100,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Relocating Crimea manually with R and {sf}

Natural Earth’s recommendation is to “mashup our countries and disputed areas themes to match their particular political outlook”, so we’ll do that here. Though we won’t use any of the point-of-view themes or features because I have no idea how to get those to work.

Instead we’ll manipulate the geometry data directly and move Crimea from the Russia shape to the Ukraine shape by extracting the Crimea POLYGON from Russia and merging it with Ukraine.

The actual geometric shapes for all the countries in world are MULTIPOLYGONs, or collections of POLYGON geometric objects. For instance, Russia is defined as a single MULTIPOLYGON:

russia |> st_geometry()
## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: 41.15 xmax: 180 ymax: 81.25
## Geodetic CRS:  WGS 84
## MULTIPOLYGON (((178.7 71.1, 180 71.52, 180 70.8...

We can split MULTIPOLYGONs into their component POLYGONs with st_cast(). Russia consists of 14 different shapes:

russia_polygons <- russia |> 
  st_geometry() |> 
  st_cast("POLYGON")

russia_polygons
## Geometry set for 14 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: 41.15 xmax: 180 ymax: 81.25
## Geodetic CRS:  WGS 84
## First 5 geometries:
## POLYGON ((178.7 71.1, 180 71.52, 180 70.83, 178...
## POLYGON ((49.1 46.4, 48.65 45.81, 47.68 45.64, ...
## POLYGON ((93.78 81.02, 95.94 81.25, 97.88 80.75...
## POLYGON ((102.8 79.28, 105.4 78.71, 105.1 78.31...
## POLYGON ((138.8 76.14, 141.5 76.09, 145.1 75.56...

The second one is the main Russia landmass:

plot(russia_polygons[2])

The last one is the Crimean peninsula:

plot(russia_polygons[14])

Identifying the Crimea POLYGON from a POINT

The only way I figured out what of these POLYGONs were was to plot them individually until I saw a recognizable shape. And if I use a different map (like the 50m or 10m resolution maps), there’s no guarantee that Russia will have 14 POLYGONs or that the 14th one will be Crimea. We need a more reliable way to find the Crimea shape.

One way to do this is to create a POINT object based somewhere in Crimea and do some geometric set math to identify which Russian POLYGON contains it. The point 45°N 34°E happens to be in the middle of Crimea:

crimea_point <- st_sfc(st_point(c(34, 45)), crs = st_crs(world))

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  geom_sf(data = crimea_point) +
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

We can use it with st_intersects() to identify the Russia POLYGON that contains it:

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons <- world |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon <- russia_polygons |>
  keep(\(x) st_intersects(x, crimea_point, sparse = FALSE))

# This is the same as russia_polygons[14]
plot(crimea_polygon)

Extracting the Crimea POLYGON from Russia

We can then remove that polygon from Russia and recombine everything back into a MULTIPOLYGON. It works!

# Remove Crimea from Russia
new_russia <- russia_polygons |>
  discard(\(x) any(st_equals(x, crimea_polygon, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = new_russia, fill = rus_red) + 
  geom_sf(data = ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Adding the Crimea POLYGON to Ukraine

Next we need to merge crimea_polygon with Ukraine. We’ll convert Ukraine to its component POLYGONs, combine those with Crimea, and recombine everything back to a MULTIPOLYGON. It also works!

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons <- world |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine <- st_union(c(ukraine_polygons, crimea_polygon)) |>
  st_cast("MULTIPOLYGON")

ggplot() +
  geom_sf(data = world, fill = clr_land) +
  geom_sf(data = new_russia, fill = rus_red) + 
  geom_sf(data = new_ukraine, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_bbox["xmin"], ukraine_bbox["xmax"]), 
    ylim = c(ukraine_bbox["ymin"], ukraine_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Updating Russia and Ukraine in the full data

The last step is to modify the full world dataset and replace the existing geometry values for the two countries with the updated boundaries:

world_un <- world |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine,
    admin == "Russia" ~ new_russia,
    .default = geometry
  ))

Now that world_un has the corrected boundaries in it, it works like normal. Here’s a map of Eastern Europe, colored by mapcolor9 (a column that comes with Natural Earth data that lets you use 9 distinct colors to fill all countries without having bordering countries share colors). Crimea is in Ukraine now:

eastern_eu_bbox <- ukraine |> 
  st_buffer(dist = 700000) |>  # Add 700,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world_un, aes(fill = factor(mapcolor9)), linewidth = 0.25, color = "white") +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(
    xlim = c(eastern_eu_bbox["xmin"], eastern_eu_bbox["xmax"]), 
    ylim = c(eastern_eu_bbox["ymin"], eastern_eu_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

The whole game

Everything above was fairly didactic, with illustrations at each intermediate step. Here’s the whole process all in one place:

world_110 <- ne_countries(scale = 110, type = "map_units")

crimea_point_110 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_110))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_110 <- world_110 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_110 <- russia_polygons_110 |>
  keep(\(x) st_intersects(x, crimea_point_110, sparse = FALSE))

# Remove Crimea from Russia
new_russia_110 <- russia_polygons_110 |>
  discard(\(x) any(st_equals(x, crimea_polygon_110, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_110 <- world_110 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")

world_un_110 <- world_110 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_110,
    admin == "Russia" ~ new_russia_110,
    .default = geometry
  ))

Moving Crimea with medium resolution (50m) data

This same approach works for other map resolutions too, like 50m:

world_50 <- ne_countries(scale = 50, type = "map_units")

crimea_point_50 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_50))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_50 <- world_50 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_50 <- russia_polygons_50 |>
  keep(\(x) st_intersects(x, crimea_point_50, sparse = FALSE))

# Remove Crimea from Russia
new_russia_50 <- russia_polygons_50 |>
  discard(\(x) any(st_equals(x, crimea_polygon_50, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_50 <- world_50 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
new_ukraine_50 <- st_union(c(ukraine_polygons_50, crimea_polygon_50)) |>
  st_cast("MULTIPOLYGON")

world_un_50 <- world_50 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_50,
    admin == "Russia" ~ new_russia_50,
    .default = geometry
  ))

Here’s a higher quality map of Eastern Europe with Crimea in Ukraine:

ggplot() +
  geom_sf(data = world_un_50, aes(fill = factor(mapcolor9)), linewidth = 0.25, color = "white") +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(
    xlim = c(eastern_eu_bbox["xmin"], eastern_eu_bbox["xmax"]), 
    ylim = c(eastern_eu_bbox["ymin"], eastern_eu_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

Using the adjusted Natural Earth data as GeoJSON in Observable JS

This updated shapefile works with Observable Plot too (see here for more about how to make nice maps with Observable), but requires one strange tweak because of weird behavior with the GeoJSON file format.

Broken GeoJSON

Let’s export the adjusted geographic data to GeoJSON:

# Save as geojson for Observable Plot
st_write(
  obj = world_un, 
  dsn = "ne_110m_admin_0_countries_un_BROKEN.geojson", 
  driver = "GeoJSON",
  quiet = TRUE,
  delete_dsn = TRUE  # Overwrite the existing .geojson if there is one
)

And then load it with Observable JS:

world_broken = FileAttachment("ne_110m_admin_0_countries_un_BROKEN.geojson").json()

clr_ocean = "#d9f0ff"
clr_land = "#facba6"
ukr_blue = "#0057b7"
ukr_yellow = "#ffdd00"
rus_red = "#d62718"

And then plot it:

Plot.plot({
  projection: "equal-earth",
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_broken, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }) 
  ]
})

lol what even. The new Ukraine shape seems to have broken boundaries that distort everything else in the map. Weirdly, Ukraine is filled with the ocean color while the rest of the globe—both the ocean and whatever countries didn’t have their borders erased—is the color of land.

Let’s zoom in on just Ukraine:

ukraine = world_broken.features.find(d => d.properties.name === "Ukraine")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_broken, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(ukraine, { fill: ukr_blue })
  ]
})

¯\_(ツ)_/¯. Now Ukraine is the color of the ocean and the whole rest of the world is the dark blue of the Ukrainian flag. And it didn’t zoom in at all.

GeoJSON and ↻ ↺ winding order ↻ ↺

This is a symptom of an issue with GeoJSON winding order. GeoJSON cares about the direction that country borders (and all LINESTRING elements) are drawn in. Exterior borders should be drawn counterclockwise; interior borders should be drawn clockwise. If a geographic shape doesn’t follow this winding order, bad things happen. Specifically:

a shape that represents a tiny speck of land becomes inflated to represent the whole globe minus that tiny speck of land, the map fills with a uniform color, the local projection explodes. (via @fil)

That’s exactly what’s happening here. Somehow the winding order is getting reversed when we combine Ukraine with Crimea. {sf} itself doesn’t care about winding order, so everything works fine within R; GeoJSON is picky about winding order, so things break.

Fixing it is tricky though! {sf} uses a bunch of different libraries behind the scenes to do its geographic calculations, including GEOS and S2, and they all have different approaches to polygon creation. Apparently GEOS goes clockwise by default while others go counterclockwise, or something. It should theoretically be possible to fix by adding st_sfc(check_ring_dir = TRUE) after making the new Ukraine shape:

# It would be cool if this worked but it doesn't :(
new_ukraine <- st_union(c(ukraine_polygons, crimea_polygon)) |>
  st_sfc(check_ring_dir = TRUE) |> 
  st_cast("MULTIPOLYGON")

But that doesn’t change anything (nor does it work for this person at GitHub).

Clean GeoJSON with correct winding order

BUT there’s another solution. We can force {sf} to not use the S2 library (which it uses by default, I guess?), since S2 seems go in the wrong direction. If we turn off S2 with sf_use_s2(FALSE), make the new Ukraine shape, and then turn S2 back on with sf_use_s2(TRUE), things work!

# Add Crimea to Ukraine
sf_use_s2(FALSE)
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")
sf_use_s2(TRUE)

Here’s the full process with the 110m map:

world_110 <- ne_countries(scale = 110, type = "map_units")

crimea_point_110 <- st_sfc(st_point(c(34, 45)), crs = st_crs(world_110))

# Extract the Russia MULTIPOLYGON and convert it to polygons
russia_polygons_110 <- world_110 |> 
  filter(admin == "Russia") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Extract the Russia polygon that has Crimea in it
crimea_polygon_110 <- russia_polygons_110 |>
  keep(\(x) st_intersects(x, crimea_point_110, sparse = FALSE))

# Extract the Ukraine MULTIPOLYGON and convert it to polygons
ukraine_polygons_110 <- world_110 |> 
  filter(admin == "Ukraine") |> 
  st_geometry() |> 
  st_cast("POLYGON")

# Add Crimea to Ukraine
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
new_ukraine_110 <- st_union(c(ukraine_polygons_110, crimea_polygon_110)) |>
  st_cast("MULTIPOLYGON")
## although coordinates are longitude/latitude, st_union assumes that they are planar
sf_use_s2(TRUE)
## Spherical geometry (s2) switched on

# Remove Crimea from Russia
new_russia_110 <- russia_polygons_110 |>
  discard(\(x) any(st_equals(x, crimea_polygon_110, sparse = FALSE))) |> 
  st_combine() |> 
  st_cast("MULTIPOLYGON")

# Add the modified Russia and Ukraine to the main data
world_un_110_fixed <- world_110 |>
  mutate(geometry = case_when(
    admin == "Ukraine" ~ new_ukraine_110,
    admin == "Russia" ~ new_russia_110,
    .default = geometry
  ))

# Save as GeoJSON
st_write(
  obj = world_un_110_fixed, 
  dsn = "ne_110m_admin_0_countries_un.geojson", 
  driver = "GeoJSON",
  quiet = TRUE,
  delete_dsn = TRUE
)

Here’s the new world map with the the correct Ukraine:

world_fixed = FileAttachment("ne_110m_admin_0_countries_un.geojson").json()

Plot.plot({
  projection: "equal-earth",
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_fixed, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }) 
  ]
})

We can zoom in on Ukraine too:

ukraine_good = world_fixed.features.find(d => d.properties.name === "Ukraine")
russia = world_fixed.features.find(d => d.properties.name === "Russia")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_good, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_fixed, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia, { fill: rus_red }),
    Plot.geo(ukraine_good, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Alternative data sources

Natural Earth isn’t the only source of geographic data online, and other sources use de jure borders instead of de facto borders, like these:

GISCO

The European Commission’s Eurostat hosts the Geographic Information System of the Commission (GISCO), which provides GIS data for the EU. They offer global shapefiles that follow EU-based de jure borders. The {giscoR} package provides a nice frontend for getting that data into R at 5 different resolutions (1:60m, 1:20m, 1:10m, 1:3m, and super detailed 1:1m!). It does not come with additional metadata for each country, though (i.e. there are no regional divisions, population values, map colors, names in other languages, and so on), so it requires some extra cleaning work if you want those details. For example, can add region information with {countrycode}:

library(giscoR)
library(countrycode)

world_gisco <- gisco_get_countries(
  year = "2024",
  epsg = "4326",
  resolution = "60"
) |> 
  # Add World Bank regions
  mutate(region = countrycode(ISO3_CODE, origin = "iso3c", destination = "region"))

world_gisco |> 
  filter(NAME_ENGL != "Antarctica") |> 
  ggplot() + 
  geom_sf(aes(fill = region), color = "white", linewidth = 0.1) +
  scale_fill_manual(values = carto_prism, guide = "none") +
  coord_sf(crs = "+proj=robin") +
  theme_void()

Since the EU doesn’t de jure-ily recognize the Russian occupation of Crimea, Crimea is in Ukraine:

ukraine_gisco <- world_gisco |> filter(NAME_ENGL == "Ukraine")
russia_gisco <- world_gisco |> filter(NAME_ENGL == "Russian Federation")

ukraine_gisco_bbox <- ukraine_gisco |> 
  st_buffer(dist = 100000) |>  # Add 100,000 meter buffer around the country 
  st_bbox()

ggplot() +
  geom_sf(data = world_gisco, fill = clr_land) +
  geom_sf(data = russia_gisco, fill = rus_red) + 
  geom_sf(data = ukraine_gisco, fill = ukr_blue, color = ukr_yellow, linewidth = 2) + 
  coord_sf(
    xlim = c(ukraine_gisco_bbox["xmin"], ukraine_gisco_bbox["xmax"]), 
    ylim = c(ukraine_gisco_bbox["ymin"], ukraine_gisco_bbox["ymax"])
  ) +
  theme_void() +
  theme(panel.background = element_rect(fill = clr_ocean))

It’s possible to use GISCO data with Observable too. We could load the data into R and clean it up there (like adding regions and other details) and the save it as GeoJSON, like we did with the Natural Earth data.

Or we can grab the original raw GeoJSON from Eurostat directly. For example, here’s the raw 60M 2024 world map using the WGS84 (4326) projection that we grabbed with gisco_get_countries() earlier.

world_gisco = await FileAttachment("https://gisco-services.ec.europa.eu/distribution/v2/countries/geojson/CNTR_RG_60M_2024_4326.geojson").json()

ukraine_gisco = world_gisco.features.find(d => d.properties.NAME_ENGL === "Ukraine")
russia_gisco = world_gisco.features.find(d => d.properties.NAME_ENGL === "Russian Federation")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_gisco, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(world_gisco, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia_gisco, { fill: rus_red }),
    Plot.geo(ukraine_gisco, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Visionscarto

There are JSON-based map files created by @fil at Observable as part of the Visionscarto project.

They’re based on Natural Earth, but with some specific adjustments like adding Crimea to Ukraine, making Gaza a little bit bigger so that it doesn’t get dropped at lower resolutions like 110m, and using UN boundaries for Western Sahara.

Like GISCO, though, these don’t have the additional columns that Natural Earth comes with (country names in a bunch of languages, region and continent designations, map coloring schemes, population and GDP estimates, etc.), and those would need to be added manually in R or Observable or whatever.

import {world110m} from "@visionscarto/geo"

countries110m = topojson.feature(world110m, world110m.objects.countries)
ukraine_visionscarto = countries110m.features.find(d => d.properties.name === "Ukraine")
russia_visionscarto = countries110m.features.find(d => d.properties.name === "Russia")

Plot.plot({
  projection: { 
    type: "equal-earth", 
    domain: ukraine_visionscarto, 
    inset: 50 
  }, 
  width: 800, 
  marks: [
    Plot.sphere({ fill: clr_ocean }),
    Plot.geo(countries110m, {
      stroke: "black",
      strokeWidth: 0.5,
      fill: clr_land
    }),
    Plot.geo(russia_visionscarto, { fill: rus_red }),
    Plot.geo(ukraine_visionscarto, { 
      fill: ukr_blue, 
      stroke: ukr_yellow, 
      strokeWidth: 3
    })
  ]
})

Less automatic sources

There are other sources too, but they require manual downloading:

Citation

BibTeX citation:
@online{heiss2025,
  author = {Heiss, Andrew},
  title = {How to Move {Crimea} from {Russia} to {Ukraine} in Maps with
    {R}},
  date = {2025-02-13},
  url = {https://www.andrewheiss.com/blog/2025/02/13/natural-earth-crimea/},
  doi = {10.59350/28kp0-nbq92},
  langid = {en}
}
For attribution, please cite this work as:
Heiss, Andrew. 2025. “How to Move Crimea from Russia to Ukraine in Maps with R.” February 13, 2025. https://doi.org/10.59350/28kp0-nbq92.