で、たまにしかやらないので、いつもQGISの使い方を忘れてしまう。そこで、Rで土地被覆データの取得をやってしまいたいと思った。
以下の方法でいけそう。
GISデータは、高解像度土地利用土地被覆図ホームページでダウンロードしたものを使用。ファイル名は「LC_N33E139.tif」。
分布データは色々な情報が入っていても良いが、とりあえず、緯度と経度の列がある。１列目が緯度、２列目が経度。ファイル名は「site.csv」。
library(raster) #raster(), SpatialPointsDataFrame()に必要
library(geobuffer) #geobuffer_pts()に必要
# library(mapview) #mapView()に必要
tochi <- raster("LC_N33E139.tif")
dis<-read.csv("site.csv",header=T)
dis_spdf<-SpatialPointsDataFrame(dis[,2:1], proj4string=tochi@crs, dis) #csvにCRSを組み込む；dis[,2:1]に東経・北緯データがある
xy<-data.frame(lon=dis_spdf$e, lat=dis_spdf$n)
buff_25m <- geobuffer_pts(xy, dist_m = 25, output = "sf") #半径25mのバッファ
# mapView(as(buff_25m, "Spatial"), alpha.regions = 0.7) #web地図で閲覧できる
buff_25_data <- extract(tochi, buff_25m, weight=TRUE, df=TRUE) #buff_25mのtochiラスタの頻度
data_25<-data.frame(tapply(buff_25_data$weight, list(buff_25_data$ID, buff_25_data$LC_N33E139), sum))
data_25[is.na(data_25)]<-0 #NAを0に変換
head(data_25)
「X1,X2,X3,,,」が凡例の「#1,#2,#3,,,」に対応するみたい。
上記の方法は、一度、バッファをgeobuffer_pts()でバッファを作成してから、extract()を実施している。
extract()にbuffer=**で指定する方法もあるみたいだけど、上手くいかなかった。
mapView()は初めて知ったけど、web地図に投影できとても便利。plot()で地図を書いて、その上に描画すると動きが遅くなるけど、mapView()はサクサク動いた。
これができるまでに数日を費やした。再来週の学会で使うつもりだったけど、そもそも解析結果がイマイチなことに本日、気づいた。
諦めて系統解析だけにするか、、、こちらもイマイチ。
RでGISに少しは慣れ、キングダム全巻読破したGWになった。あと、嬉しい連絡もあった。
