Как мне извлечь растровые значения из полигональных данных, а затем объединить их во фрейм пространственных данных?

Я хотел бы объединить данные полигонов и растровые данные в один фрейм данных для последующего использования пакета randomForests в R.
Это включает в себя сначала извлечение среднего значения растра для каждого полигона.

Пока у меня есть следующее:

#load libraries
library(raster)
library(rgdal)
library(sp)
library(maptools)

#import raster data 
r <- raster("myRasterdata.tif")

#import polygon data 
p <- readShapePoly("myPolydata.shp")

#extract mean raster value for each polygon
ExtractMyData <- extract(r, p, small=TRUE, fun=mean, na.rm=TRUE, df=FALSE,  nl=1, sp=TRUE)
# note I have also tried this with df=TRUE and sp=FALSE

На выходе получается матрица, которую я могу записать в фрейм данных. Но у него нет пространственных координат или исходных идентификаторов полигонов, поэтому я не знаю, как объединить вывод в ту же базу данных. Я думал, что аргумент sp = TRUE сделает это, но, похоже, он не работает.

Обратите внимание, что мне на самом деле придется преобразовать полигоны в точки (используя метод центроида?) Для целей RandomForests, чтобы я мог догадаться, что мне действительно нужно, - это объединить средние значения растра, соединенные с точками, а не с полигонами.

Любые предложения будут ценны. Спасибо!!


r extract raster spatial-data-frame
person user3251223    schedule 11.03.2014    source источник
comment
Похоже, что опция sp = TRUE должна делать то, что вы хотите. Что находится в p @ data? Кроме того, почему у вас nl = 2?   -  person blindjesse    schedule 11.03.2014
comment
ой, извините, в этом примере nl должно быть 1. Я использовал nl = 2, потому что у меня действительно была стопка растров, с которой я работал, но упростил ее до 1, чтобы задать свой вопрос.   -  person user3251223    schedule 11.03.2014
comment
class (p) сообщает мне, что это SpatialPolygonsDataFrame. Он имеет 12 переменных (1 из которых - POLY_ID) и 6938 наблюдений. Итак, 6938 полигонов, которые указывают наличие или отсутствие (1 или 0) каждого из 8 видов.   -  person user3251223    schedule 11.03.2014
comment
Извините, я имел в виду, что такое p @ data после запуска извлечения ... похоже, что там должны быть значения, если sp = T   -  person blindjesse    schedule 11.03.2014
comment
хммм ... Я думаю, мне нужно сделать writePolyShape (p, p-output), но он, похоже, ничего не производит, поэтому мой p в конечном итоге выглядит так же, как и ввод.   -  person user3251223    schedule 12.03.2014
comment
К сожалению, когда я использовал writePolyShape (p, p-output), он действительно создал шейп-файл, и есть добавленный столбец SP_ID, который имеет тот же идентификатор, что и значения растра извлечения, поэтому похоже, что все, что мне нужно сделать сейчас соединение основано на этом. Идеально!   -  person user3251223    schedule 12.03.2014

Ответы (1)


arrow_upward
7
arrow_downward

Это работает:

library(raster)
library(sp)
library(maptools)


#import polygon data 
data(wrld_simpl)
p <- wrld_simpl

#create raster data 
r <- raster(extent(p))
r[] <- seq_len(ncell(r))


## this does it directly, adding columns "names(r)" to "p" 
p <- extract(brick(r, r * 2), p, fun = mean, na.rm = TRUE, sp = TRUE)

Вы также можете сделать это вручную, посмотрите, как извлечение с функцией агрегирования дает один вектор-столбец:

p$ExtractData <- extract(r, p, fun = mean, na.rm = TRUE)

Или вы можете работать с многослойным растром, столбец за столбцом, например:

b <- brick(r, r * 2)
extr <- extract(b, p, fun = mean, na.rm = TRUE)
for (i in seq_len(ncol(extr))) p[[colnames(extr)[i]]] <- extr[,i]
person mdsumner    schedule 11.03.2014
comment
Я заметил, что вы используете это b <- brick(r, r * 2) ... Для чего нужен r * 2? - person Thai; 23.08.2017
comment
Этот код создает двухслойный кирпич, в первом слое есть r, а во втором - r * 2. Только так есть реалистичный фиктивный набор данных. - person mdsumner; 24.08.2017
comment
О, понятно ... Спасибо. - person Thai; 24.08.2017