3

I have a polygon and raster with the same extent:

library(rgdal)
library(raster)
library(rgeos)   
library(sp)

 # make the polygons
df <- structure(list(LAT = c(35.61226, 35.34986, 35.17794, 34.60425,34.40356, 33.94433, 33.41834, 16.89566, 16.89561, 16.89561),
                     LON = c(-9.604802, -9.803048, -9.921968, -10.30782, -10.44971,-10.76656, -11.13053, -24.99776, -24.99788, -24.99773)), 
                .Names = c("LAT","LON"), class = "data.frame", row.names = c(1L, 2L, 3L, 4L, 5L,6L, 7L, 8L, 9L, 10L))

dfSp <- SpatialPointsDataFrame(matrix(c(df$LON, df$LAT), nrow = nrow(df)), data = df)

grid <- makegrid(dfSp, n = 20)
gridSp <- SpatialPointsDataFrame(grid, data = data.frame(id = rownames(grid)))
gridSpRas <- rasterFromXYZ(gridSp)
gridPoly <- rasterToPolygons(gridSpRas, dissolve = T)

# make the raster 
grid2 <- makegrid(dfSp, n = 80)
grid2Sp <- SpatialPointsDataFrame(grid2, data = data.frame(id = rownames(grid2)))
gridSpRas2 <- rasterFromXYZ(grid2Sp)
rvals <- sample(seq(1,8,1),ncell(gridSpRas2),replace=T)
r <- setValues(gridSpRas2, rvals)
extent(r) <- extent(gridPoly)

example I now need to clip the polygons based on values in the raster. I want to remove any polygon that has a cell with a certain value within it. In this case, anytime there is a raster cell with value of 1 within a polygon I want to remove that polygon.

For a simple example such as this I can do:

values <- extract(r, gridPoly)
values2 <- lapply(values, table)

x <- do.call(rbind, values)
y <- cbind(gridPoly@data,x)
collist <- c("1","2","3","4")
sel <- apply(y[,collist],1,function(row) 1 %in% row)
z<-gridPoly[!sel,]

But the real raster has hundreds or thousands of cells per polygon. Is there a simpler way to subset a polygon by a raster? So going from this:

plot(gridPoly)

enter image description here

to this:

plot(z)

enter image description here

  • Can you expand upon "I want to remove any polygon that has a cell with a certain value within it" ? Is any raster cell that touches your polygon considered "within" the polygon, or does the center of the raster cell need to fall within the polygon, or ____ ? – dbaston Jul 10 '18 at 17:26
1

I am not exactly clear as to what you are after here. Are you wanting to speed up the raster::extract function or just optimize code? Here is an approach that collapses your code into a single line. The processing bottleneck will be the extract as, the lapply function is quite fast on large data. If the data is quite large, you will need ample RAM. If you cannot process this in memory then you will have to break up the problem so that it iterates through each polygon (which can be speeded up via multi-threading).

Add packages and create your data

library(raster)
library(sp)

e <- as(raster::extent(-25.7,-6.7,14.3,37.1), "SpatialPolygons") 
  grid <- raster(e, nrow=6, ncol=5)
    grid[] <- 1:ncell(grid)
    grid <- as(grid, "SpatialPolygonsDataFrame")
  grid2 <- raster(e, nrow=10, ncol=10) 
    grid2[] <- sample(1:8, ncell(grid2), replace=TRUE)

This code will create a Boolean based on the condition vector (does not need to be character as the raster values will always be numeric). The TRUE/FALSE Boolean will be added to the SpatialPolygonsDataFrame data.frame (grid object) and can be queried accordingly. This replaces your code block.

grid$condition <- unlist(lapply(extract(grid2, grid), 
    FUN = function(x) { any( c(1,2,3) %in% x ) } )) 

We can factorize the condition column and plot the results.

grid@data$condition <- as.factor(grid@data$condition)   
  spplot(grid, "condition") 
  • So in your example c(1,2,3) finds any cell with a value of 1, 2, or 3 within a polygon and labels the polygon TRUE, correct? – tjr Jul 10 '18 at 18:38
  • Yes, that is correct, the concatenated vector is what is evaluated. However, you could create a vector of the desired values and pass it to the function in lapply by explicit name. – Jeffrey Evans Jul 10 '18 at 19:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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