从土地覆盖栅格数据中提取

Question

我正在尝试从 USGS 创建的历史土地利用和土地覆盖数据库中提取草原值。我对Raster 包和getValues 选项有一些问题。 tiff 文件太大，无法在此帖子中添加，但可在线获取。

数据在Land-use and Land-cover Backcasting 下可用。

这是我的代码：

 install.packages("raster")
 install.packages("rastervis")
 install.packages("RCurl")
 install.packages("R.utils")
 install.packages("rgdal")
 install.packages("sp")
 install.packages("maptools")
 install.packages("tibble")
 install.packages("ggplot2")
 install.packages("gridExtra")

 library(R.utils)
 library(rgdal)
 library(sp)
 library(maptools)
 library(raster)
 library(rasterVis)
 library(RCurl)
 library(R.utils)
 library(rgdal)
 library('rgdal')
 library('raster')
 library("tibble")
 library('ggplot2')

tiff 格式的土地覆盖文件：

   Landcover1 <- raster ("CONUS_Backcasting_y1938.tif")

美国县档案：

   USA_county <- readOGR("UScounties",layer="UScounties")

这两个文件不在同一个投影中，所以投影：

 newprojection <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 
 +towgs84=0,0,0"
 projected_raster_landcover1 <- projectRaster(Landcover1, crs = 
 newprojection)

现在，我想只提取草地的土地覆盖数据（总共有 17 个土地类别，草地编码为 '11'）

 Landcover1_values <- extract(x = projected_raster_landcover1, 
                            y = USA_county) 

但是当我使用getValues提取草原时，

 Landcover1_values_count<- lapply(Landcover1_values, FUN = function(x) {      
 length(which(getValues(x) == 11)) })

显示错误：

**Error in (function (classes, fdef, mtable)  : 
 unable to find an inherited method for function ‘getValues’ for signature ‘"numeric", "missing", "missing"’** 

我以为是NA，但我不知道如何解决这个问题。

Answer 1

extract 返回一个向量或矩阵，而getValues 需要一个栅格作为输入。这就是为什么你有这个错误。

因此，这应该适用于您的情况：

Landcover1_values_count <- sum(Landcover1_values == 11, na.rm = T)

话虽如此，我不确定在您的工作流程中使用extract。我认为您正在寻找的是掩盖您的栅格。所以我建议你使用：

Landcover1_values <- mask(projected_raster_landcover1, USA_county)
Landcover1_values_count <- sum(Landcover1_values[,] == 11, na.rm = T)

编辑

根据您的评论，您真正想要的是执行区域统计（每个县标记为草原 (11) 的像素数）。以下是有关如何执行此操作的一些步骤：

library(raster)
library(plyr)

# Function for efficient zonal stats using data.table, source: https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017475.html
myZonal <- function (x, z, stat, digits = 0, na.rm = TRUE, 
                     ...) { 

  library(data.table)
  fun <- match.fun(stat) 
  vals <- getValues(x) 
  zones <- round(getValues(z), digits = digits) 
  rDT <- data.table(vals, z=zones) 
  setkey(rDT, z) 
  rDT[, lapply(.SD, fun, na.rm = TRUE), by=z] 
} 

# Add an ID field to the shapefile
USA_county@data$ID <- c(1:length(USA_county@data[,1]))

# Crop raster to 'zone' shapefile extent
r <- crop(projected_raster_landcover1, extent(USA_county))

# Reclassify raster in binary raster with 1 for grasslands and 0 for all others values
r[r != 11] <- 0
r[r == 11] <- 1 

# Rasterize shapefile using ID field
zone <- rasterize(USA_county, r, field="ID", dataType = "INT1U") # Change dataType if nrow(USA_county) > 255 to INT2U or INT4U

# Zonal stats
Zstat <- data.frame(myZonal(r, zone, "sum"))
colnames(Zstat) <- c("ID", "Grassland")

# Merge data
USA_county@data <- plyr::join(USA_county@data, Zstat, by="ID")

# Show results
USA_county@data