Extraction from Land Cover Raster Data

Question

I am trying to extract grassland values from a historical land use and land cover database created by USGS. I have some issues with the Raster package and getValues option. The tiff file is too large to add with this post, but it is available online.

The data is available under Land-use and Land-cover Backcasting.

This is my code:

 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')

Landcover file in tiff format:

   Landcover1 <- raster ("CONUS_Backcasting_y1938.tif")

USA counties file:

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

These two files are not in the same projection, so projection:

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

Now, I want to extract the land cover data only for the grassland (there are total 17 land classes, and grassland is coded as '11')

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

But when I use getValues to extract the grassland,

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

it shows error:

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

I thought it is for NA, but I could not get how to solve the problem.

DJack DJack · Accepted Answer · 2018-04-08T07:26:02

extract returns a vector or a matrix while getValues needs a raster as input. It's why you have this error.

Therefore, this should work in your case:

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

Having said that, I am not sure about the use of extract in your workflow. I think what you are looking for is to mask your raster. So I suggest you use:

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

EDIT

According to your comment, what you really want is to perform a zonal statistics (the number of pixels labeled as grasslands (11) in each county). Here are some steps on how to do it:

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

Extraction from Land Cover Raster Data

1 Answers