Rで生態学データのGIS
Rで生態学データのGIS
pointデータからバッファを発生させたり、それを植生データと交差させたり。QGISでもできますが、何回もやり直したり、再現性がないのが嫌だったので、sfパッケージでやってみました。全国データを扱っているので、UTM座標系ごとにデータを分割したり、ちょっとややこしいコードになっています。
※GISは専門外なので、自己責任でお願いします。
# 2022.7.8
# パッケージのインストール(Japanを選択)、初回のみinstall.packages("sf")install.packages("ggplot2")
# パッケージの読み込みlibrary(sf) # GIS用library(ggplot2) # 作図用
# シェープファイルの読み込み# ファイル名を適宜変えてくださいsty <- st_read("C:/Users/NK/Research/Moni1000/Site/Japan/Random500.shp"); stysty <- sty %>% st_transform(4326) # WGS84. 地理座標系（緯度経度）に変換colnames(sty) <- c("id", "X", "Y", "Site", "geometry") # 変数名を変えることもできますggplot(sty) + geom_sf() # ggplotできちんと読み込めているか確認したい場合
# UTM座標に応じてポイントデータを分割します（この後でバッファ面積をより正確に計算するため）sty_N51 <- st_crop(sty, xmin=0, xmax=126, ymin=-90, ymax=90)sty_N52 <- st_crop(sty, xmin=126, xmax=132, ymin=-90, ymax=90)sty_N53 <- st_crop(sty, xmin=132, xmax=138, ymin=-90, ymax=90)sty_N54 <- st_crop(sty, xmin=138, xmax=144, ymin=-90, ymax=90)sty_N55 <- st_crop(sty, xmin=144, xmax=180, ymin=-90, ymax=90)
# UTM座標系にデータを変換しますsty_n51 <- sty_N51 %>% st_transform(3097)sty_n52 <- sty_N52 %>% st_transform(3098)sty_n53 <- sty_N53 %>% st_transform(3099)sty_n54 <- sty_N54 %>% st_transform(3100)sty_n55 <- sty_N55 %>% st_transform(3101)
# 各ポイントから半径1kmのバッファを作成しますbuf_n51 <- st_buffer(sty_n51, dist=1000)buf_n52 <- st_buffer(sty_n52, dist=1000)buf_n53 <- st_buffer(sty_n53, dist=1000)buf_n54 <- st_buffer(sty_n54, dist=1000)buf_n55 <- st_buffer(sty_n55, dist=1000)
# 植生データを読み込みます（ファイル名を適宜変えてください）。ついでにUTM変換も。# もし植生データがUTMごとに分かれていない場合は、ポイントデータと同じようにst_cropで行けると思います（多分）。vg51 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM51.shp") %>% st_transform(3097)vg52 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM52.shp") %>% st_transform(3098)vg53 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM53.shp") %>% st_transform(3099)vg54 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM54.shp") %>% st_transform(3100)vg55 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM55.shp") %>% st_transform(3101)
# 交差を行って、バッファと重なる植生データをくりぬきます# attribute variables are assumed to be spatially constant throughout all geometries という警告が出ますが多分問題ないです# vg51ではなく、st_buffer(vg51, 0)とするのはちゃんと意味がありますが、忘れました・・・（思い出したら追記します）。int_n51 <- st_intersection(buf_n51, st_buffer(vg51, 0))int_n52 <- st_intersection(buf_n52, st_buffer(vg52, 0))int_n53 <- st_intersection(buf_n53, st_buffer(vg53, 0))int_n54 <- st_intersection(buf_n54, st_buffer(vg54, 0))int_n55 <- st_intersection(buf_n55, st_buffer(vg55, 0))
# ポリゴンの重なりなどのエラーがないか確認しますwhich(st_is_valid(int_n51)==FALSE)which(st_is_valid(int_n52)==FALSE)which(st_is_valid(int_n53)==FALSE)which(st_is_valid(int_n54)==FALSE)which(st_is_valid(int_n55)==FALSE)int_n51 <- st_make_valid(int_n51) # エラーがあれば修正を試みます
# ディゾルブをして、同じ植生のポリゴンをくっつけます（目的によってはこの作業は不要かも）int_n51s <- int_n51 %>% group_by(id, MAJOR1) %>% summariseint_n52s <- int_n52 %>% group_by(id, MAJOR1) %>% summariseint_n53s <- int_n53 %>% group_by(id, MAJOR1) %>% summariseint_n54s <- int_n54 %>% group_by(id, MAJOR1) %>% summariseint_n55s <- int_n55 %>% group_by(id, MAJOR1) %>% summarise
# 各ポリゴンの面積を計算しますint_n51s$area <- st_area(int_n51s)int_n52s$area <- st_area(int_n52s)int_n53s$area <- st_area(int_n53s)int_n54s$area <- st_area(int_n54s)int_n55s$area <- st_area(int_n55s)
# データの保存 (QGISとかで確認しておくことを強くお勧めします。目視大事)write_sf(int_n51s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn51.shp", layer_options="ENCODING=UTF-8")write_sf(int_n52s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn52.shp", layer_options="ENCODING=UTF-8")write_sf(int_n53s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn53.shp", layer_options="ENCODING=UTF-8")write_sf(int_n54s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn54.shp", layer_options="ENCODING=UTF-8")write_sf(int_n55s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn55.shp", layer_options="ENCODING=UTF-8")
# パッケージのインストール(Japanを選択)、初回のみinstall.packages("sf")install.packages("ggplot2")
# パッケージの読み込みlibrary(sf) # GIS用library(ggplot2) # 作図用
# シェープファイルの読み込み# ファイル名を適宜変えてくださいsty <- st_read("C:/Users/NK/Research/Moni1000/Site/Japan/Random500.shp"); stysty <- sty %>% st_transform(4326) # WGS84. 地理座標系（緯度経度）に変換colnames(sty) <- c("id", "X", "Y", "Site", "geometry") # 変数名を変えることもできますggplot(sty) + geom_sf() # ggplotできちんと読み込めているか確認したい場合
# UTM座標に応じてポイントデータを分割します（この後でバッファ面積をより正確に計算するため）sty_N51 <- st_crop(sty, xmin=0, xmax=126, ymin=-90, ymax=90)sty_N52 <- st_crop(sty, xmin=126, xmax=132, ymin=-90, ymax=90)sty_N53 <- st_crop(sty, xmin=132, xmax=138, ymin=-90, ymax=90)sty_N54 <- st_crop(sty, xmin=138, xmax=144, ymin=-90, ymax=90)sty_N55 <- st_crop(sty, xmin=144, xmax=180, ymin=-90, ymax=90)
# UTM座標系にデータを変換しますsty_n51 <- sty_N51 %>% st_transform(3097)sty_n52 <- sty_N52 %>% st_transform(3098)sty_n53 <- sty_N53 %>% st_transform(3099)sty_n54 <- sty_N54 %>% st_transform(3100)sty_n55 <- sty_N55 %>% st_transform(3101)
# 各ポイントから半径1kmのバッファを作成しますbuf_n51 <- st_buffer(sty_n51, dist=1000)buf_n52 <- st_buffer(sty_n52, dist=1000)buf_n53 <- st_buffer(sty_n53, dist=1000)buf_n54 <- st_buffer(sty_n54, dist=1000)buf_n55 <- st_buffer(sty_n55, dist=1000)
# 植生データを読み込みます（ファイル名を適宜変えてください）。ついでにUTM変換も。# もし植生データがUTMごとに分かれていない場合は、ポイントデータと同じようにst_cropで行けると思います（多分）。vg51 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM51.shp") %>% st_transform(3097)vg52 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM52.shp") %>% st_transform(3098)vg53 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM53.shp") %>% st_transform(3099)vg54 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM54.shp") %>% st_transform(3100)vg55 <- st_read("C:/Users/NK/Research/Moni1000/Site/vg/vg_UTM55.shp") %>% st_transform(3101)
# 交差を行って、バッファと重なる植生データをくりぬきます# attribute variables are assumed to be spatially constant throughout all geometries という警告が出ますが多分問題ないです# vg51ではなく、st_buffer(vg51, 0)とするのはちゃんと意味がありますが、忘れました・・・（思い出したら追記します）。int_n51 <- st_intersection(buf_n51, st_buffer(vg51, 0))int_n52 <- st_intersection(buf_n52, st_buffer(vg52, 0))int_n53 <- st_intersection(buf_n53, st_buffer(vg53, 0))int_n54 <- st_intersection(buf_n54, st_buffer(vg54, 0))int_n55 <- st_intersection(buf_n55, st_buffer(vg55, 0))
# ポリゴンの重なりなどのエラーがないか確認しますwhich(st_is_valid(int_n51)==FALSE)which(st_is_valid(int_n52)==FALSE)which(st_is_valid(int_n53)==FALSE)which(st_is_valid(int_n54)==FALSE)which(st_is_valid(int_n55)==FALSE)int_n51 <- st_make_valid(int_n51) # エラーがあれば修正を試みます
# ディゾルブをして、同じ植生のポリゴンをくっつけます（目的によってはこの作業は不要かも）int_n51s <- int_n51 %>% group_by(id, MAJOR1) %>% summariseint_n52s <- int_n52 %>% group_by(id, MAJOR1) %>% summariseint_n53s <- int_n53 %>% group_by(id, MAJOR1) %>% summariseint_n54s <- int_n54 %>% group_by(id, MAJOR1) %>% summariseint_n55s <- int_n55 %>% group_by(id, MAJOR1) %>% summarise
# 各ポリゴンの面積を計算しますint_n51s$area <- st_area(int_n51s)int_n52s$area <- st_area(int_n52s)int_n53s$area <- st_area(int_n53s)int_n54s$area <- st_area(int_n54s)int_n55s$area <- st_area(int_n55s)
# データの保存 (QGISとかで確認しておくことを強くお勧めします。目視大事)write_sf(int_n51s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn51.shp", layer_options="ENCODING=UTF-8")write_sf(int_n52s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn52.shp", layer_options="ENCODING=UTF-8")write_sf(int_n53s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn53.shp", layer_options="ENCODING=UTF-8")write_sf(int_n54s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn54.shp", layer_options="ENCODING=UTF-8")write_sf(int_n55s, "C:/Users/NK/Research/Moni1000/Site/Japan/int_UTMn55.shp", layer_options="ENCODING=UTF-8")
自分の記録用ですが、ご自由にお使いください (使用によって生じた如何なる損害も、当方では一切の責任を負いかねます)。もしお役に立った場合はご一報いただけると嬉しいですし、共同研究に発展することがあればさらに嬉しいです。。。
解析のご質問やご相談は katayama6 at affrc.go.jp までご連絡いただければ幸いです。