将多属性栅格中的栅格属性用于 R 中绘图中的颜色级别

Using a raster attribute from a multi-attribute raster for colour levels in a plot in R

我有一个具有大量属性的 raster 对象,我想在 R 中绘制空间数据并按特定属性对其进行颜色编码。我一直无法弄清楚如何使用特定属性的信息来实现这一点。到目前为止,我已经使用 factorValues() 成功提取了 choice 属性,但我无法确定现在如何将此信息合并到 plot() 函数中。我尝试使用栅格包文档中提到的 ratify()level() 函数,但我不明白简化的在线示例如何适用于具有多个属性的栅格。

任何关于如何实现这一目标的建议将不胜感激。

# read in shapefile
shp = readOGR(".", "grid") 

#convert to raster
r = raster(extent(shp))
res(r) = c(1,0.5) 
ra = rasterize(shp, r)

#crop raster to desired extent
rcrop = crop(ra, extent(-12, 2, 29, 51))

# extract attribute value of interest 
f = factorValues(rcrop, 1:420, layer=1, att=17, append.names=FALSE)
# here there are 420 cells in the raster and I am interested in plotting values of attribute 17 of the raster (this is currently a numeric attribute, not a factor)

#extra code to set attribute as the level to use for plotting colours???
rcrop = ratify(rcrop)
rat = levels(rcrop)[[1]] #this just extras row IDs..not what I want
#…

### plot: I want to plot the grid using 7 colours (I would ideally like to specify the breaks myself)
require(RColorBrewer)
cols = brewer.pal(7,"YlGnBu")
#set breaks 
brks = seq(min(minValue(rcrop)),max(maxValue(rcrop),7))
#plot          
plot(rcrop, breaks=brks, col=cols, axis.arg=arg)

下面的内容很老套(对于大型栅格可能表现不佳),但我不确定是否有办法 link col.regions 指定属性。

rasterVis::levelplot 在标记与因子栅格对应的颜色渐变方面做得很好,虽然它提供了一个 att 参数允许您指定您感兴趣的属性,但这似乎只是修改坡道的标签。栅格单元值控制色带如何映射到栅格,所以在我看来,我们需要自己修改单元值。也许@OscarPerpiñán 会在这里插话来证明我错了:)

我们可以创建一个简单的函数来用我们想要的任何属性替换原始单元格值:

switch_att <- function(r, att) {
  r[] <- levels(r)[[1]][values(r), att]
  r
}

让我们从 Natural Earth 下载并导入一个小示例多边形数据集:

library(rasterVis)
library(rgdal)
require(RColorBrewer)

download.file(file.path('http://www.naturalearthdata.com',
                        'http//www.naturalearthdata.com/download/110m/cultural',
                        'ne_110m_admin_0_countries.zip'),
              f <- tempfile())
unzip(f, exdir=tempdir())
shp <- readOGR(tempdir(), 'ne_110m_admin_0_countries')

rasterize矢量数据:

r <- rasterize(shp, raster(raster(extent(shp), res=c(1, 1))))

并使用 levelplot 创建一些地块:

levelplot(switch_att(r, 'continent'), col.regions=brewer.pal(8, 'Set2')) +
    layer(sp.polygons(shp, lwd=0.5))

levelplot(switch_att(r, 'economy'), par.settings=BuRdTheme) + 
  layer(sp.polygons(shp, lwd=0.5))


编辑

有了 Oscar 对 rasterVis 的更新,上面的 switch_att hack 就不再需要了。

devtools::install_github('oscarperpinan/rastervis')
levelplot(r, att='continent', col.regions=brewer.pal(8, 'Set2')) + 
  layer(sp.polygons(shp, lwd=0.5))

将生成与上面第一个相同的数字。