世界メッシュ研究所

事例

このページでは位置情報(緯度と経度)を含むデータ(ポイントデータ)から世界メッシュ統計を作成したり、作成した世界メッシュ統計を可視化する方法についてR言語用いた実際のソースコードを例に事例的に紹介する。R言語を実行するためにはRの実行環境が必要である。また、統合開発環境としてRStudioがそれぞれののOS用に提供されているので以下のソースコード実行に際して事前に実行環境と開発環境を準備されたい。


世界メッシュ統計の作り方


図1 世界メッシュ統計の作成手順の概念図

世界メッシュ統計を作成する方法について概説する。世界メッシュ統計を作成するための手順として4のステップが存在する。緯度経度情報を含むポイントデータを世界メッシュ統計に変換するためには、各レコードに含まれる緯度経度情報を基に世界メッシュコードを算出し、この世界メッシュコードを基に集計あるいは統計値計算をメッシュコードごとに行えばよい。図1にメッシュ統計を作成する概念図を示す。それぞれの手順は以下のような作業として分解される。

獲得:位置情報を含むデータ源を特定しデータのデータフォーマットの理解、ならびに、特定領域の専門知識を理解することよるデータが表現する文脈の理解が挙げられる。この時、データ源のデータ利用範囲を判明させ、許諾書や契約書を通じたデータの二次利用に関する契約作業も行われる。

収集:データの収集とは特定したデータ源からデータを繰り返し集め、csvファイルやデータベース上に格納していつでも利用できる形にデータを成型する作業を指す。データ収集では、収集されたデータのエラーや抜け、破損などがないかの検査する作業を含む。

エンコーディング:収集された位置情報付きのデータを用いることにより、データレコードに含まれる位置情報(緯度と経度)から世界メッシュコードを算出変換することができる。そして、各レコードに含まれる位置情報から計算される世界メッシュコードを各データレコードに付与する作業がこれに対応する。世界メッシュに含まれるかのラベリング作業を行う。

総和計算:エンドーディングで行われた各レコードに付与した世界メッシュコードを用いて、集計作業や統計処理の演算を行い統計データプロダクトの出力を行う。

緯度経度を含むポイントデータから世界メッシュ統計を作成することができるR言語のコードを以下に示す。ここで、このcsvファイル20170531.csv がRソースコードと同じディレクトリ上に設置されていると仮定する。

source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
mydir<-getwd()
infile<- paste0(mydir,"/20170531.csv")
ofile<- paste0(mydir,"/out_mesh3.csv")
a<-read.csv(infile,fileEncoding="UTF-8",header=F)
b<-head(a,100) # 最初の100行だけで世界メッシュ統計を算出してみる
mm<-c() for(i in 1:nrow(b)){ if(nchar(as.character(b[i,]$V13))>4){
      lat0<-unlist(strsplit(x=as.character(b[i,]$V13),split='[NS.]'))
      long0<-unlist(strsplit(x=as.character(b[i,]$V14),split='[WE.]'))
      lat<-as.numeric(lat0[2])+as.numeric(lat0[3])/60+as.numeric(lat0[4])/3600
      long<-as.numeric(long0[2])+as.numeric(long0[3])/60+as.numeric(long0[4])/3600
      wmeshcode <- cal_meshcode(lat,long)
      cat(sprintf("%s,%s,%s,%s,%s,%s,%s\n",wmeshcode,b[i,]$V1,b[i,]$V2,b[i,]$V3,b[i,]$V4,b[i,]$V5,b[i,]$V6))
      mm[i]<-wmeshcode
  }
}
wm<-unique(mm)
cat(file=ofile,sprintf("worldmeshcode,lat0,long0,lat1,long1,cnt\n"),append=F)
for(wmm in wm){
    cnt<-length(mm[mm==wmm&!is.na(mm)])
    res<-meshcode_to_latlong_grid(wmm)
    cat(file=ofile,sprintf("%s,%f,%f,%f,%f,%d\n",wmm,res$lat0,res$long0,res$lat1,res$long1,cnt),append=T)
}

Rソースコードを実行すると, 出力ファイルout_mesh3.csvには以下のように世界メッシュコード、南側経度、西側緯度、北側経度、東側緯度、求人広告数が出力される。このソースコードのhead(a,100)の数字を100より大きくすると、より多くのレコードを読み込んで求人広告数に関するメッシュ統計を作成することができる。

[出力結果例]

worldmeshcode,lat0,long0,lat1,long1,cnt
2053394526,35.691667,139.700000,35.683333,139.712500,2
2053393586,35.658333,139.700000,35.650000,139.712500,3
2053394601,35.675000,139.762500,35.666667,139.775000,2
2053394577,35.733333,139.712500,35.725000,139.725000,2
2053394611,35.683333,139.762500,35.675000,139.775000,2
2053393578,35.650000,139.725000,35.641667,139.737500,1
2053403039,35.616667,140.112500,35.608333,140.125000,2
2053394578,35.733333,139.725000,35.725000,139.737500,1
2053394507,35.675000,139.712500,35.666667,139.725000,1
2053394579,35.733333,139.737500,35.725000,139.750000,1
2053393566,35.641667,139.700000,35.633333,139.712500,1
2053394528,35.691667,139.725000,35.683333,139.737500,1
2053394672,35.733333,139.775000,35.725000,139.787500,1
2053396249,35.875000,139.362500,35.866667,139.375000,1
2054404345,36.375000,140.437500,36.366667,140.450000,1
2053393596,35.666667,139.700000,35.658333,139.712500,2
2053396522,35.858333,139.650000,35.850000,139.662500,1
2053393263,35.641667,139.287500,35.633333,139.300000,1
2053391541,35.458333,139.637500,35.450000,139.650000,1
2053393642,35.625000,139.775000,35.616667,139.787500,2
2053394353,35.716667,139.412500,35.708333,139.425000,1
2054407544,36.625000,140.675000,36.616667,140.687500,1
2053393547,35.625000,139.712500,35.616667,139.725000,1
2053394506,35.675000,139.700000,35.666667,139.712500,1
2053394709,35.675000,139.987500,35.666667,140.000000,1
2053392415,35.516667,139.562500,35.508333,139.575000,1
2053394632,35.700000,139.775000,35.691667,139.787500,9
2053407531,35.950000,140.637500,35.941667,140.650000,1
2053394576,35.733333,139.700000,35.725000,139.712500,1
2053393548,35.625000,139.725000,35.616667,139.737500,1
2053391530,35.450000,139.625000,35.441667,139.637500,1
2053391459,35.466667,139.612500,35.458333,139.625000,2
2053394435,35.700000,139.562500,35.691667,139.575000,1
2053394652,35.716667,139.775000,35.708333,139.787500,1
2053393394,35.666667,139.425000,35.658333,139.437500,1
2053393691,35.666667,139.762500,35.658333,139.775000,1
2053396463,35.891667,139.537500,35.883333,139.550000,1
2053392345,35.541667,139.437500,35.533333,139.450000,1
2053393559,35.633333,139.737500,35.625000,139.750000,1
2053393390,35.666667,139.375000,35.658333,139.387500,1
2053393330,35.616667,139.375000,35.608333,139.387500,1
2053394504,35.675000,139.675000,35.666667,139.687500,1
2053394612,35.683333,139.775000,35.675000,139.787500,1
2053394446,35.708333,139.575000,35.700000,139.587500,1
2053394543,35.708333,139.662500,35.700000,139.675000,1

世界メッシュ統計の可視化方法

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

 

install.packages(“leaflet”)

install.packages(“mapview”)

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

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

webshot::install shantomjs()

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

 


(a) オーストリア

(b) バングラデシュ

(c) 台湾

図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))
  }
}