世界メッシュ統計可視化の事例 – Research Institute for World Grid Squares

このページでは作成した世界メッシュ統計を可視化する方法についてR言語用いた実際のソースコードを例に事例的に紹介する。R言語を実行するためにはRの実行環境が必要である。また、統合開発環境としてRStudioがそれぞれののOS用に提供されているので以下のソースコード実行に際して事前に実行環境と開発環境を準備されたい。

Rのライブラリleaflet, mapview, webshot, PhantomJSを用いて世界メッシュ統計を可視化する方法について説明する。先ずは

install.packages("leaflet")
install.packages("mapview")
install.packages("webshot")

を用いてRへleaflet, mapview, webshotのライブラリをインストールする。

このRソースでは、国ごとに与えられる世界メッシュ統計データを読込み、そのファイルに含まれている世界メッシュの緯度と経度から描画範囲を決定した上で、leafletのデータを与え画像としてmapshot関数を用いて pngファイルとしてファイルを出力している。内部的には一度htmlを含むleafletを作成し、このhtmlファイルをpngファイルに変換することで画像を得ている。

このRコードを実行するためには、更にPhantomJSが必要である。PhantomJSはRコマンドラインより

webshot::install_phantomjs(force=T)

と実行することでインストールできる。

ソースコードは以下の通りある。以下のソースコードと同じディレクトリに世界メッシュ研究所のページ->データ->NASA夜間光データプロダクトからダウンロードしたtsvファイルを配置し、Rソースコードを実行することで可視化画像をpng形式で得ることができる。図2にオーストリア、バングラデシュ、台湾の2012年NASA夜間光統計データプロダクトを用いた出力画像の例を示す。これらの画像を作成するためにNASA_AUT_mesh3.tsv、NASA_BGD_mesh3.tsv、NASA_TWN_mesh3.tsvをダウンロードして使用した。

(a) オーストリア

(b) バングラデシュ

図2 NASA夜間光3次メッシュ統計の可視化画像（a）オーストリア、（b）バングラデシュ、（c）台湾

[Rソースコード（create_iamge_NASA_light.r）]

# This source code needs PhantomJS.
# PhantomJS can be installed by using webshot::install_phantomjs().
library(leaflet)
library(mapview)
#color
cols<-colorRamp(c("#000000","white","#faf500"))
pal<-colorNumeric(palette=cols, domain=(0:255))

#homedir<-"~/JST/assistant/H28/Jin/NASA_light"
homedir <- getwd()
targetdir<-dir(path=homedir,pattern="tsv")
for(d in targetdir){
  tt<-unlist(strsplit(d,"[/.]"))
  if(tt[2]=="tsv"){
    infile<-d
    ofile<-paste0(tt[1],".png")
    t<-read.csv(paste0(homedir,"/",infile), sep="\t", header=TRUE)
    longmax<-max(t$long0)
    longmin<-min(t$long1)
    latmax<-max(t$lat0)
    latmin<-min(t$lat1) 
    if(((longmax-longmin)/200*3600/45)>1){
      dlong<-(longmax-longmin)/200
    }else{dlong<-1/3600*45} 
    if(((latmax-latmin)/200*3600/30)>1){
      dlat<-(latmax-latmin)/200
    }else{dlat<-1/3600*30}
    nt<-data.frame(lat=c(), long=c(), median_light=c())

    for(i in 1:200){
      longmin<-min(t$long1)
      for(j in 1:200){
        #cat(sprintf("%f, %f    ", longmin, latmin))
        median_t<-median(t[(t$lat0<latmin+dlat)&(t$lat0>latmin)&(t$long0>longmin)&(t$long0<longmin+dlong),]$light_int)
        #cat(sprintf("%d\n",length(t[t$lat0<latmin+dlat&t$lat0>latmin&t$long0>longmin&t$long0<longmin+dlong,]$light_int)))
        if(!is.na(median_t)){
          #cat(sprintf("%f\n", median_t))
          nt <-rbind(nt,data.frame(lat=c(latmin), long=c(longmin), median_light=c(median_t)))
        }
        longmin<-longmin+dlong
      }
      latmin<-latmin+dlat
    }

    #put map
    map<-leaflet(nt) %>%
    addTiles() %>%
    addProviderTiles(providers$OpenStreetMap) %>%
    addRectangles(~long,~lat,~long+dlong,~lat+dlat,color=~pal(median_light), stroke=FALSE,fillOpacity = 1) %>%
    addLegend(position='bottomleft', pal=pal, values=~median_light) %>%
    fitBounds(min(nt$long), min(nt$lat), max(nt$long), max(nt$lat))
    mapshot(map,file=paste0(getwd(),"/", ofile))
  }
}