今回はinaturalistのデータを可視化していきます。
以下のページが参考になったので、まずはこれと同じようにやっていきます。
まずは必要なパッケージ・ライブラリをダウンロード。
#install.packages("rinat")
#install.packages("RCurl")
library(RCurl)
library(rinat)
library(sf)
library(tidyverse)
上述のサイトでは石川県のプロジェクトのデータを使っていましたが、全く同じは面白くないので、今回は大阪のデータを使用します。
#大阪のデータをスクレイピング
obs<-get_inat_obs_project("osaka-wildlife", type = "observations")
本当は20000データ以上ありましたが、サーバーの負荷を抑えるため一日10000データまでに制限されているみたいです。
列の量が非常に多くて扱いづらいので、列を絞ります。
obs_select<-obs %>%
select(id,observed_on,latitude,longitude,species_guess,iconic_taxon_name,user_login,quality_grade,taxon_geoprivacy)
#次に使用する時のためにCSVに保存しておく
write.csv(obs_select,"data_osaka.csv")
分類別に集計してデータ数を見てみます。
#分類別のデータ数
obs_taxon<-obs_select %>%
group_by(iconic_taxon_name) %>%
summarize(count = n()) %>%
as.data.frame()
obs_taxon
#分類別にデータ数をプロット
g_taxon <- ggplot(data=obs_select)+
geom_bar(aes(x = iconic_taxon_name))
g_taxon
#年の列を追加し、年毎に集計
YM<-substring(obs_select$observed_on,1,4)
obs_per_day<-obs_select %>%
bind_cols(YM) %>%
group_by(...10) %>%
summarize(count = n()) %>%
rename(year=...10) %>%
as.data.frame()
g_year<-ggplot(data=obs_per_day,aes(x = year,y = count)) +
geom_bar(stat = "identity")
g_year
#ユーザーで集計
obs_ID<-obs_select %>%
group_by(user_login) %>%
summarize(Number_of_observations = n()) %>%
as.data.frame()
#ユーザー別データ数のヒストグラムを表示
g_ID <- ggplot(obs_ID, aes(x = Number_of_observations)) +
geom_histogram()
g_ID
さあいよいよ地図にプロットしていきます。
#ポイントデータをGISデータにする
obs_shp<-st_as_sf(obs_select, coords = c("longitude", "latitude"), crs = 4326)
#CSVで保存したデータを読み込む場合
osaka_data<-read.csv("C:/Users/ユーザー名/Documents/作業ディレクトリ/data_osaka.csv")
obs_shp<-st_as_sf(osaka_data, coords = c("longitude", "latitude"), crs = 4326)
国土交通省のサイトから大阪の行政区域のデータをダウンロード。
osaka<-st_read("C:/Users/ユーザー名/Documents/作業ディレクトリ/N03-20240101_27_GML/N03-20240101_27.shp")
まずは全データをプロット。
g_point_all<-ggplot() +
geom_sf(data=osaka) +
geom_sf(data=obs_shp)
g_point_all
今度は位置情報が曖昧なデータを除いてプロット。
g_point_open<-ggplot() +
geom_sf(data=osaka) +
geom_sf(data=obs_shp %>% filter(taxon_geoprivacy != "obscured"))
g_point_open
分類別に色変えしてプロット。
g_point_open2<-ggplot() +
geom_sf(data=osaka) +
geom_sf(data=obs_shp %>%
filter(taxon_geoprivacy != "obscured"),aes(colour = iconic_taxon_name))
g_point_open2
とりあえずできました。ここまでは参考にしたサイトと同じです。
次にメッシュ毎にデータ数を集計してマッピングしていきます。
こちらのサイトが参考になりました。
メッシュデータをダウンロードできるパッケージをダウンロードします。
install.packages("jpmesh")
library(jpmesh)
2次メッシュ（10km四方）のデータを作成します。
#大阪にしたいのでコードは27、メッシュサイズ(km)は10を指定
mesh10map <- administration_mesh(code = 27, to_mesh_size = 10)
head(mesh10map)
メッシュデータとポイントデータの重なりを抽出します。
points_in_polygons <- st_join(obs_shp, mesh10map, join = st_within)
> points_in_polygons
Simple feature collection with 10000 features and 8 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 135.0343 ymin: 34.23758 xmax: 135.7996 ymax: 35.19845
Geodetic CRS: WGS 84
First 10 features:
id observed_on species_guess iconic_taxon_name user_login quality_grade
1 233032175 2024-08-01 fleabanes and horseweeds Plantae emil85738 needs_id
2 233033300 2024-07-31 多疣壁虎 Reptilia zirui1 research
3 233256192 2024-06-16 Japanese Boxer Mantis Insecta klearad research
4 233355816 2023-04-06 Typical Blues Insecta claresoong needs_id
5 233400577 2024-08-01 Clytini Insecta zabukari needs_id
6 233401158 2024-06-26 Animalia zabukari needs_id
7 233425805 2024-08-03 Serrognathus titanus pilifer Insecta thomasarendt research
8 233456469 2024-08-03 Graphium sarpedon nipponum Insecta gabriele_calabrese research
9 233483216 2024-08-03 Insecta midori_yaroo needs_id
10 233487365 2024-08-03 Insecta midori_yaroo needs_id
taxon_geoprivacy meshcode geometry
1 <NA> 523515 POINT (135.6311 34.83036)
2 open 513574 POINT (135.5768 34.63562)
3 <NA> 523523 POINT (135.4733 34.84915)
4 <NA> 523504 POINT (135.5227 34.68961)
5 <NA> 513542 POINT (135.2948 34.3994)
6 <NA> 513574 POINT (135.5194 34.61039)
7 <NA> 523525 POINT (135.6791 34.89276)
8 open 523504 POINT (135.5281 34.68449)
9 <NA> 513531 POINT (135.1815 34.31709)
10 <NA> 513531 POINT (135.1825 34.31937)
メッシュコードで集計します。
library(dplyr)
count_points <- points_in_polygons %>%
group_by(meshcode) %>%
summarise(count = n())
> count_points
meshcode count geometry
1 513530 6 MULTIPOINT ((135.1044 34.28...
2 513531 170 MULTIPOINT ((135.1464 34.28...
3 513532 23 MULTIPOINT ((135.3024 34.32...
4 513533 6 MULTIPOINT ((135.4019 34.33...
5 513541 145 MULTIPOINT ((135.1827 34.33...
6 513542 203 MULTIPOINT ((135.2573 34.35...
7 513543 345 MULTIPOINT ((135.4063 34.37...
8 513544 856 MULTIPOINT ((135.522 34.379...
9 513545 102 MULTIPOINT ((135.668 34.403...
10 513551 15 MULTIPOINT ((135.2484 34.43...
11 513552 4 MULTIPOINT ((135.2993 34.43...
12 513553 143 MULTIPOINT ((135.38 34.4937...
13 513554 216 MULTIPOINT ((135.513 34.494...
14 513555 189 MULTIPOINT ((135.6341 34.45...
15 513562 3 MULTIPOINT ((135.2904 34.50...
16 513563 112 MULTIPOINT ((135.4346 34.54...
17 513564 265 MULTIPOINT ((135.5047 34.50...
18 513565 46 MULTIPOINT ((135.6411 34.57...
19 513573 387 MULTIPOINT ((135.4147 34.60...
20 513574 610 MULTIPOINT ((135.5443 34.59...
21 513575 214 MULTIPOINT ((135.6268 34.65...
22 523503 284 MULTIPOINT ((135.4998 34.72...
23 523504 1626 MULTIPOINT ((135.5582 34.67...
24 523505 1241 MULTIPOINT ((135.6969 34.74...
25 523513 130 MULTIPOINT ((135.4876 34.77...
26 523514 487 MULTIPOINT ((135.5797 34.75...
27 523515 574 MULTIPOINT ((135.7474 34.77...
28 523523 945 MULTIPOINT ((135.4529 34.91...
29 523524 221 MULTIPOINT ((135.5104 34.90...
30 523525 267 MULTIPOINT ((135.6332 34.91...
31 523532 4 POINT (135.3696 34.9728)
32 523533 34 MULTIPOINT ((135.4945 34.95...
33 523534 38 MULTIPOINT ((135.5754 34.98...
34 523535 17 MULTIPOINT ((135.6742 34.98...
35 523543 4 MULTIPOINT ((135.458 35.046...
36 <NA> 68 MULTIPOINT ((135.4384 35.19...
ジオメトリの列が2つあると結合時にエラーとなるので、一旦データフレーム化し、ジオメトリの列を削除します。
count_points<-as.data.frame(count_points)
kansatsu<-select(.data = count_points, -geometry)
mesh10mapと集計データを、メッシュ番号を基準に結合してkansatsu_meshを作ります。
kansatsu_mesh <- left_join(kansatsu,mesh10map,by = "meshcode")
kansatsu_mesh
> kansatsu_mesh
meshcode count geometry
1 513530 6 POLYGON ((135 34.25, 135.12...
2 513531 170 POLYGON ((135.125 34.25, 13...
3 513532 23 POLYGON ((135.25 34.25, 135...
4 513533 6 POLYGON ((135.375 34.25, 13...
5 513541 145 POLYGON ((135.125 34.33333,...
6 513542 203 POLYGON ((135.25 34.33333, ...
いい感じです。グラフに描画します。
ggplot() +
geom_sf(data = osaka, fill = "white") +
#mesh10mapのメッシュを描画。alpha=0で透過させ、線はグレー。
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
#先ほど追加したメッシュ数の列のデータを基に色付け(fill)
geom_sf(data=kansatsu_mesh,alpha = 0.8,aes(fill = count,geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933")+
#x軸の目盛りを調整
scale_x_continuous(breaks=seq(135.0, 136,0.25))
今度は年別に分けたグラフを作成します。
#データに年の列を追加
points_in_polygons %>%
mutate(year=substring(observed_on,1,4))->
points_in_polygons2
#メッシュコードと年で集計
count_points2 <- points_in_polygons2 %>%
group_by(meshcode,year) %>%
summarise(count = n())
#データフレーム化し、ジオメトリの列を削除
count_points2<-as.data.frame(count_points2)
kansatsu2<-select(.data = count_points2, -geometry)
#mesh10mapと集計データを、メッシュ番号を基準に結合してkansatsu_meshを作成。
kansatsu_mesh2 <- left_join(kansatsu2,mesh10map,by = "meshcode")
#年毎に分けた図をまとめて表示
ggplot() +
geom_sf(data = osaka, fill = "white") +
geom_sf(data = mesh10map, alpha = 0, color = "grey") +
geom_sf(data=kansatsu_mesh2,alpha = 0.8,aes(fill = count,geometry = geometry)) +
scale_fill_continuous(low = "yellow", high = "#009933")+
scale_x_continuous(breaks=seq(135.0, 136,0.5))+
facet_wrap(~ year, nrow = 3) +
#凡例の位置を調整
theme(legend.position = c( .8 , .15 ))
今度は植生データと重ねあわせてみたいと思います。
環境省から5万分の1の植生図をダウンロードします。
2万5000分の1もありましたが、ダウンロード＆結合するシェープファイルが多すぎるのでやめました。
#ダウンロードした植生データを読み込み
vg1<-st_read("C:/Users/ユーザー名/Documents/作業ディレクトリ/vg_27/vg_27/vg_27a.shp")
vg2<-st_read("C:/Users/ユーザー名/Documents/作業ディレクトリ/vg_27/vg3_27/vg3_27a.shp")
vg3<-st_read("C:/Users/ユーザー名/Documents/作業ディレクトリ/vg_27/vg4_27/vg4_27a.shp")
vg4<-st_read("C:/Users/ユーザー名/Documents/作業ディレクトリ/vg_27/vg5_27/vg5_27a.shp")
#ファイルを結合
osaka_vg<-rbind(vg1,vg2,vg3,vg4)
#メッシュデータとポイントデータの重なりを抽出
#結合のために座標系を変更
osaka_vg <- st_transform(osaka_vg, crs = 4326)
#結合
points_in_polygons3 <- st_join(obs_shp, osaka_vg, join = st_within)
points_in_polygons3
#データフレーム化し必要なデータを抽出
library(dplyr)
points_shyokusei <- points_in_polygons3 %>%
as.data.frame() %>%
select(-geometry) %>%
filter(taxon_geoprivacy != "obscured")
#分類と植生で集計
points_shyokusei %>%
group_by(iconic_taxon_name,NAME) %>%
summarise(count = n())->
syukei_shyokusei
#植生別でデータ数を描画
syukei_shyokusei %>%
ggplot() +
geom_bar(aes(x =NAME,y = count,fill=iconic_taxon_name),stat = "identity") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1,size=7))+
labs(x = "植生")
市街地、公園、開放水域の順にデータが多くなってますね。
水田でのデータもそこそこ。
大阪だと大都会のイメージですが、郊外は田畑も広がっているのでデータが上がっているのかもです。