1

I have a raster (representing a polygon in black) which contains NA values (represented by white pixels) that are located inside the polygon and outside the polygon.

enter image description here

How can I get the index value for the NA values inside the polygon (delimited by the red contour) in the raster ? Next, I would like to use this index to replace the NA values as in the question Change the values of NA cells in a raster by using a geographic subset of the raster as a condition. In this question, the function "cellsFromExtent" is used to get the index value for the NA values in the extent. In my case, I would like to extract NA values within a polygon.

  • 3
    What do you mean by "find" NA values? Do you want [x,y] coordinates, ij index in the raster, cell number, ...? – Jeffrey Evans Apr 4 '16 at 23:37
  • Sorry ! I would like to get the ij index in the raster. Thanks a lot for your help. – Marine Apr 4 '16 at 23:59
  • 1
    please edit your answer to include more details (eg., that you are after the ij index of NA values of a raster within a polygon) and I can provide you an answer. Until you do this, this questions cannot be reopened, and I cannot post the answer for you. – Jeffrey Evans Apr 5 '16 at 0:16
  • @PolyGeo, would you mind reopening this post, the OP has modified the question to provide additional information as to make the question much clearer. – Jeffrey Evans Apr 5 '16 at 1:33
  • @JeffreyEvans no problem - it is now open again – PolyGeo Apr 5 '16 at 1:35
4

You can use the extract function with the cellnumbers = TRUE argument. This will return the cellnumbers and associated values for each polygon.

First, add require libraries and create some example data, raster with NAs and polygons.

library(sp)
library(raster)

set.seed(0)
r <- raster(ncols=10, nrows=10)
  r[] <- runif(ncell(r))
  r[round(runif(25,1,100),digits=0)] <- NA

cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),  
                         Polygons(list(Polygon(cds2)), 2)))
plot(r)
  plot(polys, add=TRUE)

Keep in mind that you do not actually need the row, column index and could just use the cell number(s) directly. This simplifies the code a bit, negating the need for the rowColFromCell function. For the first polygon you would simply do:

( na.cells <- as.data.frame(extract(r, polys[1,], cellnumbers = TRUE)) )
r[na.cells[is.na(na.cells$value),]$cell] <- 0

However, we can expand on things and use the row, column indexes. For all the polygons we create an empty list to store results, iterate through the polygons, use extract to pull cell values and numbers associated with each polygon and then create a row, column index using the rowColFromCell function.

na.row.col <- list()
  for(i in 1:length(polys)) {                                 
    ( na.cells <- as.data.frame(extract(r, polys[i,], cellnumbers = TRUE)) )
    ( na.idx <- na.cells[is.na(na.cells$value),]$cell )
    na.row.col[[i]] <- rowColFromCell(r, na.idx)
  }
na.row.col

If we wanted to create a single data.frame for all of the polygons we could use do.call to collapse our list (commented out) or we can subset just one polygon.

# ( na.row.col <- do.call("rbind", na.row.col) ) 
( na.row.col <- na.row.col[[1]] ) 

Here we use the row,column index to check the values in the raster (should be NA)

r[7,1]
r[5,8]

There is probably a more efficient way than a for loop but, to assign a value to each row, column associated with an NA we can pull the index from the data.frame. This is where using the cell number is more efficient.

for(j in 1:nrow(na.row.col)) { 
  r[na.row.col[j,][1],na.row.col[j,][2]] <- 0
}

If you wanted to get really tricky, you could expand the above for loop to use getValuesFocal and replace the NA value(s) with something like a focal mean.

  • Thank you very much Jeffrey for your answer. I have some problems to convert my raster into polygons. Because of NA values in my raster, I obtain polygons with gaps. I used the function gdal_rasterize {gdalUtils} to convert my raster into polygons. Thanks a lot your help. – Marine Apr 5 '16 at 15:33
  • @Marine, if this answer worked for you please mark as answered so folks know that there was a workable solution. – Jeffrey Evans Apr 5 '16 at 16:03
  • Is there a faster way to return the cell numbers and associated values for the polygon ? I'm testing the code line ( na.cells <- as.data.frame(extract(r, polys[1,], cellnumbers = TRUE)) ) and it's very long ! My raster has 48096864 cells and my polygon is large. Thanks a lot for your help. – Marine Apr 5 '16 at 22:51
  • You could mask to your polygon replace the na cells then add back to the original raster. – Jeffrey Evans Apr 5 '16 at 23:03
  • Thank you very much Jeffrey for your answer. I have NA cells in the whole raster (black polygon in the image). So, the polygon encompasses the raster and thus the 48096864 cells. – Marine Apr 5 '16 at 23:13

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.