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

『メッシュ統計』読者サポートページ

書籍情報

最終更新日:2022年11月14日 15:00:00

このページは、『メッシュ統計』（共立出版シリーズ統計学One Point15、佐藤彰洋著、2019）の読者サポートページです。お問合せは電子メールでoffice@fttsus.jpまでお願いします。

共立出版 メッシュ統計 紹介ページ:https://www.kyoritsu-pub.co.jp/bookdetail/9784320112667

内容

正誤表

頁行誤正
p.2l.20日本工業規格日本産業規格　(工業標準化法の一部改正により、標準化の対象にデータ、サービスが追加され、法律名を“産業標準化法”に改め、“日本工業規格（JIS）”が“日本産業規格（JIS）”に変わった（令和元年7月1日施行）)
p.9(1.3)式\begin{equation} \mbox{longitude}_e = \textcolor{red}{\mbox{latitude}}_w + 45 \div 3600 \end{equation} \begin{equation} \mbox{longitude}_e = \textcolor{black}{\mbox{longitude}_w} + 45 \div 3600 \end{equation}
p.17l.8得られる区画に経度方向については南から, 経度方向については得られる区画に緯度方向については南から, 経度方向については
p.18l.7北に0, 2, 4, 6,北に0, 2, 4, 6, 8
p.20(2.35)式\begin{equation} \mbox{latitude}_n=(p{\textcolor{red}{\mathbf -}}1)\times 40/60 \end{equation}\begin{equation} \mbox{latitude}_n=(p{\textcolor{black}{\mathbf +}}1)\times 40/60 \end{equation}
p.21(2.39)式\begin{eqnarray} &&\mbox{latitude}_n=p \times 40/60 \\ &+& (q{\textcolor{red}{\mathbf -}}1)\times 5/60 \end{eqnarray}\begin{eqnarray} && \mbox{latitude}_n = p \times 40/60 \\ &+& (q{\textcolor{black}{\mathbf +}}1)\times 5/60 \end{eqnarray}
p.21(2.43)式\begin{eqnarray} &&\mbox{latitude}_n=p \times 40/60 \\ &+& q\times 5/60 \\ &+& (r{\textcolor{red}{\mathbf -}}1)\times 30/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{latitude}_n=p \times 40/60 \\ &+& q\times 5/60 \\ &+& (r{\textcolor{black}{\mathbf +}}1)\times 30/3600 \end{eqnarray}
p.21(2.45)式\begin{eqnarray} &&\mbox{latitude}_s = p \times 40/60 \\ &+& q \times 5/60 \\ &+& r \times 30/3600 \\ &+& \textcolor{red}{{\mathbf (s_2 – 2 \lfloor s_2/2 \rfloor)}} \times 15/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{latitude}_s = p \times 40/60 \\ &+& q \times 5/60 \\ &+& r \times 30/3600 \\ &+& \textcolor{black}{{\mathbf \lfloor (s_2-1)/2 \rfloor}} \times 15/3600 \end{eqnarray}
p.22(2.46)式\begin{eqnarray} &&\mbox{longitude}_w=100+u \\ &+& v \times 7.5/60 \\&+& w \times 45/3600 \\ &+& \textcolor{red}{{\mathbf \lfloor s_2/2 \rfloor}} \times 22.5/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{longitude}_w=100+u \\ &+& v \times 7.5/60 \\ &+& w \times 45/3600 \\ &+& \textcolor{black}{{\mathbf (s_2-1-2\lfloor (s_2-1)/2 \rfloor )}} \\ &\times& 22.5/3600 \end{eqnarray}
p.22(2.47)式\begin{eqnarray} &&\mbox{latitude}_n = p \times 40/60 \\ &+& q \times 5/60 \\ &+& r \times 30/3600 \\ &+& \textcolor{red}{\mathbf (s_2 – 2 \lfloor s_2/2 \rfloor)-1)} \times 15/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{latitude}_n = p \times 40/60 \\ &+& q \times 5/60 \\ &+& r \times 30/3600 \\ &+& \textcolor{black}{\mathbf (\lfloor (s_2-1)/2 \rfloor +1)} \\ &\times& 15/3600 \end{eqnarray}
p.22(2.48)式\begin{eqnarray} &&\mbox{longitude}_e=100+u \\ &+& v \times 7.5/60 \\&+& w \times 45/3600 \\ &+& \textcolor{red}{\mathbf (\lfloor s_2/2 \rfloor + 1)} \times 22.5/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{longitude}_e=100+u \\ &+& v \times 7.5/60 \\ &+& w \times 45/3600 \\ &+& \textcolor{black}{\mathbf (s_2 – 2\lfloor (s_2-1)/2 \rfloor)} \\ &\times& 22.5/3600 \end{eqnarray}
p.22(2.49)式\begin{eqnarray}
&&\mbox{latitude}_s=p\times 40/60 \\ &+& q \times 5/60 + r \times 30/3600 \\ &+& \textcolor{red}{\mathbf ((s_2-1)-2\lfloor (s_2-1)/2 \rfloor)} \times 15/3600 \\ &+& \textcolor{red}{\mathbf (s_4 – 2 \lfloor s_4/2 \rfloor)} \times 7.5/3600
\end{eqnarray}		\begin{eqnarray}
&&\mbox{latitude}_s=p\times 40/60 \\ &+& q \times 5/60 + r \times 30/3600 \\ &+& \textcolor{black}{\mathbf \lfloor (s_2-1)/2 \rfloor} \times 15/3600 \\ &+& \textcolor{black}{\mathbf \lfloor (s_4-1)/2 \rfloor} \times 7.5/3600
\end{eqnarray}
p.22(2.50)式\begin{eqnarray} &&\mbox{longitude}_w=100+u \\ &+& v \times 7.5/60 \\ &+& w \times 45/3600 \\ &+& \textcolor{red}{\mathbf \lfloor (s_2-1)/2 \rfloor} \times 22.5/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{longitude}_w=100+u \\ &+& v \times 7.5/60 \\ &+& w \times 45/3600 \\ &+& \textcolor{black}{\mathbf (s_2-1-2\lfloor (s_2-1)/2 \rfloor )} \\ &\times& 22.5/3600 \\ &+& \textcolor{black}{\mathbf (s_4-1-2\lfloor (s_4-1)/2 \rfloor )} \\ &\times& 11.25/3600 \end{eqnarray}
p.22(2.51)式\begin{eqnarray} &&\mbox{latitude}_n = p \times 40/60 \\ &+& q \times 5/60 + r \times 30/3600 \\ &+& \textcolor{red}{\mathbf ((s_2-1) – 2 \lfloor (s_2-1)/2 \rfloor-1)} \\ &\times& 15/3600 \\ &+& \textcolor{red}{\mathbf (s_4 – 2 \lfloor s_4/2 \rfloor-1)} \\ &\times& 7.5/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{latitude}_n = p \times 40/60 \\ &+& q \times 5/60 + r \times 30/3600 \\ &+& \textcolor{black}{\mathbf \lfloor (s_2-1)/2 \rfloor } \\ &\times& 15/3600 \\ &+& \textcolor{black}{\mathbf (\lfloor (s_4-1)/2\rfloor+1) } \\ &\times& 7.5/3600 \end{eqnarray}
p.22(2.52)式\begin{eqnarray} &&\mbox{longitude}_e=100+u \\ &+& v \times 7.5/60 +w \times 45/3600 \\ &+& \textcolor{red}{\mathbf (\lfloor (s_2-1)/2 \rfloor)} \times 22.5/3600 \\ &+& \textcolor{red}{\mathbf (\lfloor s_4/2 \rfloor + 1)} \\ &\times& 11.25/3600 \end{eqnarray}\begin{eqnarray} &&\mbox{longitude}_e=100+u \\ &+& v \times 7.5/60 +w \times 45/3600 \\ &+& \textcolor{black}{\mathbf (s_2-1-2\lfloor (s_2-1)/2 \rfloor)} \\ &\times& 22.5/3600 \\ &+& \textcolor{black}{\mathbf (s_4-2\lfloor (s_4-1)/2 \rfloor)} \\ &\times& 11.25/3600 \end{eqnarray}
p.32最下段 (iv)(iii)で得られた時間区間\(d_k\), 地域メッシュコード\(c_l\)の直積集合\(S(d_k, c_l)\)上の統計量\(Q(d_k, c_l)\)は時間区間\(d_k\)ごとの空間\(c_l\)における(iii)で得られた時間区間\(d_k\), 地域メッシュコード\(c_l\)の直積集合\(S(d_k, c_l)\)上の統計量\(Q(d_k, c_l)\)は時間区間\(d_k\)ごとの空間\(c_l\)における地域メッシュ統計としてフォーマットしてファイルに出力される.
ソースコード2.2 (p.35)l. 25res<-meshcode_to_latlong_grid(wmm)res<-meshcode_to_latlong_grid(paste0(“20”,wmm))
p.65l.107) 環境保全・自然保護計画, 7) 土地利用計画7) 環境保全・自然保護計画, 8) 土地利用計画
ｐ.70(3.11)式\begin{eqnarray} \tilde{x}(i,j) &=& \frac{1}{5}\Bigl\{x(i,j) + x(i,j-1) \\ &+& x(i-1,j\color{red}-1\color{black}) + x(i\color{red}-\color{black}1,j) \\ &+& x(i\color{red}-1\color{black},j\color{red}-\color{black}1) \Bigr\}\end{eqnarray}\begin{eqnarray} \tilde{x}(i,j) &=& \frac{1}{5}\Bigl\{x(i,j) + x(i,j-1) \\ &+& x(i-1,j) + x(i {\bf\large +} 1,j) \\ &+& x(i,j {\bf\large +} 1) \Bigr\}\end{eqnarray}
p.80l. 12$k$-秘匿性$k$-匿名性
p.103(3.61)式$\Delta g = g – \gamma + \color{red}\delta g_B\color{black} + \delta g_F + T(\rho) + \delta g_B$$\Delta g = g – \gamma + {\bf\it\delta g_A} + \delta g_F + T(\rho) + \delta g_B$
p.137l.17$\color{red}{x}_w$${\bf\it\large y}_w$
p.140l.11GSRGRS
p.141l.27GSRGRS
p.157l.20井出町井手町
p.160表5.3井出町井手町
p.165l.19webshot::install Phantomjs()install.packages(“webshot”)
library(webshot)
webshot:install_phantomjs()

『メッシュ統計』に掲載されているURLのリンク集です。書籍に掲載されているURLを手動で入力しなくてもクリックしてページ閲覧できるようにしています。書籍内に掲載されているURLを確認して使用してください。

第１章　メッシュ統計とは

第２章　地域メッシュコードの定義と地域メッシュ統計の具体例

第３章　地域メッシュ統計の利用方法

第４章　他国におけるグリッド統計の現状

第５章　世界メッシュ統計

第６章　世界メッシュ統計の分析例とワークフロー

Rソースコード

『メッシュ統計』に掲載されているソースコードです。書籍に掲載されているRソースコードを手動で入力しなくてもコピーアンドペーストまたはダウンロードして取り込むことができるようになっています。書籍内に掲載されているRソースコード番号を確認して取り込みして使ってください。

第１章　メッシュ統計とは

Rソースコード 1.1

3次メッシュコードから3次メッシュの代表的空間座標（南側北側緯度と西側東側経度）を算出するコード. (Rソースコード 1.1)


# calculate geographical positions of 3rd level meshcode
# input: keycode: 3rd level meshcode (integer)
# output: latitude_s, longitude_w, latitude_n, longitude_e (float)
meshcode3_to_latlong <- function(keycode){
  p <- floor(keycode/1000000)
  keycode <- keycode-p*1000000
  u <- floor(keycode/10000)
  keycode <- keycode-u*10000
  q <- floor(keycode/1000)
  keycode <- keycode-q*1000
  v <- floor(keycode/100)
  keycode <- keycode-v*100
  r <- floor(keycode/10)
  keycode <- keycode-r*10
  w <- floor(keycode)
  latitude_s <- p*40/60+q*5/60+r*30/3600
  longitude_w <- 100+u+v*7.5/60+w*45/3600
  latitude_n <- latitude_s + 30/3600
  longitude_e <- longitude_w + 45/3600
  res<-data.frame(latitude_s=latitude_s,longitude_w=longitude_w,latitude_n=latitude_n,longitude_e=longitude_e)
  return(res)
}
keycode=c(52351502)
for(k in keycode){
  res <- meshcode3_to_latlong(k)
  cat(sprintf("%d,%f,%f,%f,%f\n",k,res$latitude_s,res$longitude_w,res$latitude_n,res$longitude_e))
}

第２章　地域メッシュコードの定義と地域メッシュ統計の具体例

Rソースコード 2.1

2地点の緯度経度から大圏距離を計算する. (Rソースコード 2.1)


# calculate geodesic distance based on Vincenty formula
# T. Vincenty, ``Direct and Inverse Solutions of Geodesics
# on the Ellipsoid with application of nested equations'',
# Survey Review XXIII, Vol 176 (1975) Vol. 88-93.
#
VincentyWGS84 <- function(latitude1,longitude1,latitude2,longitude2){
  # WGS84
  f <- 1/298.257223563
  a <- 6378137.0 # [m]
  b <- 6356752.314245 # [m]
  return(Vincenty(a,b,f,latitude1,longitude1,latitude2,longitude2))
}
VincentyGRS80 <- function(latitude1,longitude1,latitude2,longitude2){
  # GRS80
  f = 1/298.257222101
  a = 6378137.0 # [m]
  b = 6356752.31414 # [m]
  return(Vincenty(a,b,f,latitude1,longitude1,latitude2,longitude2))    
}
Vincenty <- function(a,b,f,latitude1,longitude1,latitude2,longitude2){
  # a : Semi-Major axis [m]
  # b : Semi-Minor axis [m]
  # f : Flattening 
  L = (longitude1 - longitude2)/180*pi
  U1 = atan((1-f)*tan(latitude1/180*pi))
  U2 = atan((1-f)*tan(latitude2/180*pi))
  lambda = L
  dlambda = 10
  while(abs(dlambda) > 1e-12){
    sinsigma = sqrt((cos(U2)*sin(lambda))**2 + (cos(U1)*sin(U2)-sin(U1)*cos(U2)*cos(lambda))**2)
    cossigma = sin(U1)*sin(U2)+cos(U1)*cos(U2)*cos(lambda)
    sigma = atan(sinsigma/cossigma)
    sinalpha = cos(U1)*cos(U2)*sin(lambda)/sinsigma
    cos2alpha = 1.0 - sinalpha**2
    if(cos2alpha==0.0){
      lambda0 = L + f*sinalpha*sigma
      C = 0.0
    }else{
      cos2sigmam = cossigma - 2*sin(U1)*sin(U2)/cos2alpha
      C = f/16*cos2alpha*(4+f*(4-3*cos2alpha))
      lambda0 = L + (1-C)*f*sinalpha*(sigma + C*sinsigma*(cos2sigmam + C*cossigma*(-1+2*cos2sigmam**2)))
    }
    dlambda = lambda0 - lambda
    lambda = lambda0
  } 
  if(C==0.0){
    A = 1.0
    dsigma = 0.0
  }else{
    u2 = cos2alpha * (a*a-b*b)/(b*b)
    A = 1 + u2/16384*(4096 + u2 * (-768 + u2*(320-175*u2)))
    B = u2/1024*(256+u2*(-128+u2*(74-47*u2)))
    dsigma = B*sinsigma*(cos2sigmam + 1/4*B*(cossigma*(-1+2*cos2sigmam**2)-1/6*B*cos2sigmam*(-3+4*sinsigma**2)*(-3+4*cos2sigmam**2)))
  }
  s = b*A*(sigma-dsigma)
  return(s)
}
#
lat1 = 35.0833333 # (degree)
long1 = 135.75 # (degree)
lat2 = 35.0833333 # (degree)
long2 = 135.875 # (degree)
cat (sprintf("WGS84: s = %f\n",VincentyWGS84(lat1,long1,lat2,long2)))
cat (sprintf("GRS80: s = %f\n",VincentyGRS80(lat1,long1,lat2,long2)))

Rソースコード 2.2

ポイントデータから3次メッシュ統計を計算. (Rソースコード 2.2)


source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
mydir<-getwd()
infile<- paste0(mydir,"/20170531.csv")
a<-read.csv(infile,fileEncoding="UTF-8",header=F)
b<-head(a,100) # The first 100 records are only used for calculation.
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 <- substring(cal_meshcode(lat,long),3,10)
      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)
for(wmm in wm){
    cnt<-length(mm[mm==wmm&!is.na(mm)])
    res<-meshcode_to_latlong_grid(paste0("20",wmm))
    cat(sprintf("%s,%f,%f,%f,%f,%d\n",wmm,res$lat0,res$long0,res$lat1,res$long1,cnt))
}

Rソースコード 2.3

ポリゴンの面積および交差領域を計算するRソースコード. (Rソースコード 2.3)


library(rgeos)
V1 = readWKT("POLYGON((0 0, 3 0, 3 3, 0 3, 0 0))")
V2 = readWKT("POLYGON((1 1, 4 0, 4 4, 0 4, 1 1))")
Q = gIntersection(V1,V2)
SV1=gArea(V1)
SV2=gArea(V2)
SQ=gArea(Q)
cat(sprintf("V1 = %f, V2 = %f, Q = %f\n",SV1,SV2,SQ))

Rソースコード 2.4

leafletを用いた土地利用詳細図の可視化. (Rソースコード 2.4)


library(sf)
library(leaflet)
library(mapview)
mydir <- getwd()
f <- paste0(mydir,"/L03-b-14_5236-jgd_GML/L03-b-14_5236.shp")
# file read
map <- st_read(f, options = "ENCODING=CP932", stringsAsFactors=F)
map5000<-tail(head(map,20000),5000)
# mapping by leaflet
nt<-data.frame(lat1=c(),long1=c(),lat2=c(),long2=c(),val=c())
for(i in 1:nrow(map5000)){
  mesh<-st_bbox(map5000$geometry[i])
  val<-map5000$土地利用種[i]
  nt<-rbind(nt,data.frame(lat1=c(mesh$ymin),long1=c(mesh$xmin),lat2=c(mesh$ymax),long2=c(mesh$xmax),val=c(val)))
}
landtypes <- unique(map$土地利用種)
# create color pallets 
pals<-colorFactor(rainbow(length(landtypes)),domain = landtypes)
#
leaflet(nt) %>%
  addTiles() %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addRectangles(~long1,~lat1,~long2,~lat2,color=~pals(val),fillOpacity=1) %>%
  addLegend(position='bottomleft',pal=pals,values=~val) %>%
  fitBounds(min(nt$long1),min(nt$lat1),max(nt$long2),max(nt$lat2))

Rソースコード 2.5

浸水想定区画のメッシュ変換. (Rソースコード 2.5)


library(sf)
library(lwgeom)
source("https://www.fttsus.jp/worldmesh/R/worldmesh.R")
#
mydir <- getwd()
f <- paste0(mydir,"/A31-12_27_GML/A31-12_27.shp")
ofile <- paste(mydir,"/A31-12_27_mesh3.csv",sep="")
#
a <- st_read(f)
header<-sprintf("#meshcode,lat0,long0,lat1,long1,label\n")
cat(file=ofile,header,append=F)
for(kk in 1:nrow(a)){
    pol <- a$geometry[kk]
    pol_label <- a$A31_001[kk]
    #
    bbox <- attr(pol,"bbox")
    minlat <- bbox$ymin
    maxlat <- bbox$ymax
    minlong <- bbox$xmin
    maxlong <- bbox$xmax
    deltalat=30
    deltalong=45
    nlat <- ceiling((maxlat-minlat)*60*60/deltalat) 
    nlong <- ceiling((maxlong-minlong)*60*60/deltalong) 
    meshlist<-c()
    plot(pol)
    for(j1 in 0:(nlat+1)){
      lat <- minlat + j1*deltalat/60/60
      for(j2 in 0:(nlong+1)){
        long <- minlong + j2*deltalong/60/60
        meshcode<-cal_meshcode3(lat,long)
        res<-meshcode_to_latlong_grid(meshcode)
        mesh <- st_sfc(st_polygon(list(rbind(c(min(res$long0,res$long1),min(res$lat0,res$lat1)),c(max(res$long0,res$long1),min(res$lat0,res$lat1)),c(max(res$long0,res$long1),max(res$lat0,res$lat1)),c(min(res$long0,res$long1),max(res$lat0,res$lat1)),c(min(res$long0,res$long1),min(res$lat0,res$lat1))))))
        st_crs(mesh) <- "+proj=longlat +ellps=GRS80 +no_defs"
        x<-st_intersection(pol,mesh)
        if(length(st_area(x))!=0){
          jmeshcode <- substring(as.character(meshcode),3,10)
          st<-sprintf("%s,%f,%f,%f,%f,%s\n",jmeshcode,res$lat0,res$long0,res$lat1,res$long1,pol_label)
          cat(st)
          cat(file=ofile,st,append=T)
        }
      }
    }
}

第４章　他国におけるグリッド統計の現状

Rソースコード 4.1

ETRS89-LAEAグリッドと世界メッシュとの比較表示. (Rソースコード 4.1)


library(raster)
source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
#
Gen_frame_e_to_w<-function(efgs_code){
  Y <- as.numeric(gsub('.*N([0-9]+)[EW].*', '\\1', efgs_code))*1000
  X <- as.numeric(gsub('.*[EW]([0-9]+)', '\\1', efgs_code))*1000
  frame <- as(extent(X,X+1000,Y,Y+1000), "SpatialPolygons")
  proj4string(frame) <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
  frame_etow <- spTransform(frame,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  return(frame_etow)
}
#
Cal_cover_efgs_to_world<-function(efgs_code){
  frame_etow<-Gen_frame_e_to_w(efgs_code)
  min_long<-xmin(frame_etow)
  min_lat<-ymin(frame_etow)
  max_long<-xmax(frame_etow)
  max_lat<-ymax(frame_etow)
  mid_lat<-(min_lat+max_lat)/2
  mid_long<-(min_long+max_long)/2
  NE=c(max_lat,max_long)
  NW=c(max_lat,min_long)
  SE=c(min_lat,max_long)
  SW=c(min_lat,min_long)
  M1=c(min_lat,mid_long)
  M2=c(max_lat,mid_long)
  M3=c(mid_lat,min_long)
  M4=c(mid_lat,max_long)
  check_points<-rbind(NE,NW,SE,SW,M1,M2,M3,M4)
  cover<-apply(check_points,1,function(p){return(cal_meshcode3(p[1],p[2]))})
  return(unique(cover))
}
#
Plot_cover_and_efgs_code<-function(efgs_code){
  frame_etow <-Gen_frame_e_to_w(efgs_code)
  cover<-Cal_cover_efgs_to_world(efgs_code)
  cover_grids<-lapply(cover,meshcode_to_latlong_grid)
  cover_frames<-lapply(cover_grids,
                       function(cover){
                         frame<-as(extent(cover$long0,cover$long1,cover$lat1,cover$lat0), "SpatialPolygons")
                         proj4string(frame) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
                         return(frame)
                       }
  )
  xmin<-xmin(frame_etow)
  xmax<-xmax(frame_etow)
  ymin<-ymin(frame_etow)
  ymax<-ymax(frame_etow)
  xd<-xmax-xmin
  yd<-ymax-ymin
  plot(frame_etow,col="blue",xlim=c(xmin-xd,xmax+xd),ylim=c(ymin-yd,ymax+yd))
  print(cover)
  print(cover_frames)
  print(frame_etow)
  t<-lapply(cover_frames,function(x){plot(x,add=TRUE)})
}
#
efgs_code <- "1kmN3211E4322"
Plot_cover_and_efgs_code(efgs_code)

第５章　世界メッシュ統計

Rソースコード 5.1

世界メッシュコード関連関数を使う. (Rソースコード 5.1)


source("https://www.fttsus.jp/worldmesh/R/worldmesh.R")
latitude <- 34.5
longitude <- 135.5
mesh3 <- cal_meshcode3(latitude,longitude)
res <- meshcode_to_latlong_grid(mesh3)
cat(sprintf("%s: NW(%f,%f), SE(%f,%f)\n",mesh3,res$lat0,res$long0,res$lat1,res$long1))

Rソースコード 5.2

GADMで提供される行政界のポリゴンデータから行政界3次メッシュデータを作成する例. (Rソースコード 5.2)


library(sp)
library(maptools)
library(rgeos)
source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
#
confirmation <- function(V1,V2){
  Q <- gIntersection(V1,V2)
  if(is.null(Q)){
    return(F)
  }
  if(gArea(Q)>0){
    return(T)
  } else {
    return(F)
  }    
}
#
prefname <- "Kyoto"
mydir <- getwd()
ifilex <- paste0(mydir,"/gadm36_JPN_shp/gadm36_JPN_2.shp")
ofile <- paste0(mydir,"/JPN_",prefname,"_mesh3.csv")
a <- readShapePoly(ifilex)
city_index<-which(a@data$NAME_1==prefname)
if(length(city_index)==0) q()
header<-sprintf("#ISO_COUTRY,country,name1,name2,meshcode,lat0,long0,lat1,long1\n")
cat(file=ofile,header,append=F)
for(i in city_index){
  for(kk in 1:length(a@polygons[[i]]@Polygons)){
      pol <- a@polygons[[i]]@Polygons[[kk]]@coords
      polym <- sprintf("POLYGON((%f %f",pol[1,1],pol[1,2])
      for(jj in 2:nrow(pol)){
         polym <- sprintf("%s, %f %f",polym,pol[jj,1],pol[jj,2])
      }
      polym <- sprintf("%s))",polym)
      V2 <<- readWKT(polym)
      minlat <- min(pol[,2])
      maxlat <- max(pol[,2])
      minlong <- min(pol[,1])
      maxlong <- max(pol[,1])
      deltalat=30
      deltalong=45
      nlat <- ceiling((maxlat-minlat)*60*60/deltalat) 
      nlong <- ceiling((maxlong-minlong)*60*60/deltalong) 
      meshlist<-c()
      for(j1 in 0:(nlat+1)){
        lat <- minlat + j1*deltalat/60/60
        for(j2 in 0:(nlong+1)){
          long <- minlong + j2*deltalong/60/60
          meshcode<-cal_meshcode3(lat,long)
          res<-meshcode_to_latlong_grid(meshcode)
          mesh <- sprintf("POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))",min(res$long0,res$long1),min(res$lat0,res$lat1),max(res$long0,res$long1),min(res$lat0,res$lat1),max(res$long0,res$long1),max(res$lat0,res$lat1),min(res$long0,res$long1),max(res$lat0,res$lat1),min(res$long0,res$long1),min(res$lat0,res$lat1))
          V1 <- readWKT(mesh)
          if(confirmation(V1,V2)){
             st<-sprintf("%s,%s,%s,%s",as.character(a@data[i,]$GID_0)[1],as.character(a@data[i,]$NAME_0)[1],as.character(a@data[i,]$NAME_1)[1],as.character(a@data[i,]$NAME_2)[1])
             st<-sprintf("%s,%s,%f,%f,%f,%f\n",st,meshcode,res$lat0,res$long0,res$lat1,res$long1)
	     cat(st)
             cat(file=ofile,st,append=T)
          }
       }
    }
  }
}

Rソースコード 5.3

オーストラリア標準グリッドから世界メッシュへの変換(1). Rソースコード5.4と共に用いる. (Rソースコード 5.3)


library(sp)
library(rgeos)
library(raster)
source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
#
Gen_frame_a_to_w<-function(aea_x,aea_y){
  X <- as.numeric(aea_x)
  Y <- as.numeric(aea_y)
  frame <- as(extent(X,X+1000,Y,Y+1000), "SpatialPolygons")
  proj4string(frame) <- CRS("+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
  frame_atow <- spTransform(frame,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  return(frame_atow)
}
Cal_cover_aea_to_world<-function(aea_x,aea_y){
  frame_atow<-Gen_frame_a_to_w(aea_x,aea_y)
  min_long<-xmin(frame_atow)
  min_lat<-ymin(frame_atow)
  max_long<-xmax(frame_atow)
  max_lat<-ymax(frame_atow)
  mid_lat<-(min_lat+max_lat)/2
  mid_long<-(min_long+max_long)/2
  NE=c(max_lat,max_long)
  NW=c(max_lat,min_long)
  SE=c(min_lat,max_long)
  SW=c(min_lat,min_long)
  M1=c(min_lat,mid_long)
  M2=c(max_lat,mid_long)
  M3=c(mid_lat,min_long)
  M4=c(mid_lat,max_long)
  check_points<-rbind(NE,NW,SE,SW,M1,M2,M3,M4)
  cover<-apply(check_points,1,function(x){return(cal_meshcode3(x[1],x[2]))})
  return(unique(cover))
}
Cal_rate_aea_to_world<-function(aea_x,aea_y,total_pop){
  cover<-Cal_cover_aea_to_world(aea_x,aea_y)
  frame_atow<-Gen_frame_a_to_w(aea_x,aea_y)
  cover_grids<-lapply(cover,meshcode_to_latlong_grid)
  cover_frames<-lapply(cover_grids,
                       function(cover){
                         frame<-as(extent(cover$long0,cover$long1,cover$lat1,cover$lat0), "SpatialPolygons")
                         proj4string(frame) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
                         return(frame)
                       }
  )
  rates<-c()
  for(i in 1:length(cover)){
    if(gIntersects(cover_frames[[i]],frame_atow)){
      #have common part
#      print(gArea(gIntersection(cover_frames[[i]],frame_atow)))      
      r<-gArea(gIntersection(cover_frames[[i]],frame_atow))/gArea(frame_atow)
    }
    else{
      #not have common part
      r<-0
    }
    rates<-c(rates,r)
  }
  total_pop_w<-total_pop*rates
  result<-data.frame(meshcode=cover,TOT_P=total_pop_w)
  return(result)
}

Rソースコード 5.4

オーストラリア標準グリッドから世界メッシュへの変換(2). Rソースコード5.3と共に用いる. (Rソースコード 5.4)

outdir<-getwd()
infile<-paste0(outdir,"/Australian_Population_Grid_2011.tif")
ofile1<-paste0(outdir,"/POP_2011_mid_AUS.tsv")
ofile2<-paste0(outdir,"/POP_2011_AUS.tsv")
if(!file.exists(ofile2)){
  sat <- raster(infile)
  xmin <- xmin(sat)
  ymax <- ymax(sat)
  cols <- ncol(sat)
  rows <- nrow(sat)
  dx <- (xmax(sat)-xmin(sat))/cols
  dy <- (ymax(sat)-ymin(sat))/rows
  kk <- 0
  val<-as.matrix(sat)
  for(x in 1:cols){
    tot_pop = val[1:rows,x]
    aea_x <- xmin + dx*(x-1)
    print(x)
    for(y in 1:rows){
      if(tot_pop[y]!=0){
        aea_y <- ymax - dy*(y-1)
        print(aea_x)
        print(aea_y)
        print(tot_pop[y])
        result<-Cal_rate_aea_to_world(aea_x,aea_y,tot_pop[y])
        if(kk == 0){
          write.table(result,file=ofile1,sep="\t",append=FALSE,col.names=TRUE,row.names=FALSE,quote=FALSE)
      	  kk <- kk + 1
        }else{
          write.table(result,file=ofile1,sep="\t",append=TRUE,col.names=FALSE,row.names=FALSE,quote=FALSE)
        }
      }
    }
  }
  #
  TOT_P_world<-read.csv(file=ofile1,header=TRUE,sep="\t",colClasses=c("character","numeric"))
  #
  TOT_P_sum<-aggregate(TOT_P_world$TOT_P,by=list(TOT_P_world$meshcode),FUN=sum)
  colnames(TOT_P_sum)<-c("#meshcode","TOT_P")
  #
  latlongs<-lapply(TOT_P_sum$'#meshcode',meshcode_to_latlong_grid) 
  TOT_P_sum<-transform(TOT_P_sum,lat0=sapply(latlongs,function(x){x$lat0})) 
  TOT_P_sum<-transform(TOT_P_sum,long0=sapply(latlongs,function(x){x$long0})) 
  TOT_P_sum<-transform(TOT_P_sum,lat1=sapply(latlongs,function(x){x$lat1}))
  TOT_P_sum<-transform(TOT_P_sum,long1=sapply(latlongs,function(x){x$long1})) 
  # 
  write.table(TOT_P_sum,file=ofile2,sep="\t",col.names = TRUE,row.names = FALSE,append=FALSE,quote=FALSE) 
}

Rソースコード 5.5

JAXA AW3D30から標高に関する3次世界メッシュ統計を作成するためのソースコード. (Rソースコード 5.5)


library(tiff)
source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
myextract<-function(t,targetdir){
  tt<-unlist(strsplit(t,"[/.]"))
  ff<-unlist(strsplit(t,"[_.]"))
  WBDfile<-gsub("DSM","MSK",t)
  cmn<-length(tt)
  out<-paste(targetdir,"/",tt[cmn-2],"/",tt[cmn-1],".tsv",sep="")
  name<-unlist(strsplit(t,"[/_.]"))
  name<-name[length(name)-3]
  name_c<-unlist(strsplit(name,""))
  if(name_c[1]=="N" && name_c[5]=="E"){
      latlong<-unlist(strsplit(name,"['N','E']"))
      lat0<-as.integer(latlong[2]) # Southern-Western location
      long0<-as.integer(latlong[3])
  }
  if(name_c[1]=="N" && name_c[5]=="W"){
      latlong<-unlist(strsplit(name,"['N','W']"))
      lat0<-as.integer(latlong[2]) # Southern-Western location
      long0<-as.integer(latlong[3])
  }
  if(name_c[1]=="S" && name_c[5]=="E"){
      latlong<-unlist(strsplit(name,"['S','E']"))
      lat0<-as.integer(latlong[2])
      long0<-as.integer(latlong[3])
  }
  if(name_c[1]=="S" && name_c[5]=="W"){
      latlong<-unlist(strsplit(name,"['S','W']"))
      lat0<-as.integer(latlong[2])
      long0<-as.integer(latlong[3])
  }
  lat0<-lat0+1
  long0<-long0
  header<-sprintf("#meshcode\talt_min\talt_mean\talt_median\talt_max\tlat0\tlong0\tlat1\tlong1\n");
  cat(file=out,header,append=F)
  dd<-data.frame(meshcode=c(),altitude=c())
  sat <- readTIFF(t,as.is=TRUE)
  wbd <- readTIFF(WBDfile,as.is=TRUE)
  for(yy in 1:120){
    lat <- lat0 - (yy-1)*30/3600 - 30/2/3600
    O <- c()
    for(xx in 1:80){
      long<-long0 + (xx-1)*45/3600 + 45/2/3600
      alt<-sat[((yy-1)*30+1):(yy*30),((xx-1)*45+1):(xx*45)]
      att<-wbd[((yy-1)*30+1):(yy*30),((xx-1)*45+1):(xx*45)]
      alt<-alt[att==0]
      O<-alt[alt!=-9999]
      O[O>32768]<-O[O>32768]-65536
      if(length(O)!=0) {
        alt_min<-min(O)
        alt_mean<-mean(O)
        alt_median<-median(O)
        alt_max<-max(O)
        meshcode3<-cal_meshcode3(lat,long)
        res<-meshcode_to_latlong_grid(meshcode3)
        st<-sprintf("%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",meshcode3,alt_min,alt_mean,alt_median,alt_max,res$lat0,res$long0,res$lat1,res$long1)
        cat(file=out,st,append=T)
        cat(st)
      }
    }
  }
}
ls<-dir(pattern="_AVE_DSM.tif")
for(f in ls){
  myextract(f,".")
}

Rソースコード 5.6

ポリゴン領域内における世界メッシュ統計の再集計. (Rソースコード 5.6)


library(sp)
library(maptools)
library(rgeos)
library(rgdal)
source("http://www.fttsus.jp/worldmesh/R/worldmesh.R")
contribution <- function(V1,V2){
  Q <- gIntersection(V1,V2)
  if(is.null(Q)){
    return(0.0)
  }
  if(gArea(Q)>0){
    return(gArea(Q)/gArea(V1))
  } else {
    return(0.0)
  }    
}
wgs.84<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
aea<-"+proj=aea +lat_1=18 +lat_2=52 +lat_0=0 +lon_0=135 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
prefname <- "Kyoto"
mydir <- getwd()
ifilex <- paste0(mydir,"/gadm36_JPN_shp/gadm36_JPN_2.shp")
ifile2 <- paste0(mydir,"/NASA_JPN_mesh3.tsv")
summary <- paste0(mydir,"/summary.csv")
cat(file=summary,sprintf("#prefname,cityname,light_int,area,light_int_per_area\n"),append=F)
a <- readOGR(ifilex)
nasa <- read.csv(file=ifile2,sep="\t",header=T)
city_index<-which(a@data$NAME_1==prefname)
if(length(city_index)==0) q()
for(i in city_index){
  area_sum = 0.0  
  light_int_sum = 0.0
  for(kk in 1:length(a@polygons[[i]]@Polygons)){
      pol <- a@polygons[[i]]@Polygons[[kk]]@coords
      polym <- sprintf("POLYGON((%f %f",pol[1,1],pol[1,2])
      for(jj in 2:nrow(pol)){
         polym <- sprintf("%s, %f %f",polym,pol[jj,1],pol[jj,2])
      }
      polym <- sprintf("%s))",polym)
      V2 <- readWKT(polym,p4s=CRS(wgs.84))
      V2 <- spTransform(V2, CRS(aea))
      minlat <- min(pol[,2])
      maxlat <- max(pol[,2])
      minlong <- min(pol[,1])
      maxlong <- max(pol[,1])
      deltalat=30
      deltalong=45
      nlat <- ceiling((maxlat-minlat)*60*60/deltalat) 
      nlong <- ceiling((maxlong-minlong)*60*60/deltalong) 
      meshlist<-c()
      for(j1 in 0:(nlat+1)){
        lat <- minlat + j1*deltalat/60/60
        for(j2 in 0:(nlong+1)){
          long <- minlong + j2*deltalong/60/60
          meshcode<-cal_meshcode3(lat,long)
          res<-meshcode_to_latlong_grid(meshcode)
          mesh <- sprintf("POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))",min(res$long0,res$long1),min(res$lat0,res$lat1),max(res$long0,res$long1),min(res$lat0,res$lat1),max(res$long0,res$long1),max(res$lat0,res$lat1),min(res$long0,res$long1),max(res$lat0,res$lat1),min(res$long0,res$long1),min(res$lat0,res$lat1))
          V1 <- readWKT(mesh,p4s=CRS(wgs.84))
          V1 <- spTransform(V1,CRS(aea))
	  rho = contribution(V1,V2)
          lval = nasa[nasa$X.meshcode==meshcode,]$light_int*rho
          aval = gArea(V1)*rho
          light_int_sum = light_int_sum + lval
          area_sum = area_sum + aval*1e-06
       }
    }
  }
  cat(file=summary,sprintf("%s,%s,%f,%f,%f\n",as.character(a@data[i,]$NAME_1)[1],as.character(a@data[i,]$NAME_2)[1],light_int_sum,area_sum,light_int_sum/area_sum),append=T)
}

第６章　世界メッシュ統計の分析例とワークフロー

Rソースコード 6.1

世界メッシュ統計の可視化. (Rソースコード6.1)


library(leaflet)
library(mapview)
cols<-colorRamp(c("#000000","white","#faf500")) #color
pal<-colorNumeric(palette=cols, domain=(0:255))
homedir <- getwd() # set homedir
targetdir<-dir(path=homedir,pattern="_mesh3.tsv")
for(d in targetdir){
  tt<-unlist(strsplit(d,"[/.]"))
  if(tt[2]=="tsv"){
    infile<-d
    ofile<-paste0(homedir,"/",tt[1],".png")
    t<-read.csv(paste0(homedir,"/",infile), sep="\t", header=TRUE, encoding='UTF-8')
    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){
        median_t<-median(t[(t$lat0 < latmin+dlat) & (t$lat0 > latmin) & (t$long0 > longmin) & (t$long0 < longmin+dlong),]$light_int)
        if(!is.na(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=ofile)
  }
}