投稿

11月, 2015の投稿を表示しています

Rでベクトル場を図示する

イメージ
ベクトル場がxy座標の式で与えられているときに、それがどんな場を表しているかは、具体的にいくつかのベクトルを図示してみると、なんとなく見えてきたりします。   f(x, y) = (x, y) みたいな単純なベクトル場でさえ、慣れてないとイメージわかなかったりしますよね。 Rで図示してみましょう。(唐突) # # ベクトル場を表す関数 # f <- function ( x , y ) { X <- x Y <- y return ( c ( X , Y ) ) }   # # 矢印を描画する関数 # write_arrow <- function ( x , y ) { a <- 0.1 # 矢の長さ(比率) b <- 0.05 # 矢尻の長さ arrow_head <- f ( x , y ) * a # aの値で長さを調整   arrows ( x , # 矢のx座標(from) y , # 矢のy座標(from) x + arrow_head [ 1 ] , # 矢のx座標(to) y + arrow_head [ 2 ] , # 矢のy座標(to) length = b ) # 矢尻のサイズ }   # # ここからメイン # n <- 10 # 格子の1辺の数   # 描画する領域を準備 plot ( 0 , 0 , xlim= c ( -n , n ) , ylim= c ( -n , n ) , type= "n" , xlab= "x" , ylab= "y" )   # 各格子点について繰り返す for ( x in -n:n ) { for ( y in -n:n ) { write_arrow ( x , y ) } } f(x, y) = (x, y)のベクトル場 「ぶ

RでFukuoka City Wi-Fiの利用状況を可視化する

イメージ
オープンデータとして、Fukuoka City Wi-Fiの利用状況のCSVファイルが公開されています↓ Fukuoka City Wi-Fi 利用状況 - データセット - 福岡市のCKAN が、例によってなかなか扱いにくそうなフォーマットになっております↓ 福岡市のオープンデータサイトで公開されていたCSVの形式 1拠点の1日が1行になっており、その行の中に0時台、1時台、2時台、~、23時台のように横に伸びてデータが格納されています。別の日は、別の行に格納されています。 そのままじゃ扱いにくそうなので、日付時刻型で表現して、形もシンプルにしましょう。 d <- read.csv ( "201410access.csv" )   d <- d [ d$X != "合計" , ] # 合計が入っているレコードは削除する   df <- data.frame ( ) # 空のデータフレームを用意しておく   for ( i in 1 : nrow ( d ) ) { tmp <- paste ( d [ i , ] $X , sprintf ( "%02d" , 0 : 23 ) ) YmdH <- strptime ( tmp , "%Y%m%d %H" )   value <- as.vector ( t ( d [ i , 4 : 27 ] ) )   df <- rbind ( df , data.frame ( 年月日時 = YmdH , 拠点名 = d [ i , ] $拠点名 , 認証回数 = value ) ) } そうすると、冗長だけどシンプルな形式になりました↓ CSVの形式を変えてみたもの ひと月の合計値で棒グラフを書いてみましょう。認証回数の少ない順にソートしています。なんとなくggplotを使ってみました。(いまいち、まだ慣れていない感が・・・) ggplotの

Rで、福岡市オープンデータサイトのCSVファイルを扱いやすいように加工する

イメージ
福岡市のオープンデータのサイト に「福岡市の大気環境測定結果(直近48時間)」というCSVデータが公開されています。 Rでこれを取り込みたいときは、 d <- read.csv("http://ckan.open-governmentdata.org/dataset/ef0a812c-538b-4e54-ae3e-b22022a7c83a/resource/c60ccc29-811d-44f8-9aba-dcb660387927/download/kankyodata48.csv") とするだけでOK。(URLが長くてごちゃごちゃしてますが、read.csvにダブルクォートで指定しているだけです) でも、中身を見てみると、いまいち使いにくそうなフォーマットなんですよね。 head(d)とやれば分かりますが、ちょっと列が多いので、列名だけ表示させてみます。 > colnames(d) [1] "年月日" "測定局名称" "測定項目名称" "緯度" [5] "経度" "単位" "測定値.1時." "測定値.2時." [9] "測定値.3時." "測定値.4時." "測定値.5時." "測定値.6時." [13] "測定値.7時." "測定値.8時." "測定値.9時." "測定値.10時." [17] "測定値.11時." "測定値.12時." "測定値.13時." "測定値.14時." [21] "測定値.15時." "測定値.16時." "測定値.17時." "測定値.18時." [25] "測定値.19時." "測定値.20

Rで、変数の値を変数名や列名にする

タイトルの意味分かりました? 例えば、x00, x01, x02, ..., x99 という名前の100個の変数を作りたいとか、それらに値を設定したいというときに使える方法です。(そんなシチュエーションがあるかどうかは別として) 上記の変数に、0, 1, 2, ..., 99 の値を設定したいとき、 x00 <- 0 x01 <- 1 x02 <- 2 ... x99 <- 99 を全部書けばいいわけですが、それは大変なので、下記のようにできます。 for(i in 0:99){   # 命令の文字列を作る   s <- paste("x", sprintf("%02d",i), " <- ", i, sep="")   # 文字列を命令として実行する   eval(parse(text=s)) } ls() # x00~x99の100個の変数ができていることを確認 このevalとparseのコンビを、列名に使いたい場合は、↓こんな感じでどうでしょうか。 ちなみにRには"letters"なる定数が用意されていて、その中身にはaからzまでの文字が格納されています。 > letters [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" [11] "k" "l" "m" "n" "o" "p" "q" "r" "s" "t" [21] "u" "v" "w" "x" "y" "z" これらを列名にもつ、データフレームを作ってみましょう。 &g

Rで10進数→16進数、16進数→10進数の変換

イメージ
10進数から16進数に変換する方法の1つとして、sprintfを使う方法があります。 > sprintf(" %X ", 0:31)    [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  [9] "8"  "9"  "A"  "B"  "C"  "D"  "E"  "F" [17] "10" "11" "12" "13" "14" "15" "16" "17" [25] "18" "19" "1A" "1B" "1C" "1D" "1E" "1F" 16進といえば、色のRGB値を指定するときに使ったりしますが、その場合は "7" とか "F" ではなく "07" とか "0F" とかの形で欲しいですよね。 そういうときは、 > sprintf(" %02X ", 0:31)  [1] "00" "01" "02" "03" "04" "05" "06" "07"  [9] "08" "09" "0A" "0B" "0C" "0D" "0E" "0

Rを使って郵便番号界の地図を作る(統計局のShape+日本郵便のCSV)

イメージ
郵便番号がどんな感じで割り当てられているか、地図上で確認したかったんです。 GPSで取った地点情報がプライバシーに配慮云々で郵便番号になっていた、なんてことがありまして、その粒度がどんなものか確認したかったってのが発端です。 検索してみても無償のものは意外に見つからないもので、国際航業さんが売っていたりします↓ PAREA-Zip郵便番号界 | 国際航業株式会社 ノースリーブな(=無い袖は振れない)わたくしとしては、自分で作るしかないわけで。 郵便番号のリストは↓ここにありました。 郵便番号データダウンロード - 日本郵便 住所文字列と郵便番号が入っていればどれでもいいんですが、とりあえず、「住所の郵便番号(CSV形式)」の「読み仮名データの促音・拗音を小書きで表記するもの」というのをダウンロードしました。 次は地図です。こういうときに使う境界データとしては、統計局のe-StatでShapeファイルが公開されているので、これが便利。入手の方法は下記を参照。 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ 今回は、我らがサラリーマンの聖地 新橋を有する港区を対象にしてみました。 とりあえず地図をプロットするところまでやってみましょう。 library(maptools)     # GIS系のパッケージ # 座標系の情報を持つオブジェクトを作る # longlatは緯度経度指定、WGS84は世界測地系であることを表す pj <- CRS("+proj=longlat +datum=WGS84") # シェープファイルの読み込み map <- readShapePoly("h22ka13103.shp", proj4string=pj) # とりあえずプロット plot(map) shapeをプロットしただけのもの この地図上に追加で郵便番号をプロットすればよさそうですね。 読み込んだシェープファイルにはポリゴンの境界情報だけでなく、それぞれのポリゴンの属性の情報も含まれています。 こうやれば確認できます head(map@data)        AREA  (略)