LINK:[GoogleEarthとR][ShapeFileライブラリ][kmlラボ][空間的なデータの分析][Rでジオコーディング][RでGPS]


目次


目的

rgdalパッケージのwriteOGR関数を使うと、ShapeFileからKMLファイルを作成することができますが、単にKMLファイルを出力するだけで、細かい設定ができず、あまり利用価値がありません。
そこでXMLパッケージを用いて、作成されたKMLファイルを操作することで、KMLファイルの利用価値を高めることができると考えました。(ついでにXMLパッケージをつかいこなすぞ!)
注)ここは、「KMLの公開」を目的としたサイトではありません。公開されているKMLファイルは、あくまでもRで加工するための基礎データにすぎません。

対応表

RとXMLパッケージのバージョンによってはうまく動かないので、注意してください。
推奨環境:Windows: R2.9.2 - XML1.96

  • Windows: R2.9.2 - XML1.9.6 O
  • Windows: R2.7.2 - XML1.94 O (CDATA内部の日本語文字化けなし)
    install.packages("XML", repos = "http://www.omegahat.org/R")
  • Windows: R2.7.2 - XML1.96 X (CDATA内部の日本語文字化け)
  • Windows: R2.6.2 - XML1.94 O
  • MacOSX10.5.5: R2.7.2 O

関連リンク

すばらしいサイトです。ここで勉強させてもらってます。
TME
EarthAtlas

kmlライブラリ

公開データの著作権等について:
*「地域境界KMLFileデータの作成に当たっては、ESRIジャパン株式会社の全国市区町村界データを使用しました。本データの著作権はESRIジャパン株式会社に帰属します。なお、地域境界KMLFileデータは営利目的には利用できません。」
*「地域メッシュKMLFileデータの作成にあたっては、財団法人 地域地盤環境研究所 村上貴志氏のメッシュデータを使用しました。本データの著作権は村上貴志氏に帰属します。また、地域メッシュKMLFileデータの利用に関してですが,著作者の表記(書式自由)をして頂ければ自由にご利用ください。公開されたデータが非営利・営利いずれの使われ方をしても構いません。」

県コード地域名称地域境界KMLFile地域メッシュ名称地域メッシュKMLFile備考
#日本全図 第1次地域区画
#日本全図第2次地域区画
01北海道 1kmメッシュ
02青森県 1kmメッシュ
03岩手県 1kmメッシュ
04宮城県 1kmメッシュ
05秋田県 1kmメッシュ
10群馬県

国勢調査のデータのダウンロード

http://www.e-stat.go.jp/SG2/toukeichiri/TopFrame.do?fromPage=init&toPage=download

コメント欄


  • XML対応表更新 -- okinawa 2009-11-30 (月) 18:21:41
  • maptoolsにkmlLine,kmlPolygon関数が追加されました。これらの関数にて生成されるkmlはwriteOGRで生成されるkmlと構造が違うためkmlFunctions.rはそのままでは使えません。 -- okinawa 2009-02-20 (金) 11:26:41
  • さて、そろそろ2次メッシュにいこか・・・。 -- okinawa 2009-02-17 (火) 17:31:12
  • メッシュkmlデータの作成:1次メッシュ生成関数作成終了。 -- okinawa 2009-02-13 (金) 16:52:40
  • メッシュkmlデータの作成更新。メッシュを連続生成できるようにしました。まだ作りかけですが。 -- okinawa 2009-02-13 (金) 09:07:23
  • メッシュkmlデータの作成を追加。現在、初期段階のコードです。 -- okinawa 2009-02-12 (木) 14:50:36
  • GE5向けにkmlFunctions.rを変更しました。今後、ここの内容はGE5向けに変更していきますのでご了承ください。 -- okinawa 2009-02-05 (木) 09:33:38
  • データを読み込んでkmlを再構成する(2) 那覇市国勢調査編のコードを若干変更(ポリゴンホールを回避します) -- okinawa 2009-02-03 (火) 21:59:56
  • 今後、TMEの技術を取り入れていこうと思います。・・・まだ自動でメッシュデータ作るのが残ってますが・・・。 -- okinawa 2008-12-19 (金) 09:22:06
  • ここの名称を「kmlラボ」に変更しました。 -- okinawa 2008-12-08 (月) 08:41:24
  • kmlFunctions.rバグ修正。 -- okinawa 2008-11-26 (水) 18:25:28
  • useInternalNodes? = TRUEを使わんといかんのか・・・。utf-8が標準じゃないwindowsの方がめんどくさいですね。 -- okinawa 2008-11-18 (火) 09:37:57
  • 那覇市の周辺市町村(浦添市、豊見城市、南風原町)を同時に表示してみた。 -- okinawa 2008-11-17 (月) 17:15:05
  • getCDATA関数復活。これでPlaceMark?名にCDATA内の指定データをnameタグに挿入できます。これをsetAltitude()と一緒に使うことで、placemarkに日本語が表示できるようになると思います。 (Macでは成功、Windowsではうまくいかないなあ・・・)-- okinawa 2008-11-16 (日) 13:50:16
  • MacOSXでrgdalを使えるようにするを追加。これでMacでもKML自由自在ですね。 -- okinawa 2008-11-15 (土) 15:26:01
  • 明日、那覇市長選がある。那覇市国勢調査編として分析したが、3Dグラフから那覇市はすでに人口が空洞化していると考えられる。このような地域は企業が集中している反面、取り残された住民に高齢者の割合が高くなる。いわゆる都市構造的な問題は、一自治体の政策だけでは解決できないと考えられるが、いかがだろうか? -- okinawa 2008-11-15 (土) 12:41:58
  • kmlFunctions.rのsetAltitude()にラベルのパラメータ指定を追加。ラベルを指定することで、KML内のPlaceMark?名をセットする。ただし、現状では日本語をセットするとエラーになる。XMLパッケージのロケール問題っぽい。 -- okinawa 2008-11-14 (金) 14:53:02
  • kmlFunctions.rバグ取り&ちょっと修正。 -- okinawa 2008-11-12 (水) 15:30:15
  • データを読み込んでkmlを再構成する(2) 那覇市国勢調査編追加 -- okinawa 2008-11-07 (金) 14:49:09
  • kmlAltitudeNode?関数修正、setAltitudeMesh?関数をsetAltitude関数に統合。MGパラメータでメッシュKMLと切り替えます。 -- okinawa 2008-11-07 (金) 11:56:58
  • 次はメッシュデータを自動でつくってみます。 -- okinawa 2008-11-06 (木) 15:20:27
  • データを読み込んで、kmlを再構成するを追加。shapefileから直接つくるのにくらべ簡素に?なりました。 -- okinawa 2008-11-06 (木) 14:50:28
  • kmlFunctions.r追加。 -- okinawa 2008-11-06 (木) 11:06:04
  • colors.kml関数追加しました。 -- okinawa 2008-11-05 (水) 16:19:15
  • getKmlNode?関数追加。ノードを指定してKMLの切り抜きができます。 -- okinawa 2008-10-24 (金) 14:47:32
  • 関数追加。あとは色の変換関数がいるな。 -- okinawa 2008-10-23 (木) 21:50:06
  • 関数を若干修正 -- okinawa 2008-10-23 (木) 08:56:03
  • 地域メッシュKMLデータの応用追加 -- okinawa 2008-10-22 (水) 18:22:05
  • kmlファイルをRから直接いじってみました。今後は、ライブラリを充実させ、色々な機能追加をしてみようと思います。 -- okinawa 2008-10-21 (火) 18:03:29

kmlデータをXMLパッケージで再構成する関数

writeOGRで出力されたkmlデータをXMLパッケージを用いて再構成し、高さ、ポリゴンの色、ポリゴンの線の色を指定できる関数を作成しました。
(構造決めうちなので、ShapeFileライブラリのデータをwriteOGRでkml出力したデータしか使えないと思います)

kmlFunctions.r

kmlを操作する関数をまとめました。下記のように、適当な場所に保存して、souceで呼び出してください。(GE5向けに新しくしました)

古いものは↓

library(XML) #宣言
source("c:/kmlFunctions.r")#ソースの読み込み
doc <- xmlTreeParse("c:/gunma.kml")#土台となるkmlファイル読み込み
data<-sample(1000:10000,39) #高度データを自動生成
fillcolor<-colors.kml(heat.colors(39),"EE")#色の割り当て(ポリゴン)
linecolor<-NA#色の割り当て(線)
x<-kmlAltitudeNode(doc,data,fillcolor,linecolor)#kmlファイルにデータを書き込み
saveXML(x,"c:/gunma2.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')#ファイルに出力


color2.JPG

kmlAltitudeNode?関数(GE5)

kmlAltitudeNode<-function(doc,data,fillcolor,linecolor,extrude="1",tessellate="1",altitudeMode="absolute",MG=T){
kmlnode <- xmlRoot(doc)
node0<-kmlnode[[1]]
kmlver5<-'http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2'
kmlnode<-xmlNode("kml",node0,attrs=c(xmlns=kmlver5))
len1<-length(node0[[1]])
for(i in 2:len1){
if(MG==T){
node1<-node0[[1]][[i]][[2]][[1]]
}else if(MG==F){
node1<-node0[[1]][[i]][[2]]
}
node1<-addChildren(node1,xmlNode("extrude",extrude),xmlNode("tessellate",tessellate),xmlNode("gx:altitudeMode",altitudeMode))
node2<-node0[[1]][[i]][[3]][[2]]
node2<-xmlNode("PolyStyle",xmlNode("color",fillcolor[i-1]))
#linecolor
node3<-node0[[1]][[i]][[3]][[1]]
node3<-xmlNode("LineStyle",xmlNode("color",linecolor))
#polygon
len2<-length(node0[[1]][[i]][[2]][[1]])
for(j in 1:len2){
   x1<-xmlValue(node1[[j]][[1]][[1]])
   x2<-rev(unlist(strsplit(x1," ")))
   x3<-paste(x2,data[[i-1]],sep=",",collapse=" ")
   node1[[j]][[1]][[1]]<-xmlNode("coordinates",x3)
}
#replace node
if(MG==T){
node0[[1]][[i]][[2]][[1]]<-node1
}else if(MG==F){
node0[[1]][[i]][[2]]<-node1
}
node0[[1]][[i]][[3]][[2]]<-node2
node0[[1]][[i]][[3]][[1]]<-node3
kmlnode[[1]]<-node0
}
return(kmlnode)
}
#sample
library(XML)
doc <- xmlTreeParse("c:/gunma.kml")
data<-sample(1000:10000,39)
fillcolor<-c("ff000000","ffff0000","ff00ff00","ff0000ff","ffffff00","ffff00ff","ff00ffff","ffffffff","00000000")
linecolor<-"NA"
x<-kmlAltitudeNode(doc,data,fillcolor,linecolor)
saveXML(x,"c:/gunma2.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
xml1.JPG

setAltitude関数(GE5)

ノード番号を指定してピンポイントでAltitudeをセットする。

setAltitude<-function(kmlnode,nodeno,data,fillcolor,linecolor,extrude="1",tessellate="1",altitudeMode="absolute",MG=T,label=""){
#kmlnode <- xmlRoot(doc)
node0<-kmlnode[[1]]
kmlver5<-'http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2'
kmlnode<-xmlNode("kml",node0,attrs=c(xmlns=kmlver5))
if(MG==T){
#altitude
node1<-node0[[1]][[nodeno]][[2]][[1]]
node1<-addChildren(node1,xmlNode("extrude",extrude),xmlNode("tessellate",tessellate),xmlNode("gx:altitudeMode",altitudeMode))
#fillcolor
node2<-node0[[1]][[nodeno]][[3]][[2]]
node2<-xmlNode("PolyStyle",xmlNode("color",fillcolor))
#linecolor
node3<-node0[[1]][[nodeno]][[3]][[1]]
node3<-xmlNode("LineStyle",xmlNode("color",linecolor))
#label
node4<-node0[[1]][[nodeno]]
node4<-addChildren(node4,xmlNode("name",label))
#polygon
len2<-length(node0[[1]][[nodeno]][[2]][[1]])
j<-1:len2
text1<-paste('x1<-xmlValue(node1[[',j,']][[1]][[1]]);x2<-rev(unlist(strsplit(x1," ")));',sep='')
text2<-paste('x3<-paste(x2,data,sep=",",collapse=" ");node1[[',j,']][[1]][[1]]<-xmlNode("coordinates",x3)',sep='')
text<-paste(text1,text2,sep="")
eval(parse(text=text))
node0[[1]][[nodeno]]<-node4
node0[[1]][[nodeno]][[2]][[1]]<-node1
node0[[1]][[nodeno]][[3]][[2]]<-node2
node0[[1]][[nodeno]][[3]][[1]]<-node3
kmlnode[[1]]<-node0
}else if(MG==F){
#altitude
node1<-node0[[1]][[nodeno]][[2]]
node1<-addChildren(node1,xmlNode("extrude",extrude),xmlNode("tessellate",tessellate),xmlNode("gx:altitudeMode",altitudeMode))
#fillcolor
node2<-node0[[1]][[nodeno]][[3]][[2]]
node2<-xmlNode("PolyStyle",xmlNode("color",fillcolor))
#linecolor
node3<-node0[[1]][[nodeno]][[3]][[1]]
node3<-xmlNode("LineStyle",xmlNode("color",linecolor))
#label
node4<-node0[[1]][[nodeno]]
node4<-addChildren(node4,xmlNode("name",label))
#polygon
len2<-length(node0[[1]][[nodeno]][[2]][[1]])
j<-1:len2
text1<-paste('x1<-xmlValue(node1[[',j,']][[1]][[1]]);x2<-rev(unlist(strsplit(x1," ")));',sep='')
text2<-paste('x3<-paste(x2,data,sep=",",collapse=" ");node1[[',j,']][[1]][[1]]<-xmlNode("coordinates",x3)',sep='')
text<-paste(text1,text2,sep="")
eval(parse(text=text))
node0[[1]][[nodeno]]<-node4
node0[[1]][[nodeno]][[2]]<-node1
node0[[1]][[nodeno]][[3]][[2]]<-node2
node0[[1]][[nodeno]][[3]][[1]]<-node3
kmlnode[[1]]<-node0
}
return(kmlnode)
}

getCDATANo関数と組み合わせて必要な場所のみ設定

#sample
doc <- xmlTreeParse("c:/gunma.kml")
x<-xmlRoot(doc)
#1つめの地域
name<-"JCODE"
value<-"10366"
nodeno<-getCDATANo(x,name,value)
data<-1000
fillcolor<-c("ffff0000")
linecolor<-"ff000000"
x<-setAltitude(x,nodeno,data,fillcolor,linecolor)
#2つめの地域
name<-"JCODE"
value<-"10201"
nodeno<-getCDATANo(x,name,value)
data<-1000
fillcolor<-c("ffff00ff")
linecolor<-"ff000000"
x<-setAltitude(x,nodeno,data,fillcolor,linecolor)
#保存
saveXML(x,"c:/gunma2.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
xml4.JPG

メッシュのKMLはMultiGeometory?ノードが無いので、MGパラメータでF指定する

#sample
library(XML)
doc <- xmlTreeParse("c:/StdMesh01.kml")
kmlnode<-xmlRoot(doc)
#1
name<-"CODE"
value<-"3622"
nodeno<-getCDATANo(kmlnode,name,value)
data<-1000000
fillcolor<-c("ffff0000")
linecolor<-"ff000000"
kmlnode<-setAltitude(kmlnode,nodeno,data,fillcolor,linecolor,MG=F)
#2
name<-"CODE"
value<-"3623"
nodeno<-getCDATANo(kmlnode,name,value)
data<-800000
fillcolor<-c("ffff00ff")
linecolor<-"ff000000"
kmlnode<-setAltitude(kmlnode,nodeno,data,fillcolor,linecolor,MG=F)
#save
saveXML(kmlnode,"c:/StdMesh012.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
color4.JPG

kmlNodeLength?関数

kmlAltitudeNode?関数を使うときに、kmlのノード数を知りたい時に使う。

#kmlNodeLength() ver 0.1.2
kmlNodeLength<-function(doc){
library(XML)
kmlnode <- xmlRoot(doc)
node0<-kmlnode[[1]]
len1<-length(node0[[1]])
return(len1)
}
#sample
library(XML)
doc <- xmlTreeParse("c:/gunma.kml")
kmlNodeLength(doc)
> #sample
> library(XML)
> doc <- xmlTreeParse("c:/gunma.kml")
> kmlNodeLength(doc)
[1] 40

注)先頭のノードに<name>gunma</name>が入っていますので、市町村数はノード数−1となります。

getCDATAList関数

CDATA(HTML)の中身を取り出す。nameに対応する項目名を指定する。(必須)

getCDATAList<-function(doc,name){
getCDATA<-function(node,name){
x<-gsub(" ","",iconv(xmlValue(node),"utf-8","cp932"))
x<-unlist(strsplit(x,"<br/>\n"))
x<-x[[grep(name,x)]]
x<-gsub("<b>","",gsub("</b>","",gsub("<i>","",gsub("</i>","",gsub("<br/>","",gsub(paste(name,":",sep=""),"",x))))))
#x<-iconv(x,"","utf-8")
return(x)
}
kmlnode <- xmlRoot(doc)
node0<-kmlnode[[1]]
len1<-length(node0[[1]])
x<-NULL
 for(i in 2:len1){
 x[i-1]<-getCDATA(node0[[1]][[i]][[1]][[1]],name)
 }
 return(x)
}
#sample
doc <- xmlTreeParse("c:/gunma.kml")
name<-"CITY1"
getCDATAList(doc,name)
> #sample
> doc <- xmlTreeParse("c:/gunma.kml")
> name<-"CITY1"
> getCDATAList(doc,name)
 [1] "前橋市"     "高崎市"     "桐生市"     "伊勢崎市"   "太田市"     "沼田市"     "館林市"     "渋川市"     "藤岡市"     "富岡市"    
[11] "安中市"     "みどり市"   "富士見村"   "榛名町"     "榛東村"     "吉岡町"     "吉井町"     "上野村"     "神流町"     "下仁田町"  
[21] "南牧村"     "甘楽町"     "中之条町"   "長野原町"   "嬬恋村"     "草津町"     "六合村"     "高山村"     "東吾妻町"   "片品村"    
[31] "川場村"     "昭和村"     "みなかみ町" "玉村町"     "板倉町"     "明和町"     "千代田町"   "大泉町"     "邑楽町"    
> #sample
> doc <- xmlTreeParse("c:/gunma.kml")
> name<-"JCODE"
> getCDATAList(doc,name)
[1] "10201" "10202" "10203" "10204" "10205" "10206" "10207" "10208" "10209" "10210" "10211" "10212" "10303" "10321" "10344" "10345" "10363"
[18] "10366" "10367" "10382" "10383" "10384" "10421" "10424" "10425" "10426" "10427" "10428" "10429" "10443" "10444" "10448" "10449" "10464"
[35] "10521" "10522" "10523" "10524" "10525"

getCDATA関数

CDATAのデータを取り出す。cdatanoは、grepで検索した場合、複数の位置が取り出された場合、初期位置からの相対位置を指定するパラメータ。

getCDATA<-function(kmlnode,name,nodeno=1,sep=":",cdatano=1){
node0<-kmlnode[[1]]
node1<-node0[[1]][[nodeno]][[1]][[1]]
#x<-gsub(" ","",xmlValue(node1))#mac
x<-gsub(" ","",iconv(xmlValue(node1),"utf-8","cp932"))
#print(x)
x<-unlist(strsplit(x,"<br/>\n"))
xno<-grep(name,x)
len<-length(xno)
if(len>1){
x<-x[[xno[cdatano]]]
}else if(len==1){
x<-x[[xno]]
}
#print(x)
x<-gsub("<b>","",gsub("</b>","",gsub("<i>","",gsub("</i>","",gsub("<br/>","",x)))))
x<-gsub(paste(name,sep,sep=""),"",x)
return(x)
}
> #sample
> doc <- xmlTreeParse("/gunma.kml")
> kmlnode<-xmlRoot(doc)
> name<-"JCODE"
> nodeno<-2
> getCDATA(kmlnode,name,nodeno)
[1] "10201" 

getCDATANo関数

nameとvalueを指定することで、CDATA(HTML)内のデータを検索し、一致したノード番号を返す。

#getCDATANo
getCDATANo<-function(kmlnode,name,value){
getCDATA2<-function(node,name){
x<-gsub(" ","",iconv(xmlValue(node),"utf-8","cp932"))
x<-unlist(strsplit(x,"<br/>\n"))
x<-x[[grep(name,x)]]
x<-gsub("<b>","",gsub("</b>","",gsub("<i>","",gsub("</i>","",gsub("<br/>","",x)))))
#x<-iconv(x,"","utf-8")
return(x)
}
#kmlnode <- xmlRoot(doc)
namevalue<-paste(name,":",value,sep="")
node0<-kmlnode[[1]]
len1<-length(node0[[1]])
for(i in 2:len1){
   if(namevalue==getCDATA2(node0[[1]][[i]][[1]][[1]],name)){
   x<-i
   break
   }else{
   x<-NA
   }
 }
 return(x)
}
#sample
doc <- xmlTreeParse("c:/gunma.kml")
kmlnode<-xmlRoot(doc)
name<-"JCODE"
value<-"10366"
getCDATANo(kmlnode,name,value)
> #sample
> doc <- xmlTreeParse("c:/gunma.kml")
> name<-"JCODE"
> value<-"10366"
> getCDATANo(kmlnode,name,value)
[1] 19

getKmlNode?関数(GE5)

ノードをnodenoで指定し、指定したノードの部分だけ抽出してKMLを再構成

#getKmlNode
getKmlNode<-function(kmlnode,nodeno){
len<-length(nodeno)
node0<-kmlnode[[1]]
kmlver5<-'http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2'
kmlnode<-xmlNode("kml",node0,attrs=c(xmlns=kmlver5))
node2<-xmlNode("Folder",xmlNode("name",as.character(nodeno)))
for(i in 1:len){
   node1<-node0[[1]][[nodeno[[i]]]]
   node2<-addChildren(node2,node1)
}
kmlnode[[1]][[1]]<-node2
return(kmlnode)
}
#sample
library(XML)
doc <- xmlTreeParse("c:/gunma.kml")
kmlnode<-xmlRoot(doc)
nodeno<-c(2,4,6)
expdata<-getKmlNode(kmlnode,nodeno)
saveXML(expdata,"c:/gunma3.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')

colors.kml関数

colors.kml("色名","アルファチャンネル")

colors.kml<-function(rgb=colors(),alpha="FF"){
         colors.hex <- function( x=colors() ) { 
         color.hex <- function(x) do.call( "rgb", as.list(col2rgb(x)/255) ) 
         sapply( x, color.hex ) 
         }
col.hex<-colors.hex(rgb)
red<-substr(col.hex,2,3)
green<-substr(col.hex,4,5)
blue<-substr(col.hex,6,7)
alpha<-alpha
col.kml<-paste(alpha,blue,green,red,sep="")
return(col.kml)
}
> colors.kml("red")
[1] "FF0000FF"
> colors.kml(rainbow(10))
 [1] "FF0000FF" "FF0099FF" "FF00FFCC" "FF00FF33" "FF66FF00" "FFFFFF00" "FFFF6600" "FFFF0033" "FFFF00CC"
[10] "FF9900FF"

参考:色見本

メッシュkmlデータの作成

メッシュkmlデータを自動作成します。生成されたkmlはもちろん他の関数で操作できます。

makeMesh1関数:第1次地域区画(1次メッシュ)生成関数

指定された地域メッシュコードを基に1次メッシュをブロック生成します。
下記のDefaultMesh1.kmlをタネにして作成しますので、ダウンロードして利用してください。

makeMesh1<-function(kmlnode,codea,codeb){
transMeshCode<-function(code){
lat0<-as.integer(substr(code,1,2))/1.5
lon0<-as.integer(substr(code,3,4))+100
lat1<-lat0;lon1<-lon0+1;lat2<-lat0+(2/3);lon2<-lon1
lat3<-lat2;lon3<-lon0
x<-paste(lon0,",",lat0," ",lon1,",",lat1," ",lon2,",",lat2," ",lon3,",",lat3," ",lon0,",",lat0,sep="")
return(x)
}
code1<-as.integer(substr(codea,1,2))
code3<-as.integer(substr(codea,3,4))
code2<-as.integer(substr(codeb,1,2))
code4<-as.integer(substr(codeb,3,4))
print(code1)
print(code2)
print(code3)
print(code4)
node0<-kmlnode[[1]]
node1<-node0[[1]]
node2<-node0[[1]][[2]]
x<<-node0[[1]]
nodex<-removeChildren(x,kids=c("Placemark"))
for(j in code1:code2){
code01<-paste(as.character(j),code3,sep="")
code02<-paste(as.character(j),code4,sep="")
print(j)
for(i in code01:code02){
print(i)
node3<-node0[[1]][[2]][[1]][[1]]
node4<-node0[[1]][[2]][[2]][[1]][[1]][[1]]
CDATA<-paste("<b>CODE:</b><i>",i,"</i><br />\n",sep="")
node3<-xmlCDataNode(CDATA)
node4<-xmlNode("coordinates",transMeshCode(i))
node0[[1]][[2]][[2]][[1]][[1]][[1]]<-node4
node0[[1]][[2]][[1]][[1]]<-node3
node2<-node0[[1]][[2]]
nodex<-addChildren(nodex,node2)
}
node0[[1]]<-nodex
}
kmlnode[[1]]<-node0
return(kmlnode)
}
#Sample
library(XML)
doc <- xmlTreeParse("c:/DefaultMesh1.kml")
kmlnode <- xmlRoot(doc)
kmlnode2<-makeMesh1(kmlnode,3622,3643)#2つの区画を指定することで、2区画間の区画をブロック生成する
saveXML(kmlnode2,"c:/meshdata.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
mesh2.JPG

1次メッシュ生成

#Sample2:少々時間がかかります
library(XML)
doc <- xmlTreeParse("c:/DefaultMesh1.kml")
kmlnode <- xmlRoot(doc)
kmlnode2<-makeMesh1(kmlnode,3622,6853)
saveXML(kmlnode2,"c:/meshtest.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
mesh4.JPG

応用編

データを読み込んでkmlを再構成する(1)

csvデータを読み込んで、kmlを再構成しています。JCODEでマッチングさせていますので、事前にデータをノード順と一致させる必要はありません。

library(XML)
source("c:/kmlFunctions.r")
#kmlファイルの読み込み
doc <- xmlTreeParse("c:/gunma.kml")
kmlnode<-xmlRoot(doc)
#csvデータの取得
dat<-read.csv("c:/agedprop.csv")
data<-dat$AP2006
Altdata<-data*1000#高度
name<-"JCODE"#マッチング用指定
#色階層の割り当て
cutnum<-16
classes <- cut(data,seq(min(data),max(data),length=cutnum+1),include.lowest=T)
cols <- colors.kml(rev(heat.colors(cutnum)))
cols2<-cols[classes]
#kmlへの書き込み
len<-kmlNodeLength(doc)-1
for(i in 1:len){
value<-dat$JCODE[[i]]
print(value)
nodeno<-getCDATANo(kmlnode,name,value)#JCODEからノードを取得
print(nodeno)
fillcolor<-cols2[[i]]
linecolor<-NA
kmlnode<-setAltitude(kmlnode,nodeno,Altdata[[i]],fillcolor,linecolor)
}
#kml保存
saveXML(kmlnode,"c:/gunma2.kml",prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')
color3.JPG

データを読み込んでkmlを再構成する(2) 那覇市国勢調査編

まず、那覇市国勢調査のデータ(平成17年)を地図で見る統計(統計GIS)のデータダウンロードから落とします。
http://www.e-stat.go.jp/SG2/toukeichiri/TopFrame.do?fromPage=init&toPage=download
(1)ダウンロードページ

naha1.JPG

(2)世界測地系緯度経度Shape形式を落とします。

naha2.JPG

(3)前準備
落としたShapeFileから、writeOGRでKMLデータ(h17ka47201.kml)を作成します。

library(maptools)
library(rgdal)
x<-readShapePoly("c:/nahashi/h17ka47201.shp")
x@data$KEN_NAME<-iconv(x@data$KEN_NAME,"","utf-8")
x@data$GST_NAME<-iconv(x@data$GST_NAME,"","utf-8")
x@data$MOJI<-iconv(x@data$MOJI,"","utf-8")
tf <- "c:/nahashi/"
writeOGR(x, paste(tf, "h17ka47201.kml", sep=""), "h17ka47201", "KML")

(4)3Dグラフの作成
作成したkmlとShapeFileに付属しているdbfを使って3Dグラフを作成します。
国勢調査データから作成したkmlはMultiGeometry?タグが無いので、setAltitude関数をMG=Fで取り込みます。
ここで、表示させるのは人口のデータです。

library(XML)
library(foreign)
#path<-"/"#Mac
path<-"c:/"
source(paste(path,"kmlFunctions.r",sep=""))
doc <- xmlTreeParse(paste(path,"nahashi/h17ka47201.kml",sep=""))
kmlnode<-xmlRoot(doc)
#データの取得
dat<-read.dbf(paste(path,"nahashi/h17ka47201.dbf",sep=""))
data<-dat$JINKO#人口データ指定(ここを変えることで他の項目も作れます)
Altdata<-data+10#高度:ポリゴンホールが出来ないように全体を10m上げる
name<-"SEQ_NO2"#マッチング用指定
#色階層の割り当て
cutnum<-16
classes <- cut(data,seq(min(data),max(data),length=cutnum+1),include.lowest=T)
cols <- colors.kml(rev(heat.colors(cutnum)))
cols2<-cols[classes]
#kmlへの書き込み
len<-kmlNodeLength(doc)-1
for(i in 1:len){
value<-dat$SEQ_NO2[[i]]#マッチング用指定
print(value)
nodeno<-getCDATANo(kmlnode,name,value)#SEQ_NO2からノードを取得
#nodeno<-i+1 #こちらを使うとかなり速いです(ノードとDBFのデータの並びが同じ場合のみ。SEQ_NO2でマッチングさせてないので速いです。)
print(nodeno)
fillcolor<-cols2[[i]]
linecolor<-NA
#Macの場合以下のコードで実行可能
#label<-getCDATA(kmlnode,"MOJI",nodeno,cdatano=2)#PlaceMark名にMOJIを指定。NMOJIも検索してくるので、相対位置パラメータcdatano=2を指定。
#Windowsの場合、いまのところlabelに日本語指定不可
label<-getCDATA(kmlnode,"KEYCODE1",nodeno)
kmlnode<-setAltitude(kmlnode,nodeno,Altdata[[i]],fillcolor,linecolor,MG=F,label=label)
}
#kml保存
saveXML(kmlnode,paste(path,"nahashi/nahashi3.kml",sep=""),prefix = '<?xml version="1.0" encoding="UTF-8"?>\n')

作成には、かなりの時間がかかりますので、気長にお待ちください。(私のPCで10分)
ただし、nodeno<-i+1の方を使うとすぐに作成できます。(1分くらい)ただし、ノードとDBFのデータの並びが同じ場合のみです。
データは、nahashi3.kmlで落ちてきます。
赤くて高い場所ほど人口が多い地域です。那覇港は手作業で着色しました。

naha5.JPG

那覇市の周辺市町村(浦添市、豊見城市、南風原町)を同時に表示してみた。
明らかに周辺市町村に人口が流れているのがわかる。

city2.JPG

なお、国勢調査データには利用の制限がありますので、ご注意ください。

利用の制限

1. 利用者は、本システムでダウンロードしたデータ及び画像データをそのまま複製(ファイル形式を変換しての複製を含む。)して第三者に譲渡することを禁じます。 なお、商用目的で複製する場合は、e-Statの「お問い合わせ」よりご連絡ください。
2.法律、政令、規則、省令その他すべての法令及び条例等の法規に違反する目的、手段又は方法で、本システムでダウンロードしたデータ及び画像データを利用することを一切禁じます。また、他人の権利を侵害する目的、手段又は方法での利用、公序良俗に反する利用についても一切禁じます。

MacOSXでrgdalを使えるようにする

rgdalパッケージはMacOSXではバイナリで提供されていませんが、下記のサイトを参考にすると、rgdalをMacOSXに導入することができます。(#GDALのインストールと#rgdalのインストールだけやってください)
http://blog.0093.tv/2008/02/mac-os-x-leopardgrassgdalrgdalspgrass6_4571.html
writeOGRもきちんと動きました。(^_^)/

MacOS X でrgdalを使えるようにするもうひとつの方法

  1. http://www.kyngchaos.com/software:frameworksより、GDAL completeをインストール
  2. /Library/Frameworks/GDAL.framework/Programs/gdal-config を /usr/local/bin などにシンボリックリンク
  3. rgdalのソースパッケージをCRANから取得
  4. $ R32 CMD INSTALL rgdal_0.6-20.tar.gz --configure-args='--with-proj-include=/Library/Frameworks/PROJ.framework/Headers --with-proj-lib=/Library/Frameworks/PROJ.Framework/unix/lib'

アクセス数:

3435 人


添付ファイル: filegunma2.kml 1748件 [詳細] fileagedprop.csv 1602件 [詳細] filemesh3.JPG 699件 [詳細] filemesh4.JPG 1441件 [詳細] filemesh1.JPG 687件 [詳細] filecity2.JPG 1462件 [詳細] filekmlFunctions.r 1661件 [詳細] filemesh2.JPG 1361件 [詳細] fileDefaultMesh1.kml 1551件 [詳細] filenaha5.JPG 1476件 [詳細] fileMesh1.kml 816件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1727d)