1

27個のラスターを含むラスタースタックがあります。空間ポリゴンデータフレームに27個の対応するポリゴンがあります。ポリゴン[i]をラスター[i]にオーバーレイし、ラスター[i]から値を抽出して合計し、ポリゴン[i]内のセル数のカウントを取得して、合計値を＃で除算します。細胞の。言い換えると、ラスターは使用率分布またはカーネル使用密度です。ポリゴンがラスターと重なっている領域で多くの使用が発生していることを知りたいです。ポリゴンのサイズを考慮して、ポリゴン内のセルの数で除算したいと思います。

これを行うスクリプトがありますが、これは、データフレーム内の任意の数の空間ポリゴンによってのみ1つのラスターからデータを抽出することを目的として作成されたものです。それは動作し、醜いです、そして私は今それをもっと流線型のものに変換したいと思います。時間がかかるかもしれないので、助けてくれる誰かが私の周りにいたらいいのにと思いますか？

これは私が与えられたコードであり、私が起こっていると思うことの要約です：

msum99Kern07 = SpatialPolygonDataFrame (many polygons)
KERNWolfPIX07m = Raster (this is a single raster, I have 27 rasters I put into a stack

）。

#Extracting value from raster to many polygons 
sRISK_Moose07m<- extract(KERNWolfPIX07m, msum99Kern07,df=FALSE,method='bilinear')

#Calculate THE SUM FOR EACH polygon#
sRISK_Moose07m<-unlist(lapply(sRISK_Moose07m, function(x) if (!is.null(x)) sum(x, na.rm=TRUE) else NA ))
sRISK_Moose07m<-as.data.frame(sRISK_Moose07m)

#Im not sure why these next commands are needed Im only guessing
#data.frame(levels) as there are many polygons creating a dataframe to put the info into
ID_SUM_07<-as.data.frame(levels(as.factor(msum07locs$ID2)))
#ADD ID TO THE risk data frame 
sRISK_Moose07m$ID<-ID_SUM_07[,1] 

#NUMBER OF CELLS WITHIN POLYGON EXTRACT CELLS/ POLYGON
NB_SUM2007m<-cellFromPolygon(KERNWolfPIX07m, msum99Kern07)
NB_SUM07m<-unlist(lapply(NB_SUM2007m, function(x) if (!is.null(x)) length(x) else NA ))

#####CONVERT TO DATA FRAME
NB_SUM07m<-as.data.frame(NB_SUM07m)

###ADD THE NB OF CELLS TO THE RISK_SUM FILE###
sRISK_Moose07m$NB_CELLS<-NB_SUM07m[,1]

###DIVIDING VALUE by NB CELLS##
sRISK_Moose07m$DIVID<-sRISK_Moose07m$sRISK_Moose07m/sRISK_Moose07m$NB_CELLS

これで、27個のポリゴンを含む空間ポリゴンデータフレームと、27個のラスターを含むラスタースタックができました。ラスター[i]とポリゴン[i]を選択し、重なり合う領域のカーネル密度を抽出、合計、計算したいと思います。ポリゴンとラスターが重なっていない可能性があるため、エラーが発生する可能性があります...Rでこれを確認する方法がまったくわかりません。

私が始めた私のスクリプト：

moose99kern = spatial polygon data frame 27 moose
Rastwtrial = stack of 27 rasters having the same unique name as the ID in moose99kern

mkernID=unique(moose99kern$id)

for (i in length(mkernID)){
           r = Rastwtrial[Rastwtrial[[i]]== mkernID[i]] #pick frm Rasterstack the raster that has the same name
            mp = moose99kern[moose99kern$id == mkernID[i]] #pick from spatialpolygondataframe the polygon that has the same name 

            RISK_MooseTrial<- extract(r, mp, df=T, method'bilinear')
            risksum = (RISK_MooseTrial, function(x) if (!is.null(x)) sum(x, na.rm=TRUE) else NA )#sum all the values that were extracted from the raster

ラスタースタックにインデックスを付ける方法がわからないため、スクリプトが機能し始めません。しかし、それでも、一度に1つのラスター/ 1ポリゴンを通過すると、コードで次に何をすべきかわかりません。これがStackOverflowにとって多すぎる場合は、お詫び申し上げます。私は真剣に立ち往生していて、どこを向くべきかわからない。
これは、ポリゴンの2人のテストデータです。

 dput(mtestpoly) 
    new("SpatialPolygonsDataFrame"
        , data = structure(list(id = structure(1:2, .Label = c("F01001_1", "F07002_1"
    ), class = "factor"), area = c(1259.93082578125, 966.364499511719
    )), .Names = c("id", "area"), row.names = c("F01001_1", "F07002_1"
    ), class = "data.frame")
        , polygons = list(<S4 object of class structure("Polygons", package = "sp")>, 
        <S4 object of class structure("Polygons", package = "sp")>)
        , plotOrder = 1:2
        , bbox = structure(c(6619693.77161797, 1480549.31292137, 6625570.48348294, 
    1485861.5586371), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"
    ), c("min", "max")))
        , proj4string = new("CRS"
        , projargs = NA_character_

dput（Rastwtest）

new("RasterStack"
    , filename = ""
    , layers = list(<S4 object of class structure("RasterLayer", package = "raster")>, 
    <S4 object of class structure("RasterLayer", package = "raster")>)
    , title = character(0)
    , extent = new("Extent"
    , xmin = 1452505.6959799
    , xmax = 1515444.7110552
    , ymin = 6575235.1959799
    , ymax = 6646756.8040201
)
    , rotated = FALSE
    , rotation = new(".Rotation"
    , geotrans = numeric(0)
    , transfun = function () 
NULL
)
    , ncols = 176L
    , nrows = 200L
    , crs = new("CRS"
    , projargs = NA_character_
)
    , z = list()
    , layernames = "Do not use the layernames slot (it is obsolete and will be removed)\nUse function 'names'"
)
4

2 に答える 2

6

何かが恋しいかもしれませんが、あなたは問題を複雑にしすぎていると思います。私にとってあなたは持っています：

  1. ラスターのスタック：ラスターのリスト：ss
  2. ssと同じサイズのポリゴンのリスト：polys

extract（ss、polys）から各ペア（layer、poly）に申し込む必要があります

sapply(1:nlayers(ss), function(i) {
     m <- extract(ss[[i]],polys[i], method='bilinear', na.rm= T)[[1]]
     d <- ifelse (!is.null(m) , sum(m)/length(m), NA)
     d
})

再現可能な例を示さないため、ここに2つの長さの例を示します。

## generate some data
library(raster)
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
## In your case you need something like SpatialPolygons(moose99kern)
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                              Polygons(list(Polygon(cds2)), 2)))
r   <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
r1   <- raster(ncol=36, nrow=18)
r1[] <- seq(-1,-2,length.out=ncell(r1))
ss <- stack(r,r1)
## density compute
sapply(1:nlayers(ss), function(i) {
         ## sum of values of the cells of a Raster ss[[i]] covered by the poly polys[i]
         m <- extract(ss[[i]],polys[i], method='bilinear', na.rm= T)[[1]]
         d <- ifelse (!is.null(m) , sum(m)/length(m), NA)

})

[1] 387.815789  -1.494714
于 2013-02-04T12:46:55.073 に答える
3

Rについて質問するときは、自分のデータではなく、常に再現可能な単純な例を使用してください。おそらくあなたがやりたいことがそのような例では機能するが、データでは機能しない場合を除いて、それでも機能する例と表示されるエラーメッセージを表示します。以下の？extractからのように、通常はヘルプファイルの例から始めることができます。

r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
s <- stack(r, r*2)
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                         Polygons(list(Polygon(cds2)), 2)))

v <- extract(s, polys, small=TRUE)

#cellnumbers for each polygon
sapply(v, NROW)

# mean for each polygon
sapply(v, function(x) apply(x, 2, mean, na.rm=T))

一部のポリゴンがラスターの外側にある場合は、sapplyの関数を改良する必要があります（つまり、NULLを返しますが、「small = TRUE」オプションを使用すると、ラスター内の非常に小さいポリゴンの問題を回避できます。また、「メソッド」がないことに注意してください。 "SpatialPolygon*オブジェクトで抽出するときの引数。

各ポリゴンに非常に多くのセルがある場合にメモリの問題を防ぐ場合を除いて、ループを使用しないでください。

于 2013-02-14T18:13:08.710 に答える