2

で次の問題を解決しようとしていますR。

x と y の 2 つのコンポーネントを持つリスト l で定義されたポリゴン オブジェクトがあります。順序は、ポリゴンのエッジを定義します。

例えば ​​：

l=list(
    x=c(-1.93400738955091,0.511747161547164,1.85047596846401,-1.4963460488281,-1.31613255558929,-0.0803828876660542,1.721752044722,-0.724002506376074,-2.08847609804132,2.13366860069641),
    y=c(-1.02967154136169,1.53216851658359,-1.39564869249673,-1.21266011692921,1.6419616619241,-1.87141898897228,0.946605074767527,1.49557080147009,0.324443917837958,-0.517303529772633)
)
plot(l,type="b",pch=16)
points(l$x[c(10,1)],l$y[c(10,1)],type="b",pch=16)

今私が興味を持っているのは、このポリゴンの外側の境界のみを保持することです (ただし、凸包は保持しません)。次の図は、私が保持したいポイントを強調しています

points(
    x=c(-1.13927707377209,-1.31613255249992,-1.3598262571216,0.511747159281619,0.264900107013767,0.671727215417383,-0.724002505140328,-1.93400738893304,-1.4811931364624,-1.45298543105533,-2.08847609804132,-1.40787406113029,-1.3598262571216,0.278826441754518,1.85047596733123,1.48615105742673,1.48615105742673,2.13366860069641,1.38016944537233,1.38016944537233,1.17232981688283,1.17232981688283,1.72175204307433,0.671727215417383,-1.496346, -0.08038289, -0.2824999),
    y=c(1.13914087952916,1.64196166071069,0.949843643913108,1.53216851597378,1.27360509238768,1.18229006681548,1.49557080106148,-1.02967154055378,-0.972634663817139,-0.525818314106921,0.324443915423533,0.188755761926866,0.949843643913108,-1.30971824545964,-1.3956486896768,-0.59886540309968,-0.59886540309968,-0.517303527559411,-0.367082245352325,-0.367082245352325,0.0874657083966551,0.0874657083966551,0.94660507315481,1.18229006681548,-1.21266,-1.871419,-1.281255),
    pch=16,
    col="red",
    cex=0.75
)

それを簡単に行うためのツールがあるかどうか、私は本当に無知です。私が見つけた最も近いpolysimplifyものは、パッケージ内の関数で、polyclip必要なすべてのポイントを識別しますが、不要なポイント (セグメントが交差する内側のポイント) も出力します。

私は実際に解決策を見つけました（以下）。次の関数は私が望むことを行いますが、なぜ機能するのか (および失敗する可能性があるかどうか) はわかりません。

実際、以下の関数は必要なポイントを正しく識別しますが、それらを間違った順序で出力するため、まだ役に立ちません...

polygon.clean<-function(poly){
  require(polyclip)
  poly.cleaned=polysimplify(poly)
  x=unlist(sapply(poly.cleaned,function(x)x$x))
  y=unlist(sapply(poly.cleaned,function(x)x$y))
  x.src=x[!x%in%x[duplicated(x)]]
  y.src=y[!y%in%y[duplicated(y)]]
  poly.cleaned=poly.cleaned[sapply(poly.cleaned,function(poly.sub,x,y){
    any(poly.sub$x%in%x&poly.sub$y%in%y)
  },x=x.src,y=y.src)]

  x=unlist(sapply(poly.cleaned,function(x){
    res=x$x
    if(length(res)==4){
      res=vector()
    }
    res
  }))
  y=unlist(sapply(poly.cleaned,function(x){
    res=x$y
    if(length(res)==4){
      res=vector()
    }
    res
  }))
  x=c(x,x.src)
  y=c(y,y.src)
  tester=duplicated(x)&duplicated(y)
  x=x[!tester]
  y=y[!tester]
  list(x=x,y=y)
}

plot(l,type="b",pch=16)
points(l$x[c(10,1)],l$y[c(10,1)],type="b",pch=16)
points(polygon.clean(l),pch=16,cex=0.75,col="red")
4

1 に答える 1

2

ルーチンを使用してrgeos、最初にラインストリングを「ノード化」してすべての交点を作成し、次にそれを「ポリゴン化」してから「結合」して内部を分解します。

最初に、SpatialLines最初/最後のポイントが重複したデータのバージョンを作成します。

library(sp)
library(rgeos)
coords = cbind(l$x, l$y); coords=rbind(coords,coords[1,])
s = SpatialLines(list(Lines(list(Line(coords)),ID=1)))

それで：

s_outer = gUnaryUnion(gPolygonize(gNode(s)))

次のようにプロットします。

plot(s,lwd=5)
plot(s_outer, lwd=2,border="red",add=TRUE)

ここに画像の説明を入力

周囲のポリゴンの座標が必要な場合、それらは返されたオブジェクトに含まれており、次の方法で抽出できます。

s_outer@polygons[[1]]@Polygons[[1]]@coords

#                x           y
# [1,]  0.27882644 -1.30971825
# [2,] -0.08038289 -1.87141899
# [3,] -0.28886517 -1.27867953

ポリゴンが 1 つしかないと仮定すると、そうではない可能性があります (線が 8 の字をトレースするとします)。すると、2 つのポリゴンが 1 点で接触することになります。あなたのぎざぎざの線がそのようなことをするのにどれほど自由であるかはわかりません...

于 2016-05-20T07:18:55.387 に答える