Tedward Tedward - 2 months ago 13
R Question

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

I would like to convert a raster layer into a dataframe and get the coordinates as well.

This works just fine, but doesn't give xy values:

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



This throws an error:

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


The error is:


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


I suspect the problem lies somewhere with the Raster Attribute Table, but am unsure how to proceed. I suppose I could convert the factors to numeric and try going from there (
xy=T
works for non factor rasters), but I'd like to find out why adding
xy=T
gives this error. So my question is really "Why does this happen, and how can I get it to work (return xy values with cell values)?"

The raster layer has factors as data (I'm not sure exactly how to call it, here is the str statement):

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"


The behavior on a dataframe with a structure like this:

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


I'd like to be able to get a column of "category", "x", and "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

'raster' expects that the first column of the raster attribute table is an integer variable named 'ID'. How did you get this layer? Perhaps this is a bug that needs to be fixed, or perhaps you made a mistake when you created it.

Either way, here is a work around that should work:

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

Here is a self-contained example

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

although in this example, you can also do:

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