Edit raster based on attribute of sf polygons

Question

I've been puzzling over this for days. I have a raster I need to edit based on polygon attributes in an sf object. The two objects have the same extent.

Raster filled with 9's

r <- raster( nrow=10, ncol=10, xmn=0, xmx=10, ymn=0, ymx=10 )
values(r) <- 9

> r
class      : RasterLayer 
dimensions : 10, 10, 100  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 9, 9  (min, max)

Two triangles with attributes 1 and NA.

# Make two triangles (sfg objects)
p1 <- matrix( c(0,0, 10,0, 0,10, 0,0), ncol=2, byrow=TRUE)
p1 <- st_polygon(list( p1 ) )
p2 <- matrix( c(10,0, 10,10, 0,10, 10,0), ncol=2, byrow=TRUE)
p2 <- st_polygon(list( p2 ) )

#Combine into a simple feature geometry column
p.sfc <- st_as_sfc( list( p1, p2 ))

# Make data frame with 1 attribute
df <- data.frame( attr=c(1,NA))

# Combine df and geometry column into sf object
p.sf <- st_sf( df, p.sfc )
# Set same CRS (has no effect on results)
p.sf <- st_set_crs( x=p.sf, value=4326 )

> p.sf
Simple feature collection with 2 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 10 ymax: 10
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  attr                          p.sfc
1    1 POLYGON ((0 0, 10 0, 0 10, ...
2   NA POLYGON ((10 0, 10 10, 0 10...

Here's what I've tried to put NA where the NA polygon is. raster::mask should "Create a new Raster* object that has the same values as x, except for the cells that are NA (or other maskvalue) in a 'mask'. These cells become NA (or other updatevalue)." Neither gives an error and neither changes the raster.

rm <- mask( x=r, mask=p.sf )
rm <- mask( x=r, mask=as_Spatial( p.sf ) )

Edit: Also not working:

Edit 2: actually, this works.

When experimenting I failed to examine the raster closely enough to see the NA values. However, it seems to work oppositely from the documentation. It is keeping the raster cells that are under NA polygon and changing everything else to NA. I'm still confused.

rm <- mask( x=r, mask=p.sf[ is.na(p.sf$attr), ] )

Yes, as I stated in my answer, it is giving an unexpected result and maintaining the raster cells' values whose center is covered by a spatial polygon feature, regardless of the feature's value. So, you are now selecting the polygon feature that has NA value; that is why it is preserving those raster values. But this is not what you wanted, so given the function's unexpected behavior, my answer gives the desired output. — ThetaFC

ThetaFC ThetaFC · Accepted Answer · 2020-01-10T17:20:39

I agree this is an unexpected result based on mask documentation.

It appears the raster cells will be converted to NA if the center of the cell is not covered by a spatial polygon. So, one of your possible solutions was close rm <- mask( x=r, mask=p.sf[p.sf$attr==1, ] ), but failed, because this logical indexing c(TRUE,NA) was trying to return an object where the second feature had an empty geometry. Here is a correct way to do the indexing.

rm <- mask(x=r, mask=p.sf[!is.na(p.sf$attr),])
rm
class      : RasterLayer 
dimensions : 10, 10, 100  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 9, 9  (min, max)
plot(rm)

#or your idea would have worked with which, because it removes NA values
mask( x=r, mask=p.sf[which(p.sf$attr==1), ] )

Interestingly, if the polygon features are converted to a raster via rasterize, then mask also gives the desired output via an expected result

rNA <- r
values(rNA) <- NA
rs <- rasterize(p.sf, rNA, field='attr', update=TRUE)
plot(rs)

mask_result <- mask(x=r, mask=rs)
plot(mask_result)

Edit raster based on attribute of sf polygons

Raster filled with 9's

Two triangles with attributes 1 and NA.

Edit 2: actually, this works.

2 Answers