LoginSignup
3

More than 5 years have passed since last update.

Rでポイントパターン分析-分析エリアのマスキング

Last updated at Posted at 2019-02-23

前回記事(Rでポイントパターン分析)ではspatstatパッケージを利用したポイントパターンのクラスタリング識別とマーク付き点過程を紹介しました。

結果として、擦文の遺跡は縄文の遺跡の近くに立地することが明らかになりました。

この結果にやや不満を感じるのは、遺跡立地そのものが地形的制約をうけるため、マクロな視点では遺跡の立地は集中してしまうと考えられることです。
分析の前提がおかしいということになります。

遺跡の存在するエリアで分析エリアをマスクする

そこで、遺跡が立地するエリアを抽出し、分析領域を限定することにします。次の手順で進めます。

  1. 遺跡から半径1000mのバッファポリゴンを作成する
  2. バッファポリゴンを結合する
  3. バッファポリゴンの範囲内領域として分析を行う

半径1000mのバッファポリゴン

まずは前回分析したところまでデータ読み込みと処理を進めます。

### パッケージの読み込み
library(tidyverse)
library(sf)
library(GISTools)
library(spatstat)
library(ggthemes)
library(maptools)
### データ読み込み
map<-st_read(dsn="data/",layer="kokudo_gyousei_utm54")
point<-st_read(dsn="data/",layer="site")
###擦文と縄文に分類した新たなフィルード(class)を追加
point<-point%>%
    mutate(class = case_when(
        grepl("擦文",period)   == TRUE ~ "sat",
        grepl("縄文",period)   == TRUE ~ "jou"
        ))%>%
    filter(class=="jou" | class=="sat") %>%
    dplyr::select(class)
point_df<-point%>%as.data.frame()   #pointをデータフレームに変換
coord<-st_coordinates(point)%>%as.data.frame()  #pointオブジェクトからxy座標を取り出す
point_join<-bind_cols(point_df,coord)   #point_dfとcoordを結合
point_join2<-point_join%>%filter(!grepl("NaN",geometry))    #NaNを除去(なぜか2行あった)

バッファ円を作成します。
オブジェクト「point」はsfクラスです。少し長いですが、spatstatパッケージで分析領域を定義するowinという種類のオブジェクトに変換するまでの手順となります。

  1. st_buffer(dist=中心点からの距離,endCapStyle=バッファの形状)関数で半径1000mのバッファ円作成
  2. st_union()関数で全てのバッファ円を結合
  3. as("Spatial")でspクラスへ変換(sfクラスのままだとspatstatパッケージで扱えない)
  4. as("owin")でspatstatパッケージのowinオブジェクトに変換
buf<-point%>%st_buffer(dist=1000,endCapStyle = "ROUND")%>%
    st_union()%>%as("Spatial")%>%as("owin")

st_buffer()で作成した半径1000mのバッファ円。スケールが小さすぎるのでドットにしか見えませんが。
01.png

st_union()でバッファ円を結合。これをspオブジェクト形式を通じてowinオブジェクトに変換します。
02.png

pppオブジェクトを作成

ppp()関数の引数windowにowinオブジェクトに変換した結合バッファ円(buf)を指定します。
これで、bufを対象領域とするpppオブジェクトになります。

point_p<-ppp(point_join2$X,point_join2$Y,
    window=buf,
    marks=factor(point_join2$class))

対象領域が設定されているので、密度図は対象領域のみ算出されます。

point_p%>%density()%>%plot()

03.png

エンベローププロットの出力

### cross-L関数によるエンベローププロット
point_p%>%envelope(Lcross,nism=99)%>%plot()

矩形の対象領域を設定したときに比べて、処理時間が大幅に増加します。
3時間ほどかかりました。
Lcross_envelope.png

対象領域が狭いうちは縄文と擦文の遺跡分布は互いに無関係に分布するようですが、対象領域が25kmぐらいになると縄文から離れたところに擦文遺跡が分布する傾向があるようです。

縄文遺跡が密集するいわゆる下海岸(南茅部地区や恵山地区)に擦文遺跡がほとんど見つかっていないことなどが影響しているのかもしれません。
06.png

density.png

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
3