根据 sf 多边形的属性编辑栅格
Edit raster based on attribute of sf polygons
这几天我一直在纠结这个问题。我有一个栅格，我需要根据 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 应该 "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)." 既不给出错误也不改变光栅。
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), ] )
我同意这是基于 mask 文档的意外结果。
如果像元的中心未被空间多边形覆盖，栅格像元似乎将被转换为 NA。因此，您可能的解决方案之一是接近 rm <- mask( x=r, mask=p.sf[p.sf$attr==1, ] )，但失败了，因为此逻辑索引 c(TRUE,NA) 试图 return 第二个特征具有空几何的对象。这是进行索引的正确方法。
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)
在我的真实数据中，这有效，但需要 6-8 个小时（约 1300 万个栅格单元和 300,000 个多边形）。由于需要用 NA 替换的栅格部分是一小部分，我认为 select 那些对应的多边形而不是其余部分可能更有效。 @ThetaFC 的回答揭示了在没有 mask.
的情况下使用 rasterize 的可能性
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
结果同上
这几天我一直在纠结这个问题。我有一个栅格，我需要根据 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 应该 "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)." 既不给出错误也不改变光栅。
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), ] )
我同意这是基于 mask 文档的意外结果。
如果像元的中心未被空间多边形覆盖，栅格像元似乎将被转换为 NA。因此，您可能的解决方案之一是接近 rm <- mask( x=r, mask=p.sf[p.sf$attr==1, ] )，但失败了，因为此逻辑索引 c(TRUE,NA) 试图 return 第二个特征具有空几何的对象。这是进行索引的正确方法。
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)
在我的真实数据中，这有效，但需要 6-8 个小时（约 1300 万个栅格单元和 300,000 个多边形）。由于需要用 NA 替换的栅格部分是一小部分，我认为 select 那些对应的多边形而不是其余部分可能更有效。 @ThetaFC 的回答揭示了在没有 mask.
rasterize 的可能性
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
结果同上