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"))
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)
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?
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:
sf::st_transform()
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