R：提取到 shapefile 中的栅格数据作为矩阵添加 - writeOGR 中的错误
R: Raster data extracted to shapefile added as a matrix - error in writeOGR
我使用以下方法将多个栅格图层的中值提取到多边形 shapefile 中：
#read shapefile
huc12 <- readOGR(dsn=gdb1, layer="HUC12_proj")
#read raster
temp.avg <- raster("projected_climate_rasters/temp_avg_copy")
#extract median of raster for each polygon
huc12$Temp.Avg.Med <- extract(temp.avg, huc12, fun = median)
我将数据转换为数字：
huc12@data <- transform(huc12@data, Temp.Avg.Med = as.numeric(Temp.Avg.Med))
数据似乎都已附加到 shapefile；但是，其中一组附加数据似乎已作为矩阵附加。例如，
Temp.Avg.Med 数据在 huc12@data 下显示为：
.. ..$ Temp.Avg. Med: num [1:877000] 962 ...
而另一个 (PET_AnCV) 显示为：
.. ..$ PET_AnCV : num [1:87700, 1] 94.4 ...
直到我尝试编写新的 shapefile 时才发现这一点，它给出了以下错误：
> write_shape(huc12, "huc12")
Error in writeOGR(shp, dir, base, driver = "ESRI Shapefile", overwrite_layer = TRUE) :
Can't convert columns of class: matrix; column names: PET_AnCV
有什么建议吗？ （考虑到所涉及的数据，我无法使其可重现，并且无法使用示例 shapefile 进行重现。）
这是因为从 PET_AnCV 中提取的值是矩阵而不是向量。当您从多层栅格中提取值时，这是正常的，您会得到一个矩阵，每个要素一行（ie. 多边形），每层一列。在这种情况下，您只有一层（和一列），因此在将其分配给 shapefile 时，您可以轻松地 select 相应的列：
vals <- extract(PET_AnCV,huc12)
huc12$PET_AnCV <- vals[,1]
但这并不是最好的方法，我将进一步解释并举个例子：
您收到的错误消息告诉您整个矩阵已分配给 shapefile 属性的一个列 table。这在 R 中很好（因为它是列表中的一个槽），但是一旦您尝试将文件写入磁盘，您 运行 就会遇到问题。
您可以使用 foreign 包的功能将结果存储在 shapefile 的 dbf 文件中，而不是尝试将提取的值直接写入 shapefile。这是一种更简洁的方法，还可以让您更好地控制自己正在做的事情。如果你的 shapefile 很大 and/or 有很多多边形，这也会快得多。
至于一个可重现的例子，怎么样：
我正在使用奥地利的 Admin-1 边界（有 11 个不同的多边形）作为您的 shapefile 的附件：
# load libraries
library(raster)
library(rgdal)
library(foreign)
# download shapefile which is used for extraction
huc12 <- getData('GADM', country='AUT', level=1)
这里我需要将它写入磁盘，因为我还没有 dbf 文件。既然你已经有了，你可以跳过这部分。
# write to disk
writeOGR(huc12,dsn = 'huc12.shp','huc12',driver = 'ESRI Shapefile')
创建一个 10 层光栅砖，其中包含一些要从中提取的随机值：
# create fake value rasterbrick
ras <- raster()
PET_AnCV <- lapply(1:10,function(x) setValues(ras,runif(ncell(ras))))
PET_AnCV <- do.call(brick,PET_AnCV)
现在我正在提取值并使用 cbind 将它们附加到 dbf 文件。变量名称将是与多层栅格的图层名称相对应的列名称。如果这些是通用的或描述性不够，您可以使用 colnames(vals) <- new_name_vector.
轻松更改它们
# extract values and append to dbf
vals <- extract(PET_AnCV,huc12,fun=median)
dbf <- read.dbf('huc12.dbf')
dbf <- cbind(dbf,vals)
write.dbf(dbf,'huc12.dbf')
HTH
我使用以下方法将多个栅格图层的中值提取到多边形 shapefile 中：
#read shapefile
huc12 <- readOGR(dsn=gdb1, layer="HUC12_proj")
#read raster
temp.avg <- raster("projected_climate_rasters/temp_avg_copy")
#extract median of raster for each polygon
huc12$Temp.Avg.Med <- extract(temp.avg, huc12, fun = median)
我将数据转换为数字：
huc12@data <- transform(huc12@data, Temp.Avg.Med = as.numeric(Temp.Avg.Med))
数据似乎都已附加到 shapefile；但是，其中一组附加数据似乎已作为矩阵附加。例如， Temp.Avg.Med 数据在 huc12@data 下显示为：
.. ..$ Temp.Avg. Med: num [1:877000] 962 ...
而另一个 (PET_AnCV) 显示为：
.. ..$ PET_AnCV : num [1:87700, 1] 94.4 ...
直到我尝试编写新的 shapefile 时才发现这一点，它给出了以下错误：
> write_shape(huc12, "huc12")
Error in writeOGR(shp, dir, base, driver = "ESRI Shapefile", overwrite_layer = TRUE) :
Can't convert columns of class: matrix; column names: PET_AnCV
有什么建议吗？ （考虑到所涉及的数据，我无法使其可重现，并且无法使用示例 shapefile 进行重现。）
这是因为从 PET_AnCV 中提取的值是矩阵而不是向量。当您从多层栅格中提取值时，这是正常的，您会得到一个矩阵，每个要素一行（ie. 多边形），每层一列。在这种情况下，您只有一层（和一列），因此在将其分配给 shapefile 时，您可以轻松地 select 相应的列：
vals <- extract(PET_AnCV,huc12)
huc12$PET_AnCV <- vals[,1]
但这并不是最好的方法，我将进一步解释并举个例子：
您收到的错误消息告诉您整个矩阵已分配给 shapefile 属性的一个列 table。这在 R 中很好（因为它是列表中的一个槽），但是一旦您尝试将文件写入磁盘，您 运行 就会遇到问题。
您可以使用 foreign 包的功能将结果存储在 shapefile 的 dbf 文件中，而不是尝试将提取的值直接写入 shapefile。这是一种更简洁的方法，还可以让您更好地控制自己正在做的事情。如果你的 shapefile 很大 and/or 有很多多边形，这也会快得多。
至于一个可重现的例子，怎么样：
我正在使用奥地利的 Admin-1 边界（有 11 个不同的多边形）作为您的 shapefile 的附件：
# load libraries
library(raster)
library(rgdal)
library(foreign)
# download shapefile which is used for extraction
huc12 <- getData('GADM', country='AUT', level=1)
这里我需要将它写入磁盘，因为我还没有 dbf 文件。既然你已经有了，你可以跳过这部分。
# write to disk
writeOGR(huc12,dsn = 'huc12.shp','huc12',driver = 'ESRI Shapefile')
创建一个 10 层光栅砖，其中包含一些要从中提取的随机值：
# create fake value rasterbrick
ras <- raster()
PET_AnCV <- lapply(1:10,function(x) setValues(ras,runif(ncell(ras))))
PET_AnCV <- do.call(brick,PET_AnCV)
现在我正在提取值并使用 cbind 将它们附加到 dbf 文件。变量名称将是与多层栅格的图层名称相对应的列名称。如果这些是通用的或描述性不够，您可以使用 colnames(vals) <- new_name_vector.
# extract values and append to dbf
vals <- extract(PET_AnCV,huc12,fun=median)
dbf <- read.dbf('huc12.dbf')
dbf <- cbind(dbf,vals)
write.dbf(dbf,'huc12.dbf')
HTH