今回は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 = "植生")

市街地、公園、開放水域の順にデータが多くなってますね。
水田でのデータもそこそこ。
大阪だと大都会のイメージですが、郊外は田畑も広がっているのでデータが上がっているのかもです。