0

I am using the following code to extract raster values over a SpatialPolygonDataFrame.

result<-extract(Stack, polygons)

>Stack
class       : RasterStack 
dimensions  : 634, 862, 546508, 13  (nrow, ncol, ncell, nlayers)
resolution  : 10, 10  (x, y)
extent      : 492660, 501280, 4378960, 4385300  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : L2A_T29SM//21_B02_10m//etc. 
min values  :             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,           119.0000000,            88.0000000,             1.0000000,   -0.9932203,   -0.9959184,    1.0000000 
max values  :          1.460000e+04,          1.547100e+04,          1.485900e+04,          1.547600e+04,          1.168700e+04,          1.131400e+04,          1.109100e+04,          1.084500e+04,          1.451000e+04,          1.105000e+04, 9.996238e-01, 9.992918e-01, 4.000000e+00 



> polygons
class       : SpatialPolygonsDataFrame 
features    : 13 
extent      : 493428.6, 499603.9, 4379481, 4383781  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : id, ID_1,    LC 
min values  :  1,    1,  Agri 
max values  :  9,    4, Water

Then I want to create a dataframe out of it but I received the following error:

resultdf<-as.data.frame(result)
Error: arguments imply differing number of rows: 225, 320, 286, 351, 68, 79, 41, 63, 70, 171, 72, 158, 497.

From what I check in the forum it is a problem with the data structure since the dataframe would not have the same length; how can I solve it?

nmtoken
  • 13,355
  • 5
  • 38
  • 87
GCGM
  • 1,118
  • 1
  • 16
  • 35
  • Please research your question before posting. There are many post on StackExchange that address the results of extracting polygons in R. I actually, just addressed this last week (https://gis.stackexchange.com/questions/237133/function-sample-code-to-extract-raster-value-per-polygon-in-r/264484#264484). Search is your friend. – Jeffrey Evans Dec 18 '17 at 16:12
  • 1
    Each polygon feature covers a different number of pixels in the raster. How do you want to make a regular data frame from an irregular list? Do you want a data frame with two columns, Polygon-ID and Raster cell number where Polygon-ID is repeated for each pixel under the polygon with that ID? – Spacedman Dec 18 '17 at 16:29
  • I don't think this is a dupe of that, depending on what format the OP wants the extraction converted to... – Spacedman Dec 18 '17 at 16:31
  • @Spacedman, Yes. the raster is a stack so my intention is to create a df with one column per layer in the stack + one extra column for polygon ID_1 repeated for each pixel under the polygon with that id_1 (not ID - according to my polygon attributes). Then one row per pixel inside the polygons. The solution proposed by @RoberH almost work although I adapated it. result<-extract(stack, polygons, df=TRUE) The problem is that it's using polygon$ID instead of polygon&ID_1 – GCGM Dec 19 '17 at 08:43
  • @JeffreyEvans just checked your link. I made a small comment on your answer. – GCGM Dec 19 '17 at 09:16
  • @JeffreyEvans this is not a dupe, please re-open. The user wants a different format than in the alleged "dupe". – Spacedman Dec 19 '17 at 10:04

2 Answers2

2

Thanks to @Spacedman for having this post reopened. This issue has more or less been solved but, there is some behavior in extract, specific to polygons, lines and buffered points, that I would like to explore in order to provide clarity to those reading this in the future. First, lets create some example data with one raster and two polygons. 

library(raster)
library(sp) 

r <- raster(ncol=36, nrow=18)
  r[] <- runif(ncell(r))
poly <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)), 
                    rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))
plot(r)
  plot(poly, add=TRUE)

Now, we we use extract with df=TRUE, the results are a data.frame with an ID column and the associated pixel values. The data is organized so the polygon index (ID) is duplicated (stacked) for each unique pixel value, commonly referred to as a long format. This is efficient since the number of pixels intersecting each polygon are different (in our example polygon 1 has 38 rows and polygon 2 has 25 rows for a total of 63 unique pixel values).

( v <- extract(r, poly, df=TRUE) )

So, one can query each unique polygon ID to pull the associated data. Or, even use a function such as tapply.

mean( v[v$ID == 1,]$layer )
tapply(v[,2], v$ID, FUN=mean)

Now, if we want to organize the data in a "wide" column format where the unequal vectors are accounted for using NA values, it gets more complicated. The call to extract with df=FALSE returns a list with each element as a vector. Using the longest vector in the list, we expand the length of every other vector with NA values. This gives us an equal length so we can then rbind the data into a data.frame. There may be times where this is a preferable format and can be leveraged using function such as apply. 

( v <- extract(r, poly) )
( v <- lapply(v, `length<-`, max(lengths(v))) )
as.data.frame(do.call(rbind, v))

However, if we are operating on a stack, the most efficient method really is to have the data stacked as we observed above with df=TRUE.

r <- stack(r, r^2, r*10) 
extract(r, poly, df=TRUE)
Jeffrey Evans
  • 31,750
  • 2
  • 47
  • 94
1

Perhaps you want something like this, if you want all the cell values for each polygon:

result <- extract(Stack, polygons, na.rm=TRUE, df=TRUE)

Or if you want a summary statistic like mean, you can do:

result <- extract(Stack, polygons, fun=mean, na.rm=TRUE, df=TRUE)
Robert Hijmans
  • 10,683
  • 25
  • 35
  • thanks for your answer. I have tried it with a small change result<-extract(stack, polygons, df=TRUE) since I do not want to use any fun, just extract the pixel values. The only issue is that it uses polygon$ID in the dataframe. I need to extract the pixel values together with polygon$ID_1 – GCGM Dec 19 '17 at 08:49
  • The ID column is the index of the feature in your polygons. You can use that to get any other column by indexing from your polygons, like: result$ID_1 = polygon$ID_1[result$ID] – Spacedman Dec 19 '17 at 10:07
  • What exactly is extract doing when the fun argument is omitted with df=TRUE on polygon data? Each polygon should have an unequal number of values so, if raw pixel values were being returned I would expect number of columns to conform to the maximum number of values (length) in the extracted polygon set and then NA values filling values in line values that do not match the maximum length eg., line 1, c(1,2,3,4,5), line 2 c(1,2,3,NA,NA), line 3 c(1,2,3,4,NA), ect... – Jeffrey Evans Dec 19 '17 at 19:47
  • Oh I got it, the ID column represents each polygon and its associated values and is duplicated for the same as length of the number of unique values in each polygon. So, I can take v[v$ID == 1,] to grab all rows (values) associated with polygon 1, correct? – Jeffrey Evans Dec 19 '17 at 20:28