[GIS] Extracting value of raster with coordinates using R

coordinatesextractrraster

i'm currently working on the distribution of different species on an area, for 23 years.

I'll show you a sample of my data, it will be easier to get.

My dataframe, let's call it df.

   x       y
4001758 3138416
3990685 3088576
4002641 3078682
3946723 3153793
3975356 2974350
4001284 3137528
3946723 3153793
3946723 3153793
4000195 3103181
4000168 3103446
3969985 3104761
3905824 3107504

I created a raster from .tif file.

Here is what i want to do :
I want to extract the value for each cell (with the extract() function) from my raster, where the coordinates corresponding to this value, are also in df.
I did it, so i'm still working on it, and for now, i didn't see any errors.

But the difficulty that i have is that i want to also get the coordinates corresponding for each value in my result of the extract.

Here is what i did :

test_extractcropland<-extract(cropland,df_coord)

(cropland is my raster, df_coord is a dataframe, you have a sample, above)
Do you have any idea of how i could do that ?
I looked on this forum, and stackoverflow (and others), but i didn't find something that i could adapt for my problem, or it's more likely that i didn't see how i could adapt it.

Here is some data : https://we.tl/3xe509KSZO
"coord_for_extract.csv" is what i want to filter with, the other are just raster if you want to test.

df_coord is coord_for_extract, r is my raster

Best Answer

Convert df to a sp object and use sp = TRUE as argument.

What you have:

library(raster)

r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)

xy <- cbind(-50, seq(-80, 80, by=20))
extract(r, xy)
## [1] 626 554 482 410 338 266 194 122  50

What you need:

extract(r, SpatialPoints(xy), sp = T)
## class       : SpatialPointsDataFrame 
## features    : 9 
## extent      : -50, -50, -80, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 1
## names       : layer 
## min values  :    50 
## max values  :   626

You'll be able to obtain coordinates from this sp object.

A workaround could be:

xy <- data.frame(xy)
names(xy) <- c('x','y')
cbind(extract(r, xy, df = T),xy)
##   ID layer   x   y
## 1  1   626 -50 -80
## 2  2   554 -50 -60
## 3  3   482 -50 -40
## 4  4   410 -50 -20
## 5  5   338 -50   0
## 6  6   266 -50  20
## 7  7   194 -50  40
## 8  8   122 -50  60
## 9  9    50 -50  80

Get coordinates of raster cells

First, use cellnumbers = T with extract():

result <- extract(r, xy, cellnumbers = T)

coordinates() retrieves coordinates by cell number, so use result to subset:

coordinates(r)[result[,2],]
##         x   y
##  [1,] -45 -85
##  [2,] -45 -65
##  [3,] -45 -45
##  [4,] -45 -25
##  [5,] -45  -5
##  [6,] -45  15
##  [7,] -45  35
##  [8,] -45  55
##  [9,] -45  75

And, with original table:

cbind(result,coordinates(r)[result[,2],])
##       cells layer   x   y
##  [1,]   626   626 -45 -85
##  [2,]   554   554 -45 -65
##  [3,]   482   482 -45 -45
##  [4,]   410   410 -45 -25
##  [5,]   338   338 -45  -5
##  [6,]   266   266 -45  15
##  [7,]   194   194 -45  35
##  [8,]   122   122 -45  55
##  [9,]    50    50 -45  75
Related Question