splancs 패키지로 마크 부착점 과정

이전 기사( R에서 포인트 패턴 분석 , R에서 포인트 패턴 분석 - 분석 영역 마스킹 )에서 spatstat 패키지를 이용한 마크 부착점 과정의 분석을 소개했습니다.
이 방법에서는 원래 분포 자체가 가지고 있는 공간 집적성(공간 분포의 편향)의 영향을 배제할 수 없기 때문에 유적 분포 영역을 추출하여 그 이외의 지역을 마스킹하는 방법으로 원하는 결과 를 계산했습니다.
이번에는 splancs 패키지로 대상 지역의 균일성을 가정하지 않는 방법을 소개합니다.

벌의 목격 지점 분포

사용하는 데이터는 2017년과 2018년의 아시나가바치와 말벌의 목격 지점의 데이터입니다. sf 패키지에서 GIS 데이터(shp 파일)를 읽고 ggplot으로 시각화합니다.

# パッケージ読み込み
library(tidyverse)
library(sf)
library(ggthemes)
library(spatstat)
library(splancs)
# 地図データ読み込み
data<-st_read(dsn="shp",layer="horne_point")
#地図データ読み込み_厚沢部町字名
assabu<-st_read(dsn="shp",layer="assabu_area")
# 座標値を抜き出す
coord<-st_coordinates(data)%>%as.data.frame()
# 目撃位置マップ
p<-data%>%
    filter(type=="スズメ" | type=="アシナガ")%>%
        ggplot()+
            geom_sf(data=assabu,fill="white")+
            geom_sf(aes(colour=type),alpha=0.4)+
            geom_density2d(data=coord,aes(X,Y),colour="grey")+
            scale_colour_ptol()+
            guides(size=FALSE)+
            theme_map() +
            theme(axis.title = element_blank(),
                axis.text= element_blank(),
                legend.title=element_blank())
print(p)

크게 3개의 클러스터를 볼 수 있습니다. 3개의 클러스터는 홋카이도 아쓰자와베쵸의 혼마치 신마치 시가지, 鶉 지구, 관 지구에 대응하고 있습니다.

아시나가바치와 말벌의 공간 집적성

spatstat 패키지를 이용한 커널 밀도 추정도를 작성합니다.

# pointをデータフレームに変換
point_df<-data%>%as.data.frame()
# ハチデータからxy座標を取り出す
coord<-st_coordinates(data)%>%as.data.frame()
# point_dfとcoordを結合
point_join<-bind_cols(point_df,coord)
# スズメとアシナガだけを選択
point_join2<-point_join%>%
    filter(type=="スズメ" | type=="アシナガ")
# pppオブジェクトを作成
# c(),c()は分析対象領域のx座標とy座標の最小値、最大値
point_p<-ppp(point_join2$X,point_join2$Y,
    c(point_join2$X%>%min(),point_join2$X%>%max()),
    c(point_join2$Y%>%min(),point_join2$Y%>%max()),
    marks=factor(point_join2$type))
# カーネル密度推定図
point_p%>%split()%>%density()%>%plot(1)

삼나무와 참새의 분포는 대개 동일하지만, 참새는 3개의 분포가 거의 같은 집적성을 나타내는 반면, 삼나무에서는 우상(鶉地区)의 집적성이 낮은 것 같습니다.

splancs 패키지에 의한 공간 점 과정

splancs 패키지의 khat() 함수에 의해 아시나가바치와 말벌의 분포(K 통계량)와 그 차이(아시나가 K-참새 K=D(h))를 반복해서 산출하는 몬테카를로 시뮬레이션을 실시합니다. 몬테카를로 시뮬레이션은 kenv.label() 함수입니다.

# splancsパッケージでマーク付点過程
# ハチデータをデータフレームに変換
point_df<-data%>%as.data.frame()
# ハチデータからxy座標を取り出す
coord<-st_coordinates(data)%>%as.data.frame()
# point_dfとcoordを結合_大文字のXYを小文字にする
point_join<-bind_cols(point_df,coord)%>%
    select(x = X,y=Y, everything())
# スズメバチだけを抽出
hor<-point_join%>%
    filter(type=="スズメ")
# アシナガバチだけを抽出
bee<-point_join%>%
    filter(type=="アシナガ")
# polyを算出
point_poly<-list(x=c(hor$x,bee$x),y=c(hor$y,bee$y))
# bboxを作成
point_bb<-bboxx(bbox(as.points(point_poly)))
# シミュレーションプロットのx軸を設定
# シミュレーション距離を10000mにして100mステップで分析
s<-seq(0,10000,100)
# khat関数適用
hor_hat<-khat(as.points(hor),point_bb,s)
bee_hat<-khat(as.points(bee),point_bb,s)
# k統計量の差分を算出
point_diff<-hor_hat-bee_hat
# エンベローブプロット描画
plot(s,point_diff,xlab="distance",ylab="D(h)",ylim=c(-8*10^7,8*10^7),type="l")
env_lab<-Kenv.label(as.points(hor),as.points(bee),point_bb,nsim=99,s=s)
lines(s,env_lab$upper,lty=2)
lines(s,env_lab$lower,lty=2)

점선은 99회 시뮬레이션의 상한값 및 하한값, 실선이 추측값입니다. 거리 3,000m 부근까지는 추측치가 하한값을 밑돌기 때문에, 말벌에 대해 아시나가바치의 공간 집적성이 높은 것을 알 수 있습니다. 3,000m 이내에서는 공간 집적성에 차이가 보이지 않게 되어, 8,000m 부근에서는 바다새벌의 공간 집적성이 오히려 저하합니다.

하우스나 농가 헛간 등 주택지 이외에서도 그 위험성으로부터 목격 정보가 늘어나는 말벌에 대해서, 비교적 위험이 적고, 주택지에서의 목격 건수가 많은 아시나가바치는, 3,000m 이내의 근거리 레인지에서는 말벌보다도 집적성이 높아질 것으로 보입니다. 3,000m를 넘는 레인지에서는, 주택 밀집지와 농지의 거리를 포함해 버려 버리기 때문인지, 집적성의 차이가 눈에 띄지 않게 되어 갑니다.