rggplot2

Force square plot area with geom_sf in ggplot

I am trying to plot with ggplot an sf object by expanding the plot area to a square area. This is the replicable code:

library(dplyr)
library(ggplot2)
library(rnaturalearth)

world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot() +
  geom_sf(data = world%>%filter(sovereignt == "Italy"))

The result is this: enter image description here

Now I would like to get the same graph but expanding the area to a square, as visible in the graph produced with the following code:

ggplot() +
  geom_sf(data = world%>%filter(sovereignt == "Italy")) +
  theme(aspect.ratio = 1)

enter image description here

The problem is that in the graph above the plot area is square but the sf object is deformed. Does anyone have any ideas to solve it?

Solution

  • As the ne_countries() coordinates are defined by a geographic coordinate system (GCS), this makes it is difficult to produce a square plot. To make it easier to produce a square plot, I recommend projecting your data to a projected coordinate system (PCS) such as WGS84 Pseudo Mercator/EPSG:3857 or the relevant UTM Zone.

    The workflow:

    1. Project data using sf::st_transform()
    2. Return points of the bounding box of your data
    3. Determine whether longitude or latitude extent is greater, and return half the difference between the two e.g the diff variable in this example
    4. Use diff to increase the extent of the lesser dimension in your plot e.g. longitude in this example

    Note that depending on where your data are in the world, the difference between the original WGS84/EPSG:3857 version and your PCS version may vary. However, there are a number of PCS available that will reduce any difference relative to where your data are.

    library(rnaturalearth)
library(dplyr)
library(ggplot2)
library(patchwork)

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

# Create projected and unprojected sf objects
italy3857 <- filter(world, sovereignt == "Italy") |> st_transform(3857)
italy4326 <- filter(world, sovereignt == "Italy")

# Function to get polygon from boundary box
bbox_polygon <- function(x) {
  bb <- sf::st_bbox(x)
  
  p <- matrix(
    c(bb["xmin"], bb["ymin"], 
      bb["xmin"], bb["ymax"],
      bb["xmax"], bb["ymax"], 
      bb["xmax"], bb["ymin"], 
      bb["xmin"], bb["ymin"]),
    ncol = 2, byrow = T
  )
  
  sf::st_polygon(list(p))
}

# Create points of bbox
sf_p <- st_sfc(bbox_polygon(italy3857)) |>
  st_as_sf(crs = st_crs(italy3857)) |>
  st_cast("POINT")

sf_p
# Simple feature collection with 5 features and 0 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 737796 ymin: 4395685 xmax: 2057834 ymax: 5955490
# Projected CRS: WGS 84 / Pseudo-Mercator
#                         x
# 1  POINT (737796 4395685)
# 2  POINT (737796 5955490)
# 3 POINT (2057834 5955490)
# 4 POINT (2057834 4395685)
# 5  POINT (737796 4395685)

# Get longitude distance
x <- as.integer(st_distance(sf_p[1,], sf_p[4,]))
x
# [1] 1320038

# Get latitude distance
y <- as.integer(st_distance(sf_p[1,], sf_p[2,]))
y
# [1] 1559805

# Subtract max from min, divide by 2
diff <- (y - x) / 2

# Plot unprojected sf (for comparison)
p1 <- ggplot() +
  geom_sf(data = italy4326) +
  labs(title = "WGS84/EPSG:4326")

p2 <- ggplot() +
  geom_sf(data = italy3857) +
  coord_sf(xlim = c(st_bbox(italy3857)[[1]] - diff, 
                    st_bbox(italy3857)[[3]] + diff)) +
  labs(title = "WGS84 Pseudo Mercator/EPSG:3857")

p1 + p2

    1