rterra

terra::extract() leads to artifacts at map boundary / date line


I want to use terra::extract() to extract values from a pixel raster into a set of polygons (for later calculations). However, this leads to artifacts at the map boundary/date line.

See some example data here:

# libraries needed for this example
library(tidyverse)
library(terra)
library(sf)

# example data: global elevation
(elevation      <- geodata::elevation_global(res = 10, download = TRUE, path = tempdir())) 

# grid made of 200x200 km quadratic polygons
grid_200km     <- terra::rast(extent = ext(-17367530.4451614, 17367530.4451614, -6045743.25594052, 7166680.28573498), 
                               crs = "EPSG:6933",  # metric WGS84 
                               resolution = 200000 # 200x200 km
                   ) %>% as.polygons(crs = "epsg:6933") %>% sf::st_as_sf()               

# plot results
plot(elevation, main="raster data")
plot(grid_200km, main="polygon grid")

original data

After extracting the values, there is a long line at the right-hand side boundary of the map. I don't think these could all be small islands that happen to be just left of the 180° meridian; it seems to be an artifact.

# extract values from raster into polygon grid
testExtract <- terra::extract(elevation, 
                              grid_200km,                        
                              fun=mean,                             
                              na.rm=TRUE,                            
                              bind=TRUE                             
) 


# from spatVector into simple features
testExtract_sf <- st_as_sf(testExtract) %>%
  drop_na(wc2.1_10m_elev)

# plotting
ggplot(testExtract_sf) +
  ggtitle("na.rm = TRUE") +
  geom_sf(aes(fill = wc2.1_10m_elev)) +
  scale_fill_viridis_c() +  # 
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), legend.position="none")

Result after extraction, with na.rm = TRUE

This doesn't happen when I set na.rm = FALSE. In that case, any grid cells with an NA value are set to NA, which is not what I want.

Result after extraction, with na.rm = FALSE

This doesn't seem to be a problem with the elevation dataset. It happened also with a dataset different from this minimal example.

What causes these artifacts and how can they be avoided?


Solution

  • The warning after terra::extract() is relevant: "[extract] transforming vector data to the CRS of the raster". There is probably something weird with the transformation of the grid, eg. the "edge" near 180 degrees go from -180 to 180.

    I tested grid_200km |> st_transform(4326) |> ggplot() +geom_sf(fill ="blue", alpha = 0.5) and its fully blue despite alpha = 0.5, there is probably multilple polygons over each others.

    My suggestion is to reproject the raster to the CRS of your grid, and everything should be good.

    elevation_6933 = terra::project(elevation , crs("EPSG:6933") )
    

    This answer is also relevant.