2017-08-22 7 views
2

私はあなたがこの問題を解決するのに役立つことを願っています。申し訳ありませんが、私はこの記事を書いている間にいくつかの間違いを犯した場合、私の英語は少し錆びています。LinearKのR関数を使ったタコの捕獲の分析R

ここに質問があります。私はRで分析したい.shpデータを持っています。.shpは、タコを捕まえるために設定したトラップ行を表す行か、それらの行の直上にあるポイントで、どこで成熟したかを表します。

Clarification image

私は答えることをしようとしている質問です:は、統計的にグループ化されたかどうかタコか?

調査の結果、RaptorとそのlinearK関数を使用して、Maptools、SpatStat、Spというライブラリを使用してその質問に答える必要があると思われます。

t1<- as.linnet(readShapeSpatial("./20170518/t1.shp")) 

は、私は次の警告を取得し、それをトラックとlinnetがオブジェクトを作成するライブラリ

library(spatstat) 
library(maptools) 
library(sp) 

をロード

:ここ

は私がRStudioで使用しているコードです。動作するように思われる

Warning messages: 
1: use rgdal::readOGR or sf::st_read 
2: use rgdal::readOGR or sf::st_read 

すべてが、私はここに同じ警告を受けるが、ときに私の実際の問題は、開始ポイント

p1<- as.ppp(readShapeSpatial("./20170518/p1.shp")) 

とのPPPオブジェクトの作成

plot(t1) 

Plot of t1 object

OKであることを確認するためにそれをプロットそれをプロットしようとしてください:

> plot(p1) 
Error in if (!is.vector(xrange) || length(xrange) != 2 || xrange[2L] < : 
    missing value where TRUE/FALSE needed 
In addition: Warning messages: 
1: Interpretation of arguments maxsize and markscale has changed (in spatstat version 1.37-0 and later). Size of a circle is now measured by its diameter. 
2: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 
3: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 
4: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 
5: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 
6: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 
7: In plot.ppp(x, ..., multiplot = FALSE, do.plot = FALSE) : 
    All mark values are NA; plotting locations only. 

今すぐ何残っていることは、LPPオブジェクト内のオブジェクトに参加するとlinearK機能

> pt1 <- lpp(p1,t1) 
> linearK(pt1) 
Function value object (class ‘fv’) 
for the function r -> K[L](r) 
...................................... 
    Math.label  Description   
r r    distance argument r 
est {hat(K)[L]}(r) estimated K[L](r) 
...................................... 
Default plot formula: .~r 
where “.” stands for ‘est’ 
Recommended range of argument r: [0, 815.64] 
Available range of argument r: [0, 815.64] 

これが今の私の状況であると、それを分析することです。私が知らないのはなぜプロット関数がpppオブジェクトで動作しないのか、そしてlinearK関数のリターンを控える方法です。ヘルプ(linearK)はヒントを提供しませんでした。私は多くのトラックを持っているので、それぞれのポイントがあり、私の希望する結果は、分析されたxトラックのような要約であり、グループ化され、分散され、未知のものです。

ありがとうございました。私がこの問題を解決するのを手伝っていただければ幸いです。

編集:ここには、1日分のshpファイル、トラックとポイント、およびコード付きのtxtファイルを含むzipファイルへのリンクがあります。 https://drive.google.com/open?id=0B0uvwT-2l4A5ODJpOTdCekIxWUU

+0

あなたははるかに例が完全に再現可能にするあなたの良い答えを得るために可能性が高いです。 – Axeman

+0

あなたは、大丈夫です。私は私の質問の一番下にリンクを追加しました。 ありがとうございます。 –

答えて

1

最初の2つの一般的なアドバイス:(1)複雑なオブジェクトを作成するたびに、それを端末で印刷して、期待どおりかどうかを確認します。 (2)エラーが発生したらすぐにtraceback()と入力して出力をコピーしてください。これにより、エラーが検出された場所が正確にわかります。

pppオブジェクトには、スタディ領域(ウィンドウ)の指定が含まれている必要があります。あなたのコードでは、オブジェクトp1は、as.ppp.SpatialPointsDataFrameによって変換された研究領域の指定を含まないクラスSpatialPointsDataFrameのデータを、pppのオブジェクトに変換して作成します。このオブジェクトでは、バウンディングボックスを取ることによってウィンドウが推測されます座標の残念ながら、あなたの例ではのデータポイントがp1にしかないので、デフォルトの境界ボックスは幅0と高さ0の矩形です。[これは印刷されたp1です。]このようなオブジェクトは、通常、 spatstatですが、この特定のオブジェクトは、ウィンドウのサイズがゼロでないことを期待する関数plot.solistのバグを引き起こします。私はあなたのケースでは

をバグを修正、しかし...う、私はあなたがすぐにp1を作成した後

Window(p1) <- Window(t1) 

を行うことをお勧めします。これにより、おそらく意図したウィンドウがp1にあることが保証されます。他のすべてが失敗した場合、シェープファイルにspatstatビネットを読ん

...

+0

あなたの答えは本当にありがとうございました。私はまた、pppオブジェクトがすべてのトラブルを始める場所であることに気づいたので、私は手作業でそれを作成しようと決心しました。それはうまくいくようです。 –

0

私は解決策を見つけるために管理しています。 Adrian Baddeleyが気づいたように、owinオブジェクトに問題がありました。私のポイントのセットを変換するのではなく、手動でpppオブジェクトを作成すると、その問題は回避されているようです(実際は解決されていないようです)。

の機能をrgdal::readOGRに変更しました。これは、最初の機能が廃止されたことが原因であり、これが私の警告の理由です。

これは私が今使っているRスクリプトで、明確にするために次のようにコメントし

#first install spatstat, maptools y sp 

#load them 

library(spatstat) 
library(maptools) 
library(sp) 


#create an array of folders, will add more when everything works fine 
folders=c("20170518") 
for(f in folders){ 

    #read all shp from that folder, both points and tracks 
    pointfiles <- list.files(paste("./",f,"/points", sep=""), pattern="*.shp$") 
    trackfiles <- list.files(paste("./",f,"/tracks", sep=""), pattern="*.shp$") 

    #for each point and track couple 
    for(i in 1:length(pointfiles)){ 


    #create a linnet object with the track 

    t<- as.linnet(rgdal::readOGR(paste("./",f,"/tracks/",trackfiles[i], sep=""))) 
    #plot(t) 



    #create a ppp object for each set of points 

    pre_p<-rgdal::readOGR(paste("./",f,"/points/",pointfiles[i], sep="")) 
    #plot(p) 
    #obtain the coordinates the current set of points 
    c<-coordinates(pre_p) 
    #create vector of x coords 
    xc=c() 
    #create vector of y coords 
    yc=c() 
    #not a very good way to fill my vectors but it works for my study area 
    for(v in c){ 
     if(v>4000000){yc<-c(yc,v)} 
     else {if(v<4000000 && v>700000){xc<-c(xc,v)}} 
    } 
    print(xc) 
    print(yc) 
    #create a ppp object using the vectors of x and y coords, and a window object 
    #extracted from my set of points 
    p=ppp(xc,yc,Window(as.ppp(pre_p))) 

    #join them into an lpp object 

    pt <- lpp(p,t) 
    #plot(pt) 

    #analize it with the linearK function, nsim=9 for testing purposes 
    #envelope.lpp is the method for analyzing linear point patterns 
    assign(paste("results",f,i,sep="_"),envelope.lpp(pt, nsim=9, fun=linearK)) 

    }#end for each points & track set 
}#end for each day of study 

このスクリプトは、CSRのためのポイントのそれぞれのカップルをテストし、それぞれの日のために追跡、正常に動作して見ることができるようにたった今。残念ながら、私はまだ報告書や報告書を作成して結果を得ていない(あるいは完全に理解している)とは思っていません。もちろん、私はあなたのアドバイスを使用することができます。これはRでの最初の試行であり、多くの新人ミスが起こるからです。

更新フォルダ構造を持つスクリプトとSHPファイルがhere(113キロバイトのサイズを)見つけることができます