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")
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")
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.
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?
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.