Add graticules to a SpatRaster in terra

Question

I want to map sample locations on a map of Antarctica.

library(terra)
r<-rast("IBCSO.tif")
#IBCSO.tif was downloaded from https://ibcso.org/current_version/
v<-vect(lonlat, crs="+proj=longlat")
p<-project(v, crs(r))
plot(r)
points(p, col="red", pch=20, cex=1)`

gives me the map.

> rast()
class       : SpatRaster 
dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 

> crs(r)
[1] "PROJCRS[\"WGS 84 / IBCSO Polar Stereographic\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"IBCSO Polar Stereographic\",\n        METHOD[\"Polar Stereographic (variant B)\",\n            ID[\"EPSG\",9829]],\n        PARAMETER[\"Latitude of standard parallel\",-65,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8832]],\n        PARAMETER[\"Longitude of origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8833]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",north,\n            MERIDIAN[90,\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (Y)\",north,\n            MERIDIAN[0,\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Hydrography and nautical charting.\"],\n        AREA[\"Southern hemisphere - south of 50°S onshore and offshore, including Antarctica.\"],\n        BBOX[-90,-180,-50,180]],\n    ID[\"EPSG\",9354]]"

I would like to add a layer of graticules with steps of 5 degrees latitude and 10 degrees longitude on top of this map. How can I do that?

Answer 1

I'm still not 100 % sure if I really got your question, but let me give shot, at least to provide a starting point for others knowing better:

terra::plot() has a grid = TRUE attribute but this option only seems to allow you to draw grid lines of the coordinate reference system of the raster to be displayed. This won't help you here, since you do not want a grid of EPGS: 9354, but EPSG: 4326 ("in steps of 5 degrees latitude and 10 degrees longitude").

Maybe there is an option to override this default which I'm not aware of. So I'm going to apply a little trick here and simply sketch graticules as an additional vector object consisting of lines, e.g. the graticules dataset from NE.

library(terra)
library(dplyr)

# This vector dataset shoud be sufficient for demonstration purposes
r <- rnaturalearth::ne_countries(country = "Antarctica", type = "countries") |> 
  vect() |> 
  project("EPSG:9354") 

# Substitution for `lonlat`
p <- spatSample(r, 30)

# Read, split, filter correspondig to your requirements (= only every second meridian), 
# and merge again; Maybe not the most elegant approach.
graticules <- sf::read_sf("ne_10m_graticules_5.shp") 

lat <- graticules |> filter((direction != "E" & direction != "W") | is.na(direction))
lon <- graticules |> filter(direction == "E" | direction == "W") |> filter(degrees %% 10 == 0)

gr_new <- bind_rows(lat, lon) |> vect() |> project("EPSG:9354")

# Inspect result
plot(r)
plot(p, col = "red", add = TRUE)
plot(gr_new, add = TRUE)
text(gr_new, gr_new$display)

As a result, we can tick off the coordinate system part, but not the "human readable" requirement, since all labels are placed on the same spot - and again, I'm not sure if it is possible to enforce a better label placement when working with terra::text().

Alternatively, you can make use of {tidyterra} to be able to use {ggplot2} with SpatRaster and SpatVector objects from {terra}. A minor tweak of scale_x_continuous seems to produce the graticules you're looking for but then again, labeling along the frame is far from perfect here.

library(tidyterra)
library(ggplot2)

ggplot() +
  geom_spatvector(data = r) +
  geom_spatvector(data = p, col = "red") + 
  scale_x_continuous(breaks = seq(-180, 180, by = 10)) 

Answer 2

Here's my take extending the approach provided by falk-env.

You can use sf::st_graticules() to get more controlled graticules. So with that you can control the steps of the longitude and latitude lines independently:

library(terra)
library(dplyr)

# I downsample here the tif for full reproducibilty and speed increase
# url: https://download.pangaea.de/dataset/937574/files/IBCSO_v2_ice-surface_RGB.tif
# 161.9 Mb
# r <- rast("IBCSO_v2_ice-surface_RGB.tif")
# r2 <- spatSample(r, size = 500000, method ="regular", as.raster = TRUE )
# writeRaster(r2, "IBCSO_downsample.tif")


r <- rast("IBCSO_downsample.tif")


# Substitution for `lonlat`
v <- rnaturalearth::ne_countries(country = "Antarctica", type = "countries") |> 
  vect() |> 
  project("EPSG:9354") 

p <- spatSample(v, 30)

# Create graticule object with sf
grat <- sf::st_graticule(lon=seq(-180,180, 10),
                      lat = seq(-90,0, 5),
                      ndiscr = 5000) %>%
  vect() %>%
  project("EPSG:9354") 


# Base plot

plotRGB(r, axes=TRUE, mar=5)
plot(p, col = "red", add = T)
plot(grat, add=TRUE)

If you want the labels to be in longitude/latitude degrees, you may try with ggplot2 + tidyterra. Basically, overlap the grat object and use the scale_x/y_continuous() trick for labelling only:

# Human readable with ggplot2

lims <- as.vector(ext(r))

library(tidyterra)
library(ggplot2)

ggplot() +
  geom_spatraster_rgb(data=r) + 
  geom_spatvector(data=p, color="red") +
  geom_spatvector(data=grat, color=alpha("grey60", 0.5)) +
  coord_sf(xlim = lims[c("xmin", "xmax")],
           ylim = lims[c("ymin", "ymax")]) +
  scale_x_continuous(breaks = seq(-180, 180, by = 10)) +
  scale_y_continuous(breaks = seq(-90, 90, by = 5))