Редактировать растр на основе атрибута sf полигонов

Question

Я ломал голову над этим уже несколько дней. У меня есть растр, который мне нужно редактировать на основе атрибутов многоугольника в объекте sf. Два объекта имеют одинаковую протяженность.

Растр, заполненный 9-ю

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)

Два треугольника с атрибутами 1 и 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...

Вот то, что я пытался поместить NA, где находится многоугольник NA. raster::mask должен «Создать новый растровый * объект, который имеет те же значения, что и x, за исключением ячеек, которые являются NA (или другим значением маски) в« маске ». Эти ячейки становятся NA (или другим значением обновления)». Ни один из них не выдает ошибку и не изменяет растр.

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

Редактировать: Также не работает:

Редактировать 2: на самом деле, это работает.

При эксперименте мне не удалось изучите растр достаточно внимательно, чтобы увидеть значения NA. Тем не менее, похоже, что это работает противоположно документации. Это сохранение растровых ячеек, находящихся под полигоном NA, и изменение всего остального на NA. Я все еще в замешательстве.

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

Answer 1

Я согласен, что это неожиданный результат, основанный на документации mask.

Похоже, растровые ячейки будут преобразованы в NA, если центр ячейки не покрыт пространственным многоугольником. Итак, одно из ваших возможных решений было близко rm <- mask( x=r, mask=p.sf[p.sf$attr==1, ] ), но не удалось, потому что это логическое индексирование c(TRUE,NA) пыталось вернуть объект, у которого второй объект имел пустую геометрию. Вот правильный способ выполнить индексацию.

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), ] )

Интересно, что объекты полигона преобразуются в растр с помощью rasterize, затем mask также дает желаемый результат через ожидаемый результат

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)

Answer 2

В моих реальных данных это работало, но это заняло 6-8 часов (~ 13 миллионов растровых ячеек и 300 000 полигонов). Так как часть растра, которая нуждается в замене на NA, - это небольшая часть, я подумал, что может быть более эффективно выбрать соответствующие полигоны, чем остальные. Ответ @ ThetaF C показал возможность использования rasterize для этого без mask.

Документация rasterize меня еще больше смущает, чем mask, так что с тоннами проб и ошибка, я наконец нашел это (используя данные примера выше). С моими реальными данными это занимает всего около 40 минут.

# Subset polys to the part that needs editing in raster
p <- subset( p.sf, is.na( attr ) )

# 'field' in this case is the replacement value 
# (can't seem to replace directly with NA)
rr <- rasterize( x=p, y=r, update=TRUE, field=500 )

# Then change it to NA
rr[ rr==500 ] <- NA

Результаты такие же, как и выше.