Search code examples
rterra

Using xy=TRUE in R terra extract


I am trying to extract the coordinates of the centroid of pixels for which I have points that fall whithin that pixel. Here is a reproducible example

library(terra)
filename <- system.file("ex/elev.tif", package="terra")
r <- rast(filename)

#Create random points in raster
points<-spatSample(r, 10, as.points=TRUE)
crds<-crds(points)
#Extract the coordinates of the random points
crds
            x        y
 [1,] 5.904167 49.97083
 [2,] 6.262500 50.12917
 [3,] 6.512500 49.63750
 [4,] 6.462500 50.03750
 [5,] 6.137500 49.90417
 [6,] 6.520833 49.96250
 [7,] 5.904167 49.52083
 [8,] 6.187500 49.68750
 [9,] 6.212500 49.68750
[10,] 5.962500 49.86250

#Extract pixel values and centroids of pixels
values<-extract(r, points, xy=TRUE)
values

ID  elevation          x         y
1 373.000000   6.462500  49.52083
2   5.904167  50.037500 367.00000
3  49.970833 368.000000   6.18750
4        NaN   6.137500  49.68750
5   6.262500  49.904167 392.00000
6  50.129167        NaN   6.21250
7        NaN   6.520833  49.68750
8   6.512500  49.962500 431.00000
9  49.637500 305.000000   5.96250
10        NaN   5.904167  49.86250

I would expect the x y coordinates extracted to be only be slightly off the original coordinates, but instead the values seems to be all mixed up. Am I making a mistake in my code or is there an alternative way of getting these values?


Solution

  • This works fine with the current CRAN version (1.3-4):

    library(terra)
    terra version 1.3.4
    filename <- system.file("ex/elev.tif", package="terra")
    r <- rast(filename)
     
    set.seed(7222021)
    points <- spatSample(r, 5, as.points=TRUE)
    crds(points)
    #            x        y
    #[1,] 5.929167 49.67083
    #[2,] 6.279167 49.94583
    #[3,] 6.487500 49.84583
    #[4,] 6.004167 49.49583
    #[5,] 6.379167 49.62083
     
    extract(r, points, xy=TRUE)
    #  ID elevation        x        y
    #1  1       323 5.929167 49.67083
    #2  2       NaN 6.279167 49.94583
    #3  3       NaN 6.487500 49.84583
    #4  4       358 6.004167 49.49583
    #5  5       283 6.379167 49.62083
    

    To update terra, you could run update.packages(). Or if you only want to update terra you should first update Rcpp.

    install.packages(c('Rcpp', 'terra'))