Извлечение из растровых данных земного покрова

avatar
AfromDhaka
8 апреля 2018 в 01:30
989
1
1

Я пытаюсь извлечь значения пастбищ из исторической базы данных землепользования и растительного покрова, созданной Геологической службой США. У меня есть некоторые проблемы с пакетом Raster и параметром getValues. Файл tiff слишком велик, чтобы добавить его к этому сообщению, но он доступен в Интернете.

Данные доступны в разделе История землепользования и растительного покрова.

Это мой код: 

 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, но я не мог понять, как решить проблему.

Источник

Ответы (1)

avatar
DJack
8 апреля 2018 в 07:26
0

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
AfromDhaka
8 апреля 2018 в 19:59
0

Привет! Большое спасибо за помощь. Это идеально подходит для варианта с маской, и я могу неправильно понять «сумму». Landcover1_values_count дает только одно число, как я могу сделать это по округам? Мне нужны данные о пастбищах округа. Большое спасибо еще раз!

DJack
8 апреля 2018 в 20:56
0

Я сделал правку. Я предполагаю, что USA_county — это шейп-файл с многоугольником для каждого интересующего округа.