渡り鳥は美しく、その存在自体で季節の移り変わりを楽しませてくれます。散歩しているときに見かけると、今年もよく来てくれたなと感慨深いものがあります。

その一方、渡り鳥は鳥インフルエンザを運んでくるため、獣医学的には怖い存在でもあります。

渡り鳥のデータ

渡り鳥の観測データが環境省のHPに掲載されているのでそのデータを用いてGIS解析をしてみたいと思います。以下のリンクにある、数の多い種の飛来数の推移というデータを使います。また、調査地の概要というところに調査地の経度緯度があるのでそのデータを使います。

http://www.env.go.jp/nature/dobutsu/bird_flu/migratory/ap_wr_transit20/index.html

観察データ

場所データ

データプロセシング

データを読み込んで、1つのデータフレームに統合します。

b1 <- readxl::read_xls("hashibiro.xls")
b2 <- readxl::read_xls("hidorigamo.xls")
b3 <- readxl::read_xls("hishikui.xls")
b4 <- readxl::read_xls("hoshihajiro.xls")
b5 <- readxl::read_xls("karugamo.xls")
b6 <- readxl::read_xls("kawaaisa.xls")
b7 <- readxl::read_xls("kinkurohajiro.xls")
b8 <- readxl::read_xls("kogamo.xls")
b9 <- readxl::read_xls("kohakucho.xls")
b10 <- readxl::read_xls("magamo.xls")
b11 <- readxl::read_xls("magan.xls")
b12 <- readxl::read_xls("ohakucho.xls")
b13 <- readxl::read_xls("okayoshi.xls")
b14 <- readxl::read_xls("onagagamo.xls")
b15 <- readxl::read_xls("suzugamo.xls")
b16 <- readxl::read_xls("yoshigamo.xls")
places <- read.csv("place.csv")

library(tidyverse)
birds <- rbind(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16)
birds <- na.omit(birds)
colnames(birds) <- c("date","place","name","number")
birds <- left_join(birds,places)
birds$date <- as.Date(as.character(birds$date),
                      format = "%Y%m%d")
birds$number <- as.numeric(birds$number)

library(sf)
birds <- st_as_sf(birds, coords=c('lon', 'lat'), crs=st_crs(4326))

最後の2行は、GISを扱うためのデータに変換するコードです。環境省のデータには経度緯度しか書かれていませんがJGD2011だと思うので、WGS84と互換性があると推測します。WGS84のEPSGコードは4326なので、CRSには4326を設定します。

birdsは以下のようなデータになっています。

ベースとなる日本地図を読み込みます。これは、以前以下のサイトから都道府県境界のみを残してdissolveしたものです。

出典:国土交通省国土数値情報ダウンロードサイト「国土数値情報 行政区域データ」(国土交通省)(https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_4.html)を加工して作成

jpn <- st_read("C:\\Users\\japan.shp")

さっそく、グラフを書いてみます。

ggplot() +
  geom_sf(data = jpn, 
          fill = "white", #線の中
          color = "cadetblue4", #線自体
          lwd = 1 ) +  #線の太さ
  geom_sf(data = birds,
          aes(size=number,color=name),
          alpha = 0.5)+
  theme_void()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("野鳥飛来数 種別 2019秋-2020春")

渡り鳥の種類が色別になっていますが、重なって隠れてしまっています。

鳥の種類の分別をなくして書いてみます。

ggplot() +
  geom_sf(data = jpn, 
          fill = "white", #線の中
          color = "grey40", #線自体
          lwd = 1 ) +  #線の太さ
  geom_sf(data = birds %>%
            group_by(place) %>% 
            summarise(total = sum(number)),
          aes(size = total, color = total)
          )+
  scale_size_continuous(range = c(1, 10))+
  scale_color_gradient(low = "lightskyblue", high = "deeppink")+
  theme_void()+
  theme(plot.title = element_text(hjust = 0.5)) +
  guides(size = FALSE)+ #legendの選択されたものを消す
  ggtitle("野鳥飛来数 総数 2019秋-2020春")

鳥の種類ごとにどこに多く来ているか知りたいです。まずは何も考えずにfacet_wrapを使います。

ggplot() +
  geom_sf(data = jpn)+
  geom_sf(data = birds, aes(size=number,color=name))+
  facet_wrap(~name)+
  theme_void()+
  theme(legend.position = "none")

種類が多くてよくわからないので、1つづつ表示します。

ggplotをforで回す

ggplotをforで回します。毎回出力させるためには、ggplotを何かにアサインして、print(アサインしたオブジェクト)をします。

library(RColorBrewer)
nb.cols <- length(unique(birds$name))
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

n = 1
for (birdname in unique(birds$name)) {
    gg <- ggplot() +
      geom_sf(data = jpn, 
              fill = "white", #線の中
              color = "grey50", #線自体
              lwd = 0.5 ) +  #線の太さ
    geom_sf(data = birds %>% filter(name == birdname),
            aes(size=number),
            color=mycolors[n])+
    theme_void()+
    ggtitle(birdname)+
    theme(legend.position = "none") 
    
    ggsave(file = paste0(birdname,'.png'), plot = gg)
    print(gg)
    rm(gg)
    n = n + 1
}

gganimateを使って、いつ飛んできているのかを分かりやすくします。なお、gganimateでgeom_sfを使うときはinstall.packages(“transformr”)をしておく必要があります。

上3行のコードは、日付と場所を使ってグループごとの合計を出しています。

total <- birds %>% 
  group_by(date,place) %>% 
  summarise(total = sum(number))

library(gganimate)

gg <- ggplot() +
  geom_sf(data = jpn, 
          fill = "white", #線の中
          color = "grey50", #線自体
          lwd = 0.5 ) +  #線の太さ
  geom_sf(data = total,
          aes(color = total,
              size = total))+
  scale_size_continuous(range = c(1, 10))+
  scale_color_gradient(low = "lightskyblue", high = "lightcoral")+
  theme_void()+
  labs(title = "野鳥飛来数 総数 2019秋-2020春",
       caption ='{frame_time}')+
  theme(legend.position = "none",
        plot.title=element_text(size=22, hjust=0.5, vjust=0.5, colour="grey40"),
        plot.caption =element_text(size=24, hjust = 0.8, vjust=10, color="grey40")) +
  transition_time(date)+
  shadow_wake(wake_length = 0.1, alpha = TRUE)

animate(gg,renderer = gifski_renderer("birds.gif"))

意外と全国均一に飛んできているのですね。もっと、北海道とかに集中しているかと思っていました。11月ころに飛来が多くなっている感じがします。

ナイチンゲールローズダイアグラム(nightingale rose diagram)のアニメーション

月ごとの飛来数を見てみたいと思います。sfの形だとgeom_barでエラーが起きるので、as_tibble()でデータフレームに戻します。

ナイチンゲールローズダイアグラムとは、ナイチンゲールが開発したデータビジュアライゼーションの方法で、イギリス議会に提出してそのデータをもって議会を説得したものです。戦時中の負傷者の数を月ごとに表したもので、医療の予算を増やすことを要求した根拠となるものです。子供の頃に絵本で呼んだナイチンゲールに、ここで出会うとは思いませんでした。

作り方は、棒グラフの横軸を円形にしたものです。普通に棒グラフを作って、coord_polarで円にします。表示したいデータに、サポート列を作ってtransition_revealを使って順番に表示させます。

total_month <- total %>% 
  group_by(strftime(date, "%m" )) %>%
  summarise(total_by_month = sum(total)) %>% 
  as_tibble() %>% 
  select(-geometry)

colnames(total_month) <- c("month","total")

total_month$month <- c("2020_01","2020_02","2020_03","2020_04","2020_05",
                       "2019_09","2019_10","2019_11","2019_12")

total_month <- total_month %>% arrange(month)
total_month$month <- as.factor(total_month$month)
total_month$support <- 1:nrow(total_month)

gg2 <- ggplot(total_month,
       aes(x=month, y=total,fill=month)) + 
  geom_bar(stat = "identity")+
  scale_fill_brewer(palette = "Set2")+
  coord_polar(start = 10)+
  theme_minimal()+
  theme(legend.position = "none",
        axis.text.x = element_text(size = 9, vjust = 5, hjust=5),
        axis.text.y = element_blank(),
        axis.title = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  ggtitle("月ごとの野鳥飛来数")+
  transition_reveal(support)

animate(gg2,renderer = gifski_renderer("birds_monthly.gif"))

10月から11月ころに増加し、そのあとは減少傾向ですね。

ナイチンゲールローズダイアグラムの伸びるアニメーション作成方法

もっと滑らかなアニメーションを作りたいという欲がでてきてしまいました。以下、完成品です。

以下ははっきり言って、作成時間がかかりすぎるのでお勧めしません。

geom_barが伸びていくようなアニメーションを作って軸を回転させれば作れるはずです。ただ、transition_reveal()にはそういった機能がないので以下のフレームワークに収まるように考えなければなりません。

transition_states() +
enter_fade() +
enter_grow()

サポート列として、transition_statesを指定してきます。stateが1づつ変化して、過去のデータは残しつつ、未来のデータを0にします。

total_month1 <- total_month
total_month2 <- total_month
total_month3 <- total_month
total_month4 <- total_month
total_month5 <- total_month
total_month6 <- total_month
total_month7 <- total_month
total_month8 <- total_month
total_month9 <- total_month

total_month1$total <- ifelse(total_month1$support < 2,
                             total_month1$total,
                             0)
total_month2$total <- ifelse(total_month2$support < 3,
                             total_month2$total,
                             0)
total_month3$total <- ifelse(total_month3$support < 4,
                             total_month3$total,
                             0)
total_month4$total <- ifelse(total_month4$support < 5,
                             total_month4$total,
                             0)
total_month5$total <- ifelse(total_month5$support < 6,
                             total_month5$total,
                             0)
total_month6$total <- ifelse(total_month6$support < 7,
                             total_month6$total,
                             0)
total_month7$total <- ifelse(total_month7$support < 8,
                             total_month7$total,
                             0)
total_month8$total <- ifelse(total_month8$support < 9,
                             total_month8$total,
                            0)
total_month9$total <- ifelse(total_month9$support < 10,
                             total_month9$total,
                             0)

total_month1$support <- 1
total_month2$support <- 2
total_month3$support <- 3
total_month4$support <- 4
total_month5$support <- 5
total_month6$support <- 6
total_month7$support <- 7
total_month8$support <- 8
total_month9$support <- 9

total_month_animation <- rbind(total_month1,
                               total_month2,
                               total_month3,
                               total_month4,
                               total_month5,
                               total_month6,
                               total_month7,
                               total_month8,
                               total_month9)

gg3 <- ggplot(total_month_animation,
              aes(x=month, y=total, fill=month)) + 
  geom_bar(stat = "identity")+
  scale_fill_brewer(palette = "Set2")+
  coord_polar(start = 10)+
  theme_minimal()+
  ggtitle("月ごとの野鳥飛来数")+
  theme(legend.position = "none",
        axis.text.x = element_text(size = 9, vjust = 5, hjust=5),
        axis.text.y = element_blank(),
        axis.title = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  transition_states(support) +
  enter_fade() +
  enter_grow()
  
animate(gg3,renderer = gifski_renderer("birds_monthly_2.gif"))
category