Actual pixel sizes of unprojected raster maps

[This article was first published on modTools, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

It is well known, though often dismissed, that the areas of spatial units (cells, pixels) based on unprojected coordinates (longitude-latitude degrees, arc-minutes or arc-seconds) are wildly inconsistent across the globe. Towards the poles, as the longitude meridians approach each other, the actual ground width of the pixels sharply decreases. So, for non-tropical regions, the actual area covered by these pixels can be quite far from the nominal “approximately 1 km at the Equator” (or 10 km, etc.) which is often referred in biogeography and macroecology studies.

The function below combines several tools from the ‘terra’ R package to compute the (mean or centroid) pixel area for a given unprojected raster map. By default, it also produces a map illustrating the latitudinal variation in the pixel areas.

pixelArea <- function(r, # SpatRaster
                       type = "mean", # can also be "centroid"
                       unit = "m",  # can also be "km"
                       mask = TRUE,  # to consider only areas of non-NA pixels
                       map = TRUE) {
  # version 1.2 (23 Jan 2024)
  stopifnot(inherits(r, "SpatRaster"),
            type %in% c("mean", "centroid"))
  r_size <- terra::cellSize(r, unit = unit, mask = mask)
  if (map) terra::plot(r_size, main = paste0("Pixel area (", unit, "2)"))
  areas <- terra::values(r_size, mat = FALSE, dataframe = FALSE, na.rm = FALSE)  # na.rm must be FALSE for areas[centr_pix] to be correct
  if (type == "mean") {
    cat(paste0("Mean pixel area (", unit, "2):\n"))
    return(mean(areas, na.rm = TRUE))
  if (type == "centroid") {
    r_pol <- terra::as.polygons(r * 0, aggregate = TRUE)
    centr <- terra::centroids(r_pol)
    if (map) {
      if (!mask) terra::plot(r_pol, lwd = 0.2, add = TRUE)
      terra::plot(centr, pch = 4, col = "blue", add = TRUE)
    centr_pix <- terra::cellFromXY(r, terra::crds(centr))
    cat(paste0("Centroid pixel area (", unit, "2):\n"))

Here’s a usage example for a central European country:

elev_switz <- geodata::elevation_30s("Switzerland", path = tempdir())

# 0.008333333 degree resolution
# i.e. half arc-minute, or 30 arc-seconds (1/60/2)
# "approx. 1 km at the Equator"

pixelArea(elev_switz, unit = "km")
# Mean pixel area (km2):
# [1] 0.5892033

As you can see, pixel area here is quite different from 1 squared km, and Switzerland isn’t even one of the northernmost countries. So, when using unprojected rasters (e.g. WorldClim, CHELSA climate, Bio-ORACLE, MARSPEC and so many others) in your analysis, consider choosing the spatial resolution that most closely matches the resolution that you want in your study region, rather than the resolution at the Equator when this is actually far away. This was done e.g. in this paper (see Supplementary Files – ExtendedMethods.docx).

You can read more details about geographic areas in the help file of the terra::expanse() function. Feedback welcome!

To leave a comment for the author, please follow the link and comment on their blog: modTools. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)