投稿

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

Rで複数のシェープファイルを結合する

イメージ
地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ 統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位(東区、博多区、中央区、南区、西区、城南区、早良区)の7つのシェープに分かれています。 こういうのをまとめて1つにしたい時ってありますよね。 QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。 ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。 でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この2つを結合させようとする際にエラーになってしまいます。 なので、東区のポリゴンはのIDは、"1.1", "1.2", "1.3", ...、博多区のポリゴンのIDは"2.1", "2.2", "2.3", ... みたいにリネームしたあとに、結合させましょうというのが、↓このスクリプトのポイントでございます。 library ( maptools )   # シェープのファイル名を取得 shape_files <- list.files ( path= "shape_folder" , pattern= ".* \\ .shp" , full.names=T )   # 緯度経度、世界測地系 pj <- CRS ( "+proj=longlat +datum=WGS84" )   # ポリゴンのIDが重複しないように前処理したあと、 # リストに格納しておく shape_li

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  (略) 

Rのデータフレームを降順でソートする(decreasing=Tを使わずに)

イメージ
ほんと、ちょっとした小ネタです。 ↓こちらの本を読んでいて知った(気付いた)方法。 RとRubyによるデータ解析入門 まず、サンプル用のデータフレームを作ります。 name <- c ( "Ana" , "Bill" , "Chris" , "Dan" , "Emily" ) age <- c ( 18 , 59 , 57 , 6 , 31 ) tbl <- data.frame ( name , age ) tbl   name age 1 Ana 18 2 Bill 59 3 Chris 57 4 Dan 6 5 Emily 31 このデータフレームを年齢順にソートしたい場合は、ソート結果のインデックスを返してくれるorder関数を使って、↓こんな風にしますよね。 tbl [ order ( tbl$age ) , ] name age 4 Dan 6 1 Ana 18 5 Emily 31 3 Chris 57 2 Bill 59 order関数のデフォルトは昇順なので、降順にしたい場合は、↓こうですよね。 tbl [ order ( tbl$age , decreasing=T ) , ] name age 2 Bill 59 3 Chris 57 5 Emily 31 1 Ana 18 4 Dan 6 でも、"decreasing=T"って打ち込むのって、けっこう面倒じゃないですか? あと、"decreasing"の綴りが思い付かなかったり(英語力・・・)。 【クイズ】マッチ棒を1本足すだけで、逆にしてください tbl [ order ( - tbl$age ) , ] name age 2 Bill 59 3 Chris 57 5 Emily 31 1 Ana 18 4 D

Rのhistとggplot(geom_histogram)を比較する

イメージ
ggplot2パッケージって使ったことありませんでした。記法がとっつきにくいというのもあったのですが、あの見た目を素直に「きれい」とは認めたくなかったというか。「ggplotの方がかっこいい」という人への反発みたいな部分も多分にあるのですが。 「白地に線画」みたいなビジュアルが好きなんです。下手に色をつけたらかえって趣味が悪くなったりするじゃないですか。自分自身センスがないことを自覚しているので、色やら塗りつぶしやら飾りっぽいものやらは極力使わないようにしていたんですよね。 でも今回、ggplot2パッケージを使ってみようと思ったのは、ベースとなるレイヤーを用意して、次に重ねたいレイヤーを(場合によっては複数)用意して、最後にバンっとplotする、っていうコンセプトがなんだか使いやすいのではないかという気がしたから。なので、徐々に試していこうかなと思っています。 前置きが長くなりました。 今回はヒストグラムですが、多くの人に馴染みのあるデフォルトのhistでの描画と、ggplot2パッケージを使った場合とで、どんな風に違うのか試してみました。 正規分布に従う乱数を発生させて、普通にヒストグラムを描いてみます。 data1 <- rnorm ( 1000 )   hist ( data1 ) 通常のhist関数で描画 ggplot2パッケージを使うと↓こんな描き方になります。 install.packages ( "ggplot2" ) library ( ggplot2 )   df1 <- data.frame ( data1 ) # データフレームにしておく必要がある   g <- ggplot ( df1 , aes ( x=data1 ) ) # 対象がdata1列であることを指定 g <- g + geom_histogram ( ) # ヒストグラムを描画することを指定 plot ( g ) # 描画 ggplot2パッケージを使ったもの 先ほどの、デフォルトのhistと階級幅を合わせる(=0.5)には、binwidthを指定します↓ g <- g