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

实例

本页面通过使用R语言的代码的实例，介绍如何从包含位置信息（经度和纬度）在内的数据（记点值数据）制作世界网格统计，以及将世界网格统计可视化的方法。为了运行R语言，需要安装R的运行环境。并且，作为综合开发环境的RStudio提供了可在各OS环境下运行的版本，在运行以下代码之前请先准备好运行环境与开发环境

世界网格统计的制作方法


图1 世界网格统计的制作步骤的概念图

总体介绍一下世界网格统计的制作方法。制作世界网格统计有4个步骤。在将含有经度纬度信息的记点值数据转变为世界网格统计时，从各记录的经纬度信息计算世界网格编码，再在各网格编码中收集计算统计值。图１是制作网格统计的概念图。各个步骤可分解为如下操作。

获取：特定包含位置信息的情报的数据源并取得提供平台的理解，以及，了解特定领域的专业知识以充分理解数据所表现的含义等等。同时，辨明数据源所允许的利用范围，以及进行通过承诺书，合同等取得数据的二次利用的签约。

收集：数据收集是指从特定的数据源中收集数据，并整理归纳到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”)

将leaflet包和mapview包安装到R。

在此R代码中，读取各个国家对应的世界网格统计数据，通过文件中所包含的世界网格的经度与纬度决定绘制区域，使用leaflet的mapshot函数将数据绘制成png文件并输出。在内部先生成一次含有html文件的leaflet，再将其转换为png文件从而获得图像。为了运行此R代码，还需要PhantomJS。PhantomJS可以通过在R的命令栏输入

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