I see exactly where the issue lies: when you export the ggplot to a TIFF first and then slap on geographic metadata, you're only mapping the panel's data range to the raster—this misses the alignment between the full plot (including north arrow, scale bar, axis labels, and margins) and the geographic bounds. Let's fix this step by step to get a GeoTIFF that matches exactly what plot(gg) shows, with every element intact and proper georeferencing.

Key Problems in Your Current Workflow #

• You're using the panel's x/y range (lat_long$x_range, lat_long$y_range) which only covers the map data area, not the entire image including margins and annotations.
• The TIFF from ggsave has no inherent link between pixel positions and geographic coordinates, so manually setting the extent misaligns plot elements with the geospatial bounds.

Solution: Align Pixel Coordinates to Geographic Bounds Precisely #

Here's a revised workflow that preserves all ggplot elements and correctly georeferences the output:

Step 1: Prepare Your Data & Plot (Keep Your Existing Code) #

First, retain your data prep and ggplot code—we'll only modify the export and georeferencing steps:

# Load packages (we'll use terra instead of raster for better modern support)
library(rgdal)
library(terra)
library(tiff)
library(ggplot2)
library(ggspatial)
library(sf)

# Get data set
all.stands.predict <- read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/prediction__bug_2021-03-18.csv")
all.stands.predict <- all.stands.predict[all.stands.predict[,3]=="VILA PALMA",]

# Convert to sf object
sites <- st_as_sf(all.stands.predict, coords = c("x", "y"), crs = 4326, agr = "constant")

# Create your ggplot map
gg <- ggplot() + 
  geom_sf(data=sites, color="red") + 
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.3, "in"), pad_y = unit(0.5, "in"), 
                         style = north_arrow_fancy_orienteering) +
  annotation_scale(location = "bl", width_hint = 0.55) +
  xlab("Latitude") + ylab("Longitude") + 
  coord_sf() + 
  theme_bw()

Step 2: Export ggplot to TIFF with Fixed Dimensions #

Export the plot to a TIFF with fixed width/height (in inches) and high DPI. This ensures we can calculate the exact pixel-to-geographic ratio later:

# Define output dimensions (adjust these to match your desired plot size)
plot_width <- 8
plot_height <- 6
dpi_val <- 600

# Save the full plot (including all margins/annotations) to TIFF
ggsave(plot = gg, filename = "gg_full.tiff", device = "tiff",
       width = plot_width, height = plot_height, dpi = dpi_val)

Step 3: Calculate Exact Geographic Extent for the Entire Image #

We need to map the entire image's pixel space to the geographic bounds. Here's how:

1. Grab the panel's geographic range (the actual map area)
2. Calculate the proportion of the image taken up by the panel (excluding margins)
3. Extrapolate to get the full image's geographic bounds

# Extract panel parameters from ggplot
panel_params <- ggplot_build(gg)$layout$panel_params[[1]]
panel_x_range <- panel_params$x_range  # Geographic min/max x
panel_y_range <- panel_params$y_range  # Geographic min/max y

# Extract plot layout dimensions (in inches)
plot_layout <- ggplot_build(gg)$layout
panel_width_in <- plot_layout$panel_widths[1]  # Panel width in inches
panel_height_in <- plot_layout$panel_heights[1]  # Panel height in inches

# Calculate pixel per geographic unit (for x and y axes)
x_pixel_per_unit <- (panel_width_in * dpi_val) / diff(panel_x_range)
y_pixel_per_unit <- (panel_height_in * dpi_val) / diff(panel_y_range)

# Calculate total image's geographic extent
total_x_pixels <- plot_width * dpi_val
total_y_pixels <- plot_height * dpi_val

# Calculate extra pixels outside the panel (left/right, top/bottom)
extra_x_left <- (total_x_pixels - (panel_width_in * dpi_val)) / 2
extra_x_right <- extra_x_left
extra_y_top <- (total_y_pixels - (panel_height_in * dpi_val)) / 2
extra_y_bottom <- extra_y_top

# Extrapolate to full image's geographic bounds
full_x_range <- c(
  panel_x_range[1] - (extra_x_left / x_pixel_per_unit),
  panel_x_range[2] + (extra_x_right / x_pixel_per_unit)
)
full_y_range <- c(
  panel_y_range[1] - (extra_y_bottom / y_pixel_per_unit),
  panel_y_range[2] + (extra_y_top / y_pixel_per_unit)
)

Step 4: Load TIFF & Assign Geographic Metadata #

Use the terra package (modern replacement for raster) to load the TIFF and assign the calculated extent and CRS:

# Load the TIFF as a SpatRaster (terra's equivalent of RasterStack)
raster_obj <- rast("gg_full.tiff")

# Assign the full geographic extent
ext(raster_obj) <- c(full_x_range, full_y_range)

# Assign the CRS (WGS84, EPSG:4326)
crs(raster_obj) <- "EPSG:4326"

# Export as GeoTIFF
writeRaster(raster_obj, filename = "my_geotiff_gg.tif", 
            options = "PHOTOMETRIC=RGB", datatype = "INT1U", overwrite = TRUE)

Step 5: Verify the Result #

Load and plot the GeoTIFF to confirm it matches your original ggplot:

# Load the final GeoTIFF
final_geotiff <- rast("my_geotiff_gg.tif")

# Plot it to check alignment and elements
plotRGB(final_geotiff)

Why This Works #

• We explicitly map every pixel of the full plot (including margins, north arrow, scale bar) to the correct geographic coordinates.
• Fixed plot dimensions and DPI ensure consistent pixel-to-geographic ratios.
• Using terra instead of raster gives more reliable georeferencing and modern support.

Bonus Tip: Asymmetric Margins #

If your plot has uneven margins (e.g., more space on the left than right), use ggplot_build(gg)$layout$panel_ranges to get exact panel positions relative to the full image, instead of assuming equal margins.

内容的提问来源于stack exchange，提问作者Leprechault