6

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

3
  • 2
    your extract comes out ordered, just bind it to the df, something like df$raster_value <- extract(cropland, df_coord)
    – Elio Diaz
    Commented Apr 11, 2018 at 17:08
  • The problem that i get is that it just create two times the same column by doing that Commented Apr 11, 2018 at 18:08
  • try the cell indexing cells <- cellFromXY(r, as.matrix(df_coord)); data.frame(xyFromCell(r, cells), value = extract(r, cells)) - it's neat because you separate the geometric task from the extraction (so it can be applied repetitively)
    – mdsumner
    Commented Apr 20, 2018 at 23:02

1 Answer 1

12

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
11
  • Damn it, that's really good, working perfectly ! I had a little question. The problem is : i have 314181 rows with my extract, but when i transform my raster into point, i get less than 300 000 rows. Do you know where it can come from ? (I know without the data it's not easy to tell, but i don't think i did a mistake in my code) Commented Apr 11, 2018 at 18:44
  • Which output gives you dim(unique(df_coord))?
    – aldo_tapia
    Commented Apr 11, 2018 at 18:54
  • 265409 2 for my df_coord, this the one i want to "filter" the raster with, and 314181 4 for the result that i get, it should be max 265 409 Commented Apr 11, 2018 at 19:26
  • Wait, why you should convert raster into points? Points represent only one location per feature, raster represents a specific area per feature
    – aldo_tapia
    Commented Apr 11, 2018 at 19:29
  • Raster have coordinates, x, y and a value for each cell no ? The idea is to extract the coordinates of the raster that only appears in the df_coord Commented Apr 11, 2018 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.