一般社団法人世界メッシュ研究所

實例

本頁面將介紹的是，根據包含位置信息(緯度和經度)在內的數據（點數據）製作網格統計的方法，並且以使用R言語所做成的世界網格統計可視化方法之實際原始編碼為例介紹給大家。為了執行R語言，R語言的執行環境是必要的。並且、由於RStudio是作為統合開發環境提供給每個操作系統使用的，所以，在執行以下的原始編碼時希望能預先準備好執行環境和開發環境。

世界網格統計的做法


圖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)
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原始代碼所輸出的結果如下所示，形成世界網格編碼、南側經度、西側緯度、北側經度、東側緯度、徵才廣告數。這個原始代碼的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執行了

webshot::install shantomjs()

後就可以安裝了。原始代碼如下所示。與以下的原始代碼在同一個目錄上的世界網格研究所的頁面->數據->NASA夜間光照的數據產品下載tsv檔案並配置，執行R原始代碼可以得出png形式的可視化圖像。圖2為使用奧地利、孟加拉、台灣的2012年NASA夜間光照統計數據產品而得出的圖像例。為了做成這些圖像，我們下載並使用NASA_AUT_mesh3.tsv、NASA_BGD_mesh3.tsv、NASA_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))
  }
}