ラベル R の投稿を表示しています。 すべての投稿を表示
ラベル R の投稿を表示しています。 すべての投稿を表示

2014年11月22日土曜日

plotKMLのメモ

plotKML(http://plotkml.r-forge.r-project.org/)が面白い。

RでKMLを処理するためのパッケージ、簡単にKMLを生成できてしまう。

たとえば、こんな感じの可視化が簡単にできる。



Rスクリプトはこんな感じ。


library(maptools)
library(plotKML)
shp<-readShapePoly("h22ka18207.shp")
proj4string(shp) <- CRS("+proj=longlat +datum=WGS84")
z<-shp@data$JINKOz[z==0<-1
kml(shp,labels=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),altitude=z,colour="#ff00ff",alpha=0.75,plot.lab=TRUE)




Shapefileはe-Statからダウンロードしてきたもの使った。

バブルチャートっぽく表示する場合はこんな感じ。

data <- data.frame(x=shp@data$X_CODE,y=shp@data$Y_CODE,id=shp@data$KEY_CODE,label=iconv(shp@data$MOJI,to="UTF8",from="SJIS"),jinko=shp@data$JINKO)
coordinates(data) <- ~x+y
proj4string(data)<-CRS("+proj=longlat +datum=WGS84")
kml(data,shape="http://maps.google.com/mapfiles/kml/pal2/icon18.png",color="#ff0000",size=jinko,labels=label)




公式ページをみるとアニメーションにも対応してる模様。


KMLを書き起こすのは面倒だけど、これなら簡単。

何かの時に役立つかもしれない。










2013年11月10日日曜日

RでSPARQL、結果はKML

RでSPARQLが使えないかと検索してみたら、あった。

パッケージ名はそのままのSPARQL。

さっそく、dbpedia.orgのEndpointを使ってテストしてみた。

クエリは以前に使った日本の都道府県の人口を検索するもの。

せっかくなので、検索結果はKMLで出力する。

Google Earthで表示してみると、こんな感じ。


ちょっとデータのウェブを使っている感じがある。

昨今のオープンデータをSPARQLで探してこれるようにすると、面白いかもしれない。(read.csv等で読みこめばいいのだけど、ファイル自体を探してこないといけない)

むかしExcelをRDFに変換して、SPARQLで検索、グラフで表示ということをやっていたけど、あの苦労が嘘のよう。

以下はRのソース。

library(SPARQL)
library("sp")
library("rgdal")
url<-"http://dbpedia.org/sparql"
query<-"PREFIX rdf: <http://www.w3.org/1999/02/22-rdf-syntax-ns#>
PREFIX rdfs: <http://www.w3.org/2000/01/rdf-schema#>
PREFIX wgs84_pos: <http://www.w3.org/2003/01/geo/wgs84_pos#>
PREFIX dbpedia-owl: <http://dbpedia.org/ontology/>
PREFIX yago: <http://dbpedia.org/class/yago/>
SELECT DISTINCT ?title ?lat ?lon ?populationTotal where{
?s rdf:type yago:PrefecturesOfJapan .
?s rdfs:label ?title .
? S dbpedia-owl: populationTotal ?populationTotal .
?s wgs84_pos:lat ?lat .
?s wgs84_pos:long ?lon .
filter langMatches(lang(?title),\"en\") .
}"
res<-SPARQL(url,query)
d<-res$result
coordinates(d)<-c('lon','lat')
proj4string(d)<-CRS("+init=epsg:2451")
writeOGR(d, dsn="hoge.kml", layer="sample", driver="KML")

2013年6月22日土曜日

Rでアニメーション

Rのanimationライブラリを使ってみた.

plotで出力した画像をくつけてアニメーションにする.

こんな感じ,

library(animation)
f<-4
wave<-function(n){w<-0;t<-seq(0,1,0.01);for(k in 1:n){w<-w+sin((2*k-1)*2*pi*f*t)/(2*k-1)};w<-4/pi*w}
wave2<-function(x){for(n in 1:x){plot(x=seq(0,1,0.01),y=wave(n),type="l",main=paste("Square wave(n=",n,")"),xlab="Time(sec)",ylab="Amplitude")}}
saveMovie(wave2(30),interval=0.5,moviename="wave",movietype="gif",outdir=getwd(),width=640,height=480)

で,できたアニメーションはこちら.

グラフの変化をみせるにはいいツールだね.

2013年5月31日金曜日

RでKMLファイルを出力する

RとGoogle Earthを連携させることができるらしい.

検索してみると,縫村氏の「RとGoogle Earth を使った空間情報の可視化」と題したスライド資料を発見.

スライド資料ではRを使ってアメリカ地質調査所(USGS)の地震動データ(CSV)からKMLファイルを作成するやり方が紹介されていた.

こりゃ面白そうということで早速トライ.

DataFrameをSpatialPointsDataFrameに変換するところでハマったけど,proj4string(qk)<-CRS("+init=epsg:4326")とすることで解決.

無事,作成したKMLファイルをGoogle Earthに読み込ませることができました.

ちょっと見にくいけど,地震の発生箇所に目印が表示されています.


うーん.こんなことが簡単にできてしまうなんてよい世の中になったものです.

Google Earthの使い方がわかってきました.こりゃ迫力あるわー.

Rのコンソール出力は以下のとおり.

> library(sp)
> library(rgdal)
> qk<-read.csv("http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.csv")
> head(qk)
                      time latitude longitude  depth mag magType nst  gap
1 2013-05-31T13:40:41.000Z  60.6428 -147.5946   3.60 3.4      Ml  NA   NA
2 2013-05-31T13:40:28.000Z  61.5455 -148.6714  28.50 2.9      Ml  NA   NA
3 2013-05-31T13:25:55.190Z -28.2755 -178.6027 262.93 5.2      mb  46 66.0
4 2013-05-31T13:24:51.800Z  33.6878 -119.1308   0.40 3.1      Ml  24 75.6
5 2013-05-31T13:06:50.150Z -20.3071  169.0243  34.03 5.1      mb  57 53.0
6 2013-05-31T11:32:41.100Z  40.1692 -121.1097   0.00 3.1      Md  NA 46.8
       dmin  rms net         id                  updated
1        NA 1.02  ak ak10727802 2013-05-31T13:49:58.971Z
2        NA 1.12  ak ak10727799 2013-05-31T13:48:21.339Z
3 1.1300000 0.94  us usb000hacp 2013-05-31T13:50:51.623Z
4 0.2245788 0.23  ci ci15352009 2013-05-31T13:32:07.619Z
5 2.7200000 1.25  us usb000hacc 2013-05-31T13:42:49.516Z
6 0.2155957 0.38  nc nc72003416 2013-05-31T13:15:16.232Z
> coordinates(qk)<-c("latitude","longitude")
> class(qk)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
> summary(qk)
Object of class SpatialPointsDataFrame
Coordinates:
                min      max
longitude -178.6027 169.0243
latitude   -28.2755  66.2616
Is projected: NA 
proj4string : [NA]
Number of points: 33
Data attributes:
                       time        depth             mag       
 2013-05-30T14:09:44.000Z: 1   Min.   :  0.00   Min.   :2.500  
 2013-05-30T14:38:25.650Z: 1   1st Qu.: 10.00   1st Qu.:3.000  
 2013-05-30T14:54:01.000Z: 1   Median : 28.50   Median :3.300  
 2013-05-30T17:46:40.000Z: 1   Mean   : 46.29   Mean   :3.739  
 2013-05-30T18:37:00.620Z: 1   3rd Qu.: 67.70   3rd Qu.:4.600  
 2013-05-30T20:41:29.000Z: 1   Max.   :262.93   Max.   :5.300  
 (Other)                 :27                                   
  magType        nst              gap             dmin         
 mb   :13   Min.   :  5.00   Min.   : 33.0   Min.   : 0.01797  
 mb_Lg: 1   1st Qu.: 16.00   1st Qu.: 75.6   1st Qu.: 0.73033  
 Md   : 6   Median : 28.00   Median :109.0   Median : 1.21000  
 Ml   :12   Mean   : 36.15   Mean   :136.1   Mean   : 3.43408  
 Mw   : 1   3rd Qu.: 46.50   3rd Qu.:174.1   3rd Qu.: 2.47250  
            Max.   :109.00   Max.   :324.0   Max.   :32.25000  
            NA's   :13       NA's   :9       NA's   :9         
      rms         net              id                         updated  
 Min.   :0.1200   ak: 9   ak10726951: 1   2013-05-30T14:22:58.908Z: 1  
 1st Qu.:0.3100   ci: 1   ak10727009: 1   2013-05-30T16:41:03.895Z: 1  
 Median :0.8200   hv: 1   ak10727095: 1   2013-05-30T16:56:39.447Z: 1  
 Mean   :0.7391   nc: 3   ak10727182: 1   2013-05-30T18:00:03.809Z: 1  
 3rd Qu.:1.0200   pr: 5   ak10727211: 1   2013-05-30T20:39:31.117Z: 1  
 Max.   :1.2900   us:14   ak10727538: 1   2013-05-30T22:47:39.888Z: 1  
                          (Other)   :27   (Other)                 :27  
> proj4string(qk)<-CRS("+init=epsg:4326")
> summary(qk)
Object of class SpatialPointsDataFrame
Coordinates:
                min      max
longitude -178.6027 169.0243
latitude   -28.2755  66.2616
Is projected: FALSE 
proj4string :
[+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0]
Number of points: 33
Data attributes:
                       time        depth             mag       
 2013-05-30T14:09:44.000Z: 1   Min.   :  0.00   Min.   :2.500  
 2013-05-30T14:38:25.650Z: 1   1st Qu.: 10.00   1st Qu.:3.000  
 2013-05-30T14:54:01.000Z: 1   Median : 28.50   Median :3.300  
 2013-05-30T17:46:40.000Z: 1   Mean   : 46.29   Mean   :3.739  
 2013-05-30T18:37:00.620Z: 1   3rd Qu.: 67.70   3rd Qu.:4.600  
 2013-05-30T20:41:29.000Z: 1   Max.   :262.93   Max.   :5.300  
 (Other)                 :27                                   
  magType        nst              gap             dmin         
 mb   :13   Min.   :  5.00   Min.   : 33.0   Min.   : 0.01797  
 mb_Lg: 1   1st Qu.: 16.00   1st Qu.: 75.6   1st Qu.: 0.73033  
 Md   : 6   Median : 28.00   Median :109.0   Median : 1.21000  
 Ml   :12   Mean   : 36.15   Mean   :136.1   Mean   : 3.43408  
 Mw   : 1   3rd Qu.: 46.50   3rd Qu.:174.1   3rd Qu.: 2.47250  
            Max.   :109.00   Max.   :324.0   Max.   :32.25000  
            NA's   :13       NA's   :9       NA's   :9         
      rms         net              id                         updated  
 Min.   :0.1200   ak: 9   ak10726951: 1   2013-05-30T14:22:58.908Z: 1  
 1st Qu.:0.3100   ci: 1   ak10727009: 1   2013-05-30T16:41:03.895Z: 1  
 Median :0.8200   hv: 1   ak10727095: 1   2013-05-30T16:56:39.447Z: 1  
 Mean   :0.7391   nc: 3   ak10727182: 1   2013-05-30T18:00:03.809Z: 1  
 3rd Qu.:1.0200   pr: 5   ak10727211: 1   2013-05-30T20:39:31.117Z: 1  
 Max.   :1.2900   us:14   ak10727538: 1   2013-05-30T22:47:39.888Z: 1  
                          (Other)   :27   (Other)                 :27  
> writeOGR(qk, dsn="qk.kml", layer="sample", driver="KML")

2012年11月1日木曜日

Rで波形を描画する




Rを使って色々な波形を描画しようと思う.

まずは基本パラメータを設定する.

サンプリング回数は44.1kHz,
描画する波形の周波数は440Hzとする.
作成する波形は1秒間とする.
freq=44100;
nn=1:freq;
tt=nn/freq;
f1=440;
描画用の関数とクリッピング用の関数を定義する.

draw<-function(s,sec){
freq=44100;
nn=1:(freq*sec);
tt=nn/freq;    
plot(tt,s[1:(freq*sec)],type="l",col="blue");
}
cliping<-function(s){
for(i in 1:length(s)){
           if(s[i] > 1.0 ) s[i]=1.0;
           if(s[i] < -1.0 ) s[i]=-1.0;
}
return(s);
}

これで準備は整った.

まずは正弦波を描画する.
描画する波形は0.01秒分としている.
http://en.wikipedia.org/wiki/Sine_wave
s=sin(2*pi*tt*f1);
s=cliping(s);
draw(s,0.01);
つぎは三角波.
http://en.wikipedia.org/wiki/Triangle_wave
s=tt*0;
for(k in 0:10){
s=(-1)^k*(sin((2*k+1)*2*pi*f1*tt)/(2*k+1)^2)+s;
}
s=8/pi^2*s;
s=cliping(s);
draw(s,0.01);
つぎは矩形波.
http://en.wikipedia.org/wiki/Square_wave
s=tt*0;
for(k in 1:1000){
      s=s+sin((2*k-1)*pi*f1*tt)/(2*k-1);
}
s=4/pi*s;
s=cliping(s);
draw(s,0.01);
最後はノコギリ波.

http://en.wikipedia.org/wiki/Sawtooth_wave

s=tt*0;
for(k in 1:1000){
s=s+sin(2*pi*k*f1*tt)/k;
}
s=2/pi*s;
s=cliping(s);
draw(s,0.01);

こんなことが簡単にできてしまうのだからRは面白い.





波形データは以下で書き出せる.
write(s,"wav.txt",ncolumns=1);
書き出したファイルはtxt2wavでwavファイルに変換できる.

txt2wav.cのコンパイルは
$ gcc -m32 -o txt2wav txt2wav.c 
な感じで(64bit環境の場合ね).

2011年12月4日日曜日

ダイクストラ法

ダイクストラ法をRとigraphを使って表現する.

igraphのshortest.paths関数で一発でできるのだけど,お勉強のためということで.

まずは最短経路を求めるグラフを用意する.
library(igraph)
g<-graph.empty(5)
V(g)$name<-c("A","B","C","D","E")
g<-add.edges(g,c(0,1,0,2,0,3,1,2,1,3,1,4,2,3,3,4))
g<-as.undirected(g)
E(g)$weight<-c(2,1,7,5,6,1,1,8)
lay<-rbind(c(1,3),c(2,3),c(1,1),c(2,1),c(3,2))
plot(g, layout=lay,vertex.label=V(g)$name,vertex.size=12,edge.label=E(g)$weight)


こんな感じ.今回はノード「E」をゴールとする.

つぎに各ノードのコストを初期化する.
V(g)$cost<-c(Inf,Inf,Inf,Inf,0)
V(g)$fixed=FALSE
V(g)$color<-"skyblue"
plot(g, layout=lay,vertex.label=V(g)$cost,vertex.size=12,edge.label=E(g)$weight,vertex.color=V(g)$color)
つぎはコスト計算.
V(g)[nei(V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]) & fixed==FALSE][V(g)[nei(V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]) & fixed==FALSE]$cost>(E(g)[V(g)[fixed==FALSE] %->% V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]]$weight)+(V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]$cost)]$cost<-(E(g)[V(g)[fixed==FALSE] %->% V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]]$weight)+(V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]$cost)
ながい...もっと短くできるかもしれない.
V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]$color="yellow"
V(g)[fixed==FALSE & cost==min(V(g)[fixed==FALSE]$cost)]$fixed=TRUE
計算が終わったノードは黄色に塗りつぶす.
plot(g, layout=lay,vertex.label=V(g)$cost,vertex.size=12,edge.label=E(g)$weight,vertex.color=V(g)$color)

これを計算するノードが無くなるまで繰り返す.
もう一回.

もう一回.
もう一回.

もう一回.


これで終わり.

実際にグラフを描きながらだと楽しいね.








2011年12月1日木曜日

R + igraph + Cytoscape


まずはR+igraphでネットワークグラフを作成する.

> library(igraph)
> d<-matrix(c(0,0,1,1,1 ,0,1,1,0,2,0,3,1,0,1,6,0,0,0,2,1,8,5,2,0),ncol=5,nrow=5)
> ls()
[1] "d"
> d
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    6    1
[2,]    0    1    3    0    8
[3,]    1    1    1    0    5
[4,]    1    0    0    0    2
[5,]    1    2    1    2    0
> g<-graph.adjacency(d,weighted=TRUE)
> summary(g)
Vertices: 5
Edges: 15
Directed: TRUE
No graph attributes.
No vertex attributes.
Edge attributes: weight.
> degree(g)
[1] 5 6 7 4 8
> png("c:/tmp/myplot.png",width=600, height=600,pointsize=14)
> plot(g,layout=layout.fruchterman.reingold,vertex.label=V(g)$smr,vertex.size=degree(g)*4)
> dev.off()
null device
          1
なかなか,いい感じ.

つぎはCytoscape(http://www.cytoscape.org/)を使って描画してみる.

まずは,Rからグラフをエキスポートする.

このときファイルタイプとして”ncol”を指定するのを忘れずに!

> write.graph(g,"c:/tmp/mygraph.ncol","ncol")



Cytoscapeを起動して[File]>[Import]>[Network from Table(Text/MS Excel...)]を選択する.

先ほどエキスポートした「c:/tmp/mygraph.ncol」を読み込む.

[Source Interaction]に「Column 1」,[Interaction Type]に「Column 3」,[Target Interaction]に「Column 2」を設定する.


[Import]ボタンをクリックする.

これでグラフをインポートできた.

あとはCytoscpae上でグラフを色々と見やすくしてやればいい.





2011年11月27日日曜日

R+igraphのお勉強 その3

ノードの次数を確認する.
> g
Vertices: 10
Edges: 9
Directed: TRUE
Edges:
         
[0] 0 -> 1
[1] 0 -> 2
[2] 1 -> 3
[3] 1 -> 4
[4] 2 -> 5
[5] 2 -> 6
[6] 3 -> 7
[7] 3 -> 8
[8] 4 -> 9
次数を確認する.

> degree(g)
 [1] 2 3 3 3 2 1 1 1 1 1

入次数を確認する.
> degree(g,mode="in")
 [1] 0 1 1 1 1 1 1 1 1 1
出次数を確認する.
> degree(g,mode="out")
 [1] 2 2 2 2 1 0 0 0 0 0
 当然,無向グラフに入出の区別はない.
> g2
Vertices: 10
Edges: 9
Directed: FALSE
Edges:
         
[0] 0 -- 1
[1] 0 -- 2
[2] 1 -- 3
[3] 1 -- 4
[4] 2 -- 5
[5] 2 -- 6
[6] 3 -- 7
[7] 3 -- 8
[8] 4 -- 9
次数を確認する.
> degree(g2)
 [1] 2 3 3 3 2 1 1 1 1 1
> degree(g2)
 [1] 2 3 3 3 2 1 1 1 1 1
> degree(g2,mode="in")
 [1] 2 3 3 3 2 1 1 1 1 1
> degree(g2,mode="out")
 [1] 2 3 3 3 2 1 1 1 1 1
ノードの次数をノードのサイズに設定してグラフを描画する.
> png("plot3.png",width=400, height=400,pointsize=12)
次数の10倍をノードサイズとする.
> plot(g,layout=lay,vertex.size=degree(g)*10)
> dev.off()
null device
          1


R+igraphのお勉強 その2

前回の続き.

有向グラフを無向グラフに変換してみる.

有向グラフの確認.
> is.directed(g)
[1] TRUE
無向グラフに変換する.
> g2<-as.undirected(g)
変換結果を確認する.
> g2
Vertices: 10
Edges: 9
Directed: FALSE
Edges:
         
[0] 0 -- 1
[1] 0 -- 2
[2] 1 -- 3
[3] 1 -- 4
[4] 2 -- 5
[5] 2 -- 6
[6] 3 -- 7
[7] 3 -- 8
[8] 4 -- 9
> is.directed(g2)
[1] FALSE
グラフを描画する.
> png("plot2.png",width=400, height=400,pointsize=12) 
> plot(g2,layout=lay)
> dev.off()
null device
          1
作業経過をファイル保存する.
 > save(g,g2,lay,file="test.dat")

R+igraphのお勉強 その1


igraphライブラリを読み込む.
> library(igraph)
10ノードのツリーグラフを作成する.
> g<-graph.tree(10)
> summary(g)
Vertices: 10
Edges: 9
Directed: TRUE
No graph attributes.
No vertex attributes.
No edge attributes.
> g
Vertices: 10
Edges: 9
Directed: TRUE
Edges:
         
[0] 0 -> 1
[1] 0 -> 2
[2] 1 -> 3
[3] 1 -> 4
[4] 2 -> 5
[5] 2 -> 6
[6] 3 -> 7
[7] 3 -> 8
[8] 4 -> 9
グラフをレイアウトする. 
> lay<-layout.fruchterman.reingold(g)
> lay
            [,1]       [,2]
 [1,]  0.5750165 -0.7339218
 [2,]  4.2972271 -0.9157740
 [3,] -2.8644359  0.1233019
 [4,]  5.7313203  2.7760283
 [5,]  5.3732637 -4.4586431
 [6,] -3.4715509  3.4824008
 [7,] -4.1404984 -2.9629157
 [8,]  4.2804038  5.9832150
 [9,]  8.9190178  2.4068283
[10,]  4.1109849 -7.2698722
 ここまでの作業をファイル保存する.
> save(g,lay,file="test.dat")
グラフを描画する.
> plot(g,layout=lay)
> png("plot.png",width=400, height=400,pointsize=12)
> plot(g,layout=lay)
> dev.off()
null device
          1

2011年11月20日日曜日

グラフレイアウト

とりあえずライブラリをロードする.
> library(igraph)
とりあえずツリーグラフを作成する.
> g<-graph.tree(4)
> plot.igraph(g)
グラフをレイアウトする.
> lay<-layout.kamada.kawai(g)
> lay
             [,1]       [,2]
 [1,] -1.95786789 -0.6726298
 [2,] -0.97281249 -1.2272012
 [3,] -2.94060962 -0.1762201
 [4,]  0.09102538 -0.7586657
 [5,] -1.03996471 -2.3633110
 [6,] -3.22145269  0.8896139
 [7,] -3.96005388 -0.5896342
 [8,]  0.59269741  0.2250630
 [9,]  1.10354845 -1.2044752
[10,] -1.05079882 -3.4245322
再描画.
> plot.igraph(g,layout=lay)



ちゃんとレイアウトされた.

2011年11月19日土曜日

Rでネットワーク図を描く3

 今度は行ラベル,列ラベルを付けてみる.
> d<-matrix(c(0,0,1,1,1 ,0,1,1,0,2,0,3,1,0,1,6,0,0,0,2,1,8,5,2,0),ncol=5,nrow=5
> d
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    6    1
[2,]    0    1    3    0    8
[3,]    1    1    1    0    5
[4,]    1    0    0    0    2
[5,]    1    2    1    2    0
 まずは列ラベルを設定する.
> colnames(d)<-c("A","B","C","D","E")
> d
     A B C D E
[1,] 0 0 0 6 1
[2,] 0 1 3 0 8
[3,] 1 1 1 0 5
[4,] 1 0 0 0 2
[5,] 1 2 1 2 0
お次は行ラベルを設定する.
> rownames(d)<-c("A","B","C","D","E")
> d
  A B C D E
A 0 0 0 6 1
B 0 1 3 0 8
C 1 1 1 0 5
D 1 0 0 0 2
E 1 2 1 2 0
グラフに変換する.
> g<-graph.adjacency(d,weighted=TRUE)
画像に書き出す.
> png("c:/tmp/myplot.png",width=400, height=400,pointsize=12)
> plot(g)
> dev.off()
null device
          1 

エッジの重みを表示する.
> E(g)$weight
 [1] 6 1 1 3 8 1 1 1 5 1 2 1 2 1 2
一応,グラフ情報もみておく.
> summary(g)
Vertices: 5
Edges: 15
Directed: TRUE
No graph attributes.
Vertex attributes: name.
Edge attributes: weight.

今度は特定のエッジを削除してみる.

重みが5以下のエッジを削除してみよう.
> g2<-delete.edges(g,E(g)[weight<=5])
グラフ情報をみると,
> summary(g2)
Vertices: 5
Edges: 2
Directed: TRUE
No graph attributes.
Vertex attributes: name.
Edge attributes: weight.
描画してみよう.
>  png("c:/tmp/myplot.png",width=400, height=400,pointsize=12)
> plot(g2,layout=layout.fruchterman.reingold,vertex.color="white",vertex.label=V(g2)$name,edge.label=E(g2)$weight,vertex.size=10)
> dev.off()
null device
          1



次数0のノードができてしまった.

こいつを削除しよう.
> g2<-delete.vertices(g2,which(degree(g2)<1)-1)
> summary(g2)
Vertices: 4
Edges: 2
Directed: TRUE
No graph attributes.
Vertex attributes: name.
Edge attributes: weight.
もう一回,書き出してみる.
>  png("c:/tmp/myplot.png",width=400, height=400,pointsize=12)
> plot(g2,layout=layout.fruchterman.reingold,vertex.color="white",vertex.label=V(g2)$name,edge.label=E(g2)$weight,vertex.size=10)
> dev.off()
null device
          1

 ちゃんと削除できた.

2011年11月18日金曜日

Rでネットワーク図を描く2

こんどはもう少し大きめのグラフを書いてみる.

> library(igraph)
> d<-matrix(c(0,0,1,1,1 ,0,1,1,0,2,0,3,1,0,1,6,0,0,0,2,1,8,5,2,0),ncol=5,nrow=5
> d
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    6    1
[2,]    0    1    3    0    8
[3,]    1    1    1    0    5
[4,]    1    0    0    0    2
[5,]    1    2    1    2    0
> g<-graph.adjacency(d,weighted=TRUE)
> png("c:/tmp/myplot.png",width=400, height=400,pointsize=12)
> plot(g,,layout=layout.fruchterman.reingold,vertex.color="white",vertex.label=V(g)$name,edge.color="black",edge.label=E(g)$weight,vertex.size=10)
> dev.off()
null device
          1
重みをエッジの幅に設定する.
 > E(g)$weight
 [1] 6 1 1 3 8 1 1 1 5 1 2 1 2 1 2
> E(g)$width<-E(g)$weight
美しくないのでやっぱり戻そう.

> E(g)$width<-1

Rでネットワーク図を描く

Rが面白い.簡単に色々とできてしまう.

今回はigraphライブラリでネットワーク図の作成にトライしてみる.

まずライブラリを読み込む.
> library(igraph)
つぎに隣接行列を作る.
> d<-matrix(c(0,0,1,1,0,1,1,0,0),ncol=3,nrow=3)
行列のサイズを確認する.
> dim(d)
[1] 3 3
行列を表示する.
> d
     [,1] [,2] [,3]
[1,]    0    1    1
[2,]    0    0    0
[3,]    1    1    0
 行列から重み付き有向グラフを作成する.
> g<-graph.adjacency(d,mode="directed",weighted=TRUE)
グラフ情報を確認する.
> summary(g)
Vertices: 3
Edges: 4
Directed: TRUE
No graph attributes.
No vertex attributes.
Edge attributes: weight.
最後にグラフを画像ファイルに書き出す.
> png("c:/tmp/myplot.png",width=400, height=400,pointsize=12)
> plot(g,vertex.color="white",vertex.label=V(g)$name,edge.color="black",edge.label=E(g)$weight,vertex.size=10)
> dev.off()
null device
          1