为什么 as.data.frame 在转换具有因子作为数据的栅格并指定 xy=T 时出错？

Question

我想将栅格图层转换为数据框并获取坐标。

这很好用，但没有给出 xy 值：

as.data.frame(raster_layer)

as.data.frame(x = naip_svm_cropped)

               category
1   forest_broadleafdark
2   forest_broadleafdark
3   forest_broadleafdark
4              grassland
5              grassland
6              grassland

这会引发错误：

as.data.frame(raster_layer, xy = T)

错误是：

Error in match(round(v), rat$ID) : error in evaluating the argument 'x' in selecting a method for function 'match': Error in v[, i] : incorrect number of dimensions

我怀疑问题出在栅格属性 Table 的某个地方，但不确定如何进行。我想我可以将因子转换为数字并尝试从那里开始（xy=T 适用于非因子栅格），但我想找出为什么添加 xy=T 会出现此错误。所以我的问题真的是 "Why does this happen, and how can I get it to work (return xy values with cell values)?"

栅格层有因子作为数据（我不知道具体怎么称呼它，这里是str语句）：

str(naip_svm_cropped@data)
Formal class '.SingleLayerData' [package "raster"] with 13 slots
  ..@ values    : logi(0) 
  ..@ offset    : num 0
  ..@ gain      : num 1
  ..@ inmemory  : logi FALSE
  ..@ fromdisk  : logi TRUE
  ..@ isfactor  : logi TRUE
  ..@ attributes:List of 1
  .. ..$ :'data.frame': 5 obs. of  2 variables:
  .. .. ..$ ID      : num [1:5] 0 1 2 3 4
  .. .. ..$ category: Factor w/ 5 levels "forest_broadleafdark",..: 4 1 2 3 5
  ..@ haveminmax: logi TRUE
  ..@ min       : num 0
  ..@ max       : num 4
  ..@ band      : int 1
  ..@ unit      : chr ""
  ..@ names     : chr "madison_classcombine"

具有如下结构的数据帧的行为：

str(mad_veg_cropped@data)
Formal class '.SingleLayerData' [package "raster"] with 13 slots
 ..@ values    : logi(0) 
 ..@ offset    : num 0
  ..@ gain      : num 1
  ..@ inmemory  : logi FALSE
  ..@ fromdisk  : logi TRUE
  ..@ isfactor  : logi FALSE
  ..@ attributes: list()
  ..@ haveminmax: logi TRUE
  ..@ min       : num -0.719
  ..@ max       : num 1.04
  ..@ band      : int 1
  ..@ unit      : chr ""
  ..@ names     : chr "PercentVeg"

head(as.data.frame(mad_veg_cropped, xy = T))
       x       y   PercentVeg
1 291855.5 4775116  0.7182595
2 291856.5 4775116  0.7402779
3 291857.5 4775116  0.7601378
4 291858.5 4775116  0.7702084
5 291859.5 4775116  0.7774438
6 291860.5 4775116  0.7574666

我希望能够获得 "category"、"x" 和 "y" 的列。

sessionInfo() R version 3.1.2 (2014-10-31) Platform: x86_64-apple-darwin13.4.0 (64-bit) other attached packages: [1] raster_2.3-24 sp_1.0-17

dput(naip_svm_cropped)
    new("RasterLayer"
        , file = new(".RasterFile"
        , name = "/private/var/folders/yj/vjkj1yyx1n510rf_rggqdb640000gr/T/R_raster_tedward/2015-06-04_122215_4840_06850.grd"
        , datanotation = "INT2S"
        , byteorder = structure("little", .Names = "value")
        , nodatavalue = -32768
        , NAchanged = FALSE
        , nbands = 1L
        , bandorder = structure("BIL", .Names = "value")
        , offset = 0L
        , toptobottom = TRUE
        , blockrows = 0L
        , blockcols = 0L
        , driver = "raster"
        , open = FALSE
    )
        , data = new(".SingleLayerData"
        , values = logical(0)
        , offset = 0
        , gain = 1
        , inmemory = FALSE
        , fromdisk = TRUE
        , isfactor = TRUE
        , attributes = list(structure(list(ID = c(0, 1, 2, 3, 4), category = structure(c(4L, 
    1L, 2L, 3L, 5L), .Label = c("forest_broadleafdark", "grassland", 
    "shadow1_tree", "Unclassified", "urban_buildings"), class = "factor")), .Names = c("ID", 
    "category"), row.names = c(NA, -5L), class = "data.frame"))
        , haveminmax = TRUE
        , min = 0
        , max = 4
        , band = 1L
        , unit = ""
        , names = "madison_classcombine"
    )
        , legend = new(".RasterLegend"
        , type = character(0)
        , values = logical(0)
        , color = logical(0)
        , names = logical(0)
        , colortable = c("#000000", "#008B00", "#FF0000", "#FFFF00", "#FF00FF")
    )
        , title = character(0)
        , extent = new("Extent"
        , xmin = 291855
        , xmax = 311023
        , ymin = 4768423
        , ymax = 4775116
    )
        , rotated = FALSE
        , rotation = new(".Rotation"
        , geotrans = numeric(0)
        , transfun = function () 
    NULL
    )
        , ncols = 19168L
        , nrows = 6693L
        , crs = new("CRS"
        , projargs = "+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs"
    )
        , history = list()
        , z = list()
    )

Answer 1

'raster' 期望栅格属性的第一列 table 是一个名为 'ID' 的整型变量。你是怎么得到这一层的？也许这是一个需要修复的错误，或者也许是您在创建它时犯了错误。

无论哪种方式，这里有一个应该可行的解决方法：

xy <- xyFromCell(raster_layer, 1:ncell(raster_layer))
v <- as.data.frame(raster_layer) 
xyv <- data.frame(xy, v)

这是一个独立的例子

library(raster)
r <- raster(nrow=10, ncol=10)
r[] = 1
r[51:100] = 2
r[3:6, 1:5] = 3
r <- ratify(r)

rat <- levels(r)[[1]]
rat$landcover <- c('Pine', 'Oak', 'Meadow')
rat$code <- c(12,25,30)
levels(r) <- rat

xy <- xyFromCell(r, 1:ncell(r))
v <- as.data.frame(r) 
xyv <- data.frame(xy, v)

head(xyv)
#     x  y landcover code
#1 -162 81      Pine   12
#2 -126 81      Pine   12
#3  -90 81      Pine   12

尽管在此示例中，您还可以这样做：

vv <- as.data.frame(r, xy=TRUE)

为什么 as.data.frame 在转换具有因子作为数据的栅格并指定 xy=T 时出错？

Why does as.data.frame give error when converting a raster with factors as data and specifying xy=T?

r

r-raster