Как да извлека растерни стойности от полигонални данни, след което да ги обединя в рамка с пространствени данни?

Бих искал да обединя полигонални данни и растерни данни в един кадър с данни за целите на след това да използвам пакета 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