3

私は R spatstat パッケージの新しいユーザーであり、owin() を使用して多角形の観測ウィンドウを作成する際に問題があります。コードは次のとおりです。

library("maps")
library ("sp")` 
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T) # This returns a data frame includding x and y components that form a polygon of massachusetts mainland`

mass.win <- owin(poly=data.frame(x=mass.map$x, y=mass.map$y)

if (w.area < 0) stop(paste("Area of​​ polygon is negative -", "maybe traversed in >wrong direction?")) のエラー: TRUE/FALSE が必要な場所に値がありません

ポリゴンの順序を逆にしてみましたが、同じエラーが発生しました。

 mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))

ポリゴンに重複した頂点が含まれています

Polygon is self-intersecting Error in owin(poly = data.frame(x = rev(mass.map$x), y = rev(mass.map$y))) : ポリゴン データに重複した頂点と自己交差が含まれています

次に、map() によって返されるポリゴンは、owin() に渡されることを意図していない可能性があると考えました。そこで、マサチューセッツ州のシェイプ ファイルを読み込んでみました (この時点では完全に推測しています)。

x <- readShapePoly("../Geog/OUTLINE25K_POLY") ## The shape file for MASS, loaded from MassGIS website
mass.poly <- x <- readShapePoly("../Geog/OUTLINE25K_POLY", force_ring=T, delete_null_obj=T) ## I got following error whether or not I used force_ring

mass.owin <- as(mass.poly, "owin") 1006 個のポリゴンをチェックしています...1、ポリゴン 1 には重複した頂点が含まれています [91844 個のエッジを持つポリゴンをチェックしています...] 2、3、.. [etd 1:21: 52] ....10 [etd 36:12] ..... [etd 23:10] ....20 [etd 16:59] ..... [etd 13:22] .... 30 [etd 11:01] ..... [etd 9:21] ....40 [etd 8:06] ..... [etd 7:09] ....50 [etd 6:23] ] ..... [etd 5:46] ....60 [etd 5:15] ...[2449 個のエッジを持つポリゴンをチェック中...] .. [etd 4:49] ....70 [ etd 4:27] ..... [etd 4:07] ....80 [etd 3:50] ..... [etd 3:36] ....90 [etd 3:22] . .... [etd 3:11] ....100 [など。

交差する頂点などについて不平を言うメッセージが表示され、ポリゴンの構築に失敗しました。

問題に関するいくつかのコンテキスト：空間相対リスクの計算、つまり、ケースとコントロールの密度の空間比にspatstatの関数を使用しようとしています。そのためには、観測ウィンドウとそのウィンドウ内のポイント プロットが必要です。マサチューセッツ州周辺の長方形の観測ウィンドウをごまかすこともできますが、それはおそらく海岸近くの値を歪めるでしょう. いずれにせよ、このパッケージで今後行う作業のために、これを正しく行う方法を学びたいと思います。ご協力いただきありがとうございます。

4

3 に答える 3

2

編集 2021-02-04: この時計回り/反時計回りの問題が再び発生し、spatstat の問題に対処する唯一の答えを見つけました。答えはあまり役に立ちませんでした。より幅広い用途でお使いいただけるよう改良しました。

spatstatパッケージは反時計回りowinに指定されたオブジェクト座標を望んでいます。バージョン 1.36-0 のドキュメントからの引用:

単一の多角形: poly が 2 つの列をもつ行列またはデータ フレーム、または同じ長さの 2 つのコンポーネント ベクトル x と y をもつ構造体である場合、これらの値は、ウィンドウに外接する多角形の頂点のデカルト座標として解釈されます。頂点は反時計回りにリストする必要があります。頂点を繰り返さないでください (つまり、最初の頂点を繰り返さないでください)。

library("maps")
library ("sp")
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T)

ここに画像の説明を入力

まず、ポリゴンが時計回りか反時計回りかを調べる必要があります。ここでのアプローチは、答えを見つけるために使用できます。それを関数として定式化します。

#' @title Check whether points for an owin are clockwise
#' @param x a dataframe with x coordinates in the first column and y coordinates in the second. 
#' @details Similarly to owin, the polygon should not be closed
#' @return A logical telling whether the polygon is arranged clockwise.
#' @author The idea has been scavenged from https://stackoverflow.com/a/1165943/1082004

clockwise <- function(x) {
    
  x.coords <- c(x[[1]], x[[1]][1])
  y.coords <- c(x[[2]], x[[2]][1])
  
  double.area <- sum(sapply(2:length(x.coords), function(i) {
    (x.coords[i] - x.coords[i-1])*(y.coords[i] + y.coords[i-1])
  }))
  
  double.area > 0
}

mass.mapここで、座標が時計回りかどうかを確認します。

clockwise(data.frame(x=mass.map$x, y=mass.map$y))
#> TRUE

彼らです。とのrevベクトルを逆にする functionを使用して、座標を反時計回りに回します。xy

mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))
plot(mass.win)

ここに画像の説明を入力

于 2014-03-18T10:27:40.863 に答える
0

私はあなたが抱えている問題に精通しています。私の理解が正しければ、あなたはクラス owin の多角形オブジェクトを取得したいと考えています。次のコードは、必要なものを取得するのに役立ちます。

#Loading the required packages
library (sp)
library(maps)

# Reading the dataset
mass.map <- map("state", "massachusetts:main", fill=T)

IDs <- sapply(strsplit(mass.map$names, ":"), function(x) x[1])
# converting can be done by loading the powerful package "maptools"
library(maptools)
# Converting into SpatialPolygons-object
mass.poly <- map2SpatialPolygons(mass.map, IDs = IDs)

# Converting by coercing into owin-object
mass.owin <- as.owin.SpatialPolygons(mass.poly)
于 2014-01-17T15:40:53.960 に答える