I am trying to create georeferenced pdf versions of maps created in ggplot so they can be used in the field on tablets. Based on this post sf: Create GeoPDF using ggplot2 output object it seems like is should be possible, however the end result is a mirrored image of the original ggplot map. The code below replicates the issue. I had originally thought that this issue arises from my use a basemap which necessitated doing things with CRS 3857, but even replicating in a 4326 (WGS84) projection I still get the mirror image. I have tried changing the order extents are supplied, changing datatype on the export to geotiff function. The disparity does not occur until the export to geotiff. It seems very close to working, I am wondering if my data, and the nc data being in the North America with negative longitude (and x when in 3857) values are causing the issue. Any insight would be appreciated.
library(sf)
library(raster)
library(gdalUtilities)
library(tidyverse)
nc <- st_read(system.file("shape/nc.shp", package="sf"),
quiet=TRUE) %>%
st_transform(4326)
nc_points <- nc %>%
st_centroid()
GGPLOT_MAP <-
# basemaps::basemap_ggplot(
# ext = nc %>%
# st_buffer(5000),
# map_service = "esri" ,
# map_type = "world_topo_map",
# map_res = 1,
# dpi = 300,
# interpolate = T
# )+
ggplot()+
geom_sf(
data = nc,
inherit.aes = F,
aes(colour = NAME),
fill = NA, lwd=1
) +
geom_sf(
data = nc_points,
inherit.aes = F,
aes(colour = as.character(NAME))
) +
geom_sf_text(
data = nc_points,
inherit.aes = F,
aes(label =NAME
),size = 1)+
coord_sf(expand = F )+
theme(legend.position = "none")
GGPLOT_MAP
ggsave(plot=GGPLOT_MAP, "gg.tiff", device = "tiff", dpi = 600)
# Create a StackedRaster object from the saved plot
stackedRaster <- stack("gg.tiff")
# Get the GeoSpatial Components
lat_long <- ggplot_build(GGPLOT_MAP)$layout$panel_params[[1]][c("x_range","y_range")]
# Supply GeoSpatial data to the StackedRaster
extent(stackedRaster) <- c(lat_long$x_range,lat_long$y_range)
projection(stackedRaster) <- CRS("+proj=longlat +datum=WGS84")
# Create the GeoTiff
writeRaster(stackedRaster, "myGeoTiffgg.tif", options="PHOTOMETRIC=RGB", datatype="INT1U",overwrite=TRUE)
#using gdalutilities vs gdalUtil
#Save raster as GeoPDF
gdalUtilities::gdal_translate("myGeoTiffgg.tif","myGeoTiffgg.pdf",of="PDF", ot="Byte",
co="TILED=YES", a_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
After some testing determined that the mirroring of the ggplot map was occurring when re-importing the tiff file created by ggsave. Using either terra::flip() or raster::flip() fixed this issue ex,
stackedRaster <- rast("gg_base.tif") %>%
flip()
However as L Tyrone pointed out the previous approach would give wrong locations on the final pdf due to the calculated extent being applied across the entire plot including axis and legend space. This incorporates his approach to limit the plot to only the spatial extent and correctly georeference the pdf.
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) %>%
st_transform(4326)
nc_points <- nc %>%
st_centroid()
# Define plot ratio
pr <- st_as_sfc(st_bbox(nc)) %>%
st_as_sf(crs = 4326) %>%
get_asp_ratio()
# Generate initial plot (modified for better output quality)
ggplot() +
geom_sf(
data = nc,
aes(colour = NAME),
fill = NA,
lwd = 1
) +
geom_sf(
data = nc_points,
aes(colour = NAME),
size = 4
) +
geom_sf_text(
data = nc_points,
aes(label = NAME),
size = 5)+
coord_sf(expand = F)+
theme_void() +
theme(legend.position = "none")
# Save plot using pr value
ggsave("gg_fix.tif",
bg = "white",
width = 10 * pr,
height = 10,
dpi = 600)
# Load gg1.tif, assign extent and CRS values
stackedRaster <- rast("gg_fix.tif") %>%
flip()
plot(stackedRaster)
ext(stackedRaster) <- ext(nc)
crs(stackedRaster) <- "EPSG:4326"
# Save modified stackedRaster
writeRaster(stackedRaster, "gg2_fix.tif", overwrite = TRUE)
# Save as georeferenced PDF
gdal_translate("gg2_fix.tif",
"gg2_fix.pdf",
of = "PDF",
ot = "Byte",
co = "TILED=YES",
a_srs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
To get to the original intent of spatial data with a basemap requires a little more manipulation as my method of getting basemaps requires espg 3857 projection. This creates a ggplot with basemap and then reprojects the result back into lat longs that seem to be required for a georeferenced pdf.
# with basemap ---------------------------------------------------------
#do everything in 3857 for basemaps
nc <- nc %>%
st_transform(3857)
nc_points <- nc_points %>%
st_transform(3857)
BOUND <- st_bbox(nc ) %>%
st_as_sfc() %>%
st_as_sf(crs = 3857) %>%
st_buffer(5000)
EXTENT <- st_bbox(BOUND)
# Define plot ratio
pr <- BOUND %>%
get_asp_ratio()
# Generate initial plot (modified for better output quality)
GGPLOT_MAP<- basemaps::basemap_ggplot(
ext = EXTENT,
map_service = "esri" ,
map_type = "world_topo_map",
map_res = 6,
dpi = 600,
interpolate = T
)+
geom_sf(
data = nc,
aes(colour = NAME),
fill = NA,
lwd = 1
) +
geom_sf(
data = nc_points,
aes(colour = NAME),
size = 4
) +
geom_sf_text(
data = nc_points,
aes(label = NAME),
size = 5)+
coord_sf(expand = F)+
theme_void() +
theme(legend.position = "none")
# Save plot using pr value
ggsave(plot=GGPLOT_MAP,filename = "gg_base.tif",
bg = "white",
width = 5 * pr,
height = 5,
dpi = 600)
# Load gg1.tif, assign extent and CRS values
stackedRaster <- rast("gg_base.tif") %>%
flip()
#plot(stackedRaster)
ext(stackedRaster) <- ext(BOUND)
crs(stackedRaster) <- "EPSG:3857"
#reproject back to 4326 as needed for pdf
stackedRaster <-terra::project(stackedRaster, crs("EPSG:4326"))
# Save modified stackedRaster
writeRaster(stackedRaster, "gg2_base.tif", overwrite = TRUE)
# Save as georeferenced PDF
gdal_translate("gg2_base.tif",
"gg2_base.pdf",
of = "PDF",
ot = "Byte",
co = "TILED=YES",
a_srs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")