R기반 SDM패키지를 이용한 생물종 분포모델 - 2편

1. SDM 패키지 설치 및 입력데이터와 환경데이터 입력하기

1) SDM 패키지에서 받아들이는 입력파일 형태로 GIS shp파일을 변경한다.

  - 다음과 같은 구조가 아니면 오류가 발생한다.

 

2) SDM 패키지를 R프로그램에서 설치한다.

# sdm 패키지 설치
install.packages('sdm')
# or
# devtools::install_github('babaknaimi/sdm')

# 라이브러리 불러오기
library(sdm)

 

3) SDM 패키지를 원할하게 실행하기 위해 부수적인 패키지를 InstallAll()함수로 설치한다.

   SDM패키지를 설치한 후 처음 한번만 실행한다..

InstallAll() 

4) SDM과 생물종 분포모델을 개발하는데 필요한 패키지 불러오기

library(raster)
library(dismo)
library(dplyr)
library(tidyr)
library(mapview)
library(rgeos)

5) 생물종 위치자료 불러오기

# 남색이마잠자리 shp 파일 경로 지정
file <- "./남색이마잠자리_pt.shp"

# shp 파일 불러오기
sp <- shapefile(file)
head(sp)
class(sp)

 

5) 현재 기후환경 데이터 불러오기

# bioclim data 해상도 10분(해상도 300km) 짜리 다운로드하기
bio <- raster::getData('worldclim', var='bio', res=10)

# bioclim data 이름 확인하기(bio1~bio19)
bio names(bio)

# bio 변수 1번인 평균기온 출력하고 그 위에 남색이마잠자리 위치포인터 올리기
plot(bio[[1]])
points(sp)

# bio변수 자료형태 확인하기
class(bio[[1]])

bioclimate data 중 bio1번 자료 위에 생물종 위치자료 중첩하기

 

6) 잠자리 위치가 있는 지역으로 대상지역을 한정하기

# 대상지역 한정하기

# 대상지역 shp파일 경로 지정하기
ext <- "./남색이마잠자리_AOI.shp"

# shape 파일 불러오기
sp_ext <- shapefile(ext) 
head(sp_ext)
class(sp_ext)

# 직접 좌표값을 xmin, xmax, ymin, ymax 순으로 입력해서 사각형 영역 설정
#e <- extent(-25,60,-50,50)
#e

# 마우스를 이용해서 사각형 영역 설정방법
# e <- drawExtent()

# GIS shp파일을 이용해서 사각형 영역 설정방법
e <- extent(sp_ext)
e

# 대상지역 안에 있는 입력데이터만 선택하기
spg <- crop(sp, e)
points(sp,col='red')

# 대상지역 안에 있는 기후환경 자료 생성하기
bioc <- crop(bio, e)

# 기후환경자료와 입력데이터 중첩하여 출력하기
plot(bioc[[1]])
points(sp)

# 대상지역으로 한정된 bioclimate 기후자료 19개를 QGIS에서 사용하기 위해 GIS 래스터 asc파일로 저장하기
for (x in 1:19) {
   print(names(bioc[[x]]))
    writeRaster(bioc[[x]], paste('wc10_asc/','bio_',x,'.asc'), NAflag=-9999, overwrite=TRUE)
}

2. SDM 분포모델 실행하기

3편으로 >>>> 2023.03.03 - [머신러닝 & 딥러닝] - R기반 SDM패키지를 이용하여 생물종 분포모델 - 3편

 

6) 생물종 분포모델 실행하기

 

7) 미래 기후환경 데이터 기반 생물종 분포모델 실행하기

 

8) 현재와 미래의 분포가능성 비교 분석하기

 

 

3. SDM 전체 실행 코드

 

# 남색이마잠자리 분포변화 예측
# Error in .updateGDAL(object, v, cell, band, setminmax) : 
# no longer supported
# devtools::install_url('https://cran.r-project.org/src/contrib/Archive/raster/raster_3.5-15.tar.gz')
# devtools::install_url('https://cran.r-project.org/src/contrib/Archive/dismo/dismo_1.3-5.tar.gz')

# sdm 패키지 설치
# install.packages('sdm')
# or
# devtools::install_github('babaknaimi/sdm')

# 라이브러리 불러오기
library(sdm)
#installAll() # only firest time after installing the sdm package


library(dismo)
library(dplyr)
library(tidyr)
library(mapview)
library(raster)
library(rgeos)


# 작업폴더 설정
setwd("c:/R_work")
getwd()


#남색이마잠자리 shp 파일 지정
file <- "./남색이마잠자리_pt.shp" 

# shape 파일 불러오기
sp <- shapefile(file) 
head(sp)
class(sp)


# bioclim data 해상도 10분(해상도 300km) 짜리 다운로드하기

bio <- raster::getData('worldclim', var='bio', res=10)
bio
names(bio)

# bio 변수 1번인 평균기온 출력하고 그 위에 남색이마잠자리 위치포인터 올리기
plot(bio[[1]])  
points(sp)
class(bio[[1]])


# 대상지역 한정하기
ext <- "./남색이마잠자리_AOI.shp" # shape 파일 지정

sp_ext <- shapefile(ext) # shape 파일 불러오기
head(sp_ext)
class(sp_ext)
# 직접 좌표값을 xmin, xmax, ymin, ymax 순으로 입력해서 사각형 영역 설정
# e <- extent(-25,60,-50,50)
# e
e <- extent(sp_ext)
e

# 마우스를 이용해서 사각형 영역 설정
# e <- drawExtent()

spg <- crop(sp, e)
points(sp,col='red')
bioc <- crop(bio, e)
plot(bioc[[1]])
points(sp)

# 기후자료 GIS 래스터 asc파일로 저장하기
for (x in 1:19) {
  print(names(bioc[[x]]))
  writeRaster(bioc[[x]], paste('wc10_asc/','bio_',x,'.asc'), NAflag=-9999, overwrite=TRUE)
}


#------------
library(usdm)

vif(bioc)

ex <- raster::extract(bioc,spg)
class(ex)
head(ex)
head(spg)

ex_df = as.data.frame(ex) 
ex_df
coordinates(ex_df) <- c('lon','lat')
raster::shapefile(ex_df, "ex_val.shp")

v <- vifstep(ex)

#vifcor
v
bioc <- exclude(bioc, v)

#--------------------
library(sdm)

head(sp)

d <- sdmData(sp~., sp, predictors= bioc, bg = list(method='gRandom',n=1000))
d

getmethodNames()

m <- sdm(sp~., d, methods=c('glm','brt','rf','tree','mars','maxent','svm'), replication=c('sub','boot'),
         test.p=30,n=3, parallelSetting=list(ncore=4,method='parallel'))

m

roc(m[[c(1,7,13,19,25,31,37)]])
roc(m[[c(4,10,16,22,28,34,40)]])
#m@models$species$rf$`13`@object

gui(m)

p <- predict(m, bioc,filename='pr_Brachydiplax_bioc.img' ,overwrite=TRUE)
p
names(p)
plot(p[[c(1,7,13,19,25,31,37)]])
plot(p[[c(4,10,16,22,28,34,40)]])

p1 <- predict(m, bio,filename='pr_Brachydiplax.img' ,overwrite=TRUE)
p1
names(p1)
plot(p1[[c(1,7,13,19,25,31,37)]])

p2 <- predict(m, bio,filename='pr_Brachydiplax.img' ,mean=T ,overwrite=TRUE)
p2
names(p2)
plot(p2)

#en1 <- ensemble(m, bio, filename='en.img',setting=list(method='weighted',stat='tss',opt=2))
en <- ensemble(m, p, filename='en_Brachydiplax_bioc.img',setting=list(method='weighted',stat='tss',opt=2))
cl <- colorRampPalette(c('#3E49BB','#3498DB','yellow','orange','red','darkred'))
plot(en , col=cl(200))


##################
biof <- raster::getData('CMIP5', var='bio', res=10, rcp=85, model='CN', year=70)
biof
plot(biof[[1]])
names(biof) <- names(bio)

biofc <- crop(biof, e)


for (x in 1:19) {
  print(names(biofc[[x]]))
  writeRaster(biofc[[x]], paste('future_CMIP5_asc/','bio_',x,'.asc'), NAflag=-9999, overwrite=TRUE)
}


en1 <- ensemble(m, biofc, filename='enf_future.img',setting=list(method='weighted',stat='tss',opt=2))
#--------------
plot(en1, col=cl(200))



cl <- colorRampPalette(c('#3E49BB','#3498DB','yellow','orange','red','darkred'))
#------
plot(en, col=cl(200))
plot(en1, col=cl(200))

proj4string(spg) <- projection(en)
library(mapview)
mapview(en,col.regions=cl(200)) + spg

#-----------
ch <- en1 - en
cl <- colorRampPalette(c('red','orange','yellow','gray','green','blue'))
plot(ch,col=cl(200))


writeRaster(ch, filename = "chagne_results.tif", format="GTiff", overwrite=T)