Rプログラミングの小ネタ

福岡県の市単位（政令指定都市は区単位）の人口密度単位のランキングは、↓こちらにあって、 福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】 　　1位 福岡市 中央区（11770人/㎢） 　　2位 福岡市 城南区（ 8031人/㎢） 　　3位 福岡市 南区　（ 7976人/㎢） となっております。トップは予想通りの福岡市中央区。 ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。 ということで、Rを使って町丁目単位（博多駅前１丁目とか、天神２丁目とか）でのランキングを出してみました。 とりあえず結果をどうぞ↓ 中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸２丁目で人口密度を計算すると 49294人/㎢ と、4倍以上の数字を叩き出しました。 また、上位5位の荒戸、平尾、薬院、高砂はいずれも中央区の町名ですが、6位に入っている愛宕浜は西区にある町です。 西区と言えば福岡市の中では辺境ですし（注：著者の個人的なイメージです）、前述の市区単位のランキングでは17位とかなり下位にランキングされてました、 実は、この愛宕浜2丁目は室見川をはさんで早良区のすぐ対岸にあって、比較的都市部に近い場所なんですよね。そして、ほとんどマンションしかありません（他には中学校とマルキョウがある）。大抵のマンションは10階以上あるので、おのずと人口密度が高くなるわけですね。 ↓Google Earthで見てみるとこんな感じ では、最後にRスクリプトを載せておきます。 シェープファイルは↓こちらに書いてある方法で入手して、 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ ↓このやり方でマージしました。 Rで複数のシェープファイルを結合する - Rプログラミングの小ネタ library ( maptools ) library ( ggplot2 )   # シェープファイルを読み込む pj <- CRS ( "+proj=longlat +datum=WGS84" ) shp <-r...

地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ 統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位（東区、博多区、中央区、南区、西区、城南区、早良区）の7つのシェープに分かれています。 こういうのをまとめて１つにしたい時ってありますよね。 QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。 ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。 でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この２つを結合させようとする際にエラーになってしまいます。 なので、東区のポリゴンはの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...

郵便番号がどんな感じで割り当てられているか、地図上で確認したかったんです。 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)    ...

↓こちらで公開したKMLファイルは、シェープファイルをRで変換処理したものなのですが、その手順です。 神奈川県の土砂災害危険箇所ハザードマップ（Google Earth閲覧用KML形式）: 主張 元にしたシェープは↓こちらで入手しました。 国土数値情報ダウンロードサービス ＜災害・防災＞のカテゴリに、「土砂災害危険箇所」というのがあります。シェープファイルは都道府県別に分かれて公開されています。 データ形式のタブ（？）で、「JPGIS2.1」のところを選択すると、JPGIS2.1形式とともになぜかShapeファイルも置いてあるという、微妙なトラップがあるので要注意です。 このShapeファイルをRでKMLに変換するんですが、スクリプトはごく簡単です↓ # 必要なパッケージの読み込み library ( maptools ) library ( plotKML )   # 座標系の情報を持つオブジェクトを作る # longlatは緯度経度指定、WGS84は世界測地系であることを表す pj <- CRS ( "+proj=longlat +datum=WGS84" )   # シェープファイルの読み込み hazard.area <-readShapePoly ( "xxx.shp" , proj4string=pj )   # もし、表示して確認したいときはplotすればOK # plot(hazard.area)   # KMLファイルとして出力する kml ( hazard.area , file = "xxx.kml" ) あとは、GoogleEarthで開けば閲覧できます。 ↓こんな感じですね。 国交省が公開している「土砂災害危険箇所」シェープをRでKMLに変換後、Google Earthで表示したもの 閲覧時の色の変更などは、 冒頭のリンク先 を参照してみてください。■ [2014年11月23日追記] 東京都、埼玉県、千葉県のKMLファイルも追加しました。あと、Google Earthの設定で色を変更しなくて済むよう、KMLファイル内のタグで色を赤に変えときました↓ 土砂災害危険箇所のKMLファ...

表題は正確に言うと、shapeファイルを読み込んで、加工やいろいろな可視化を施した後、画像(的なもの)として出力するときは、PDFにすると便利ですよ、という感じの趣旨です。 shapeファイルを読み込んで、plotする方法は↓以前に書きました。 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ ただ、shapeファイルってものによっては、ものすごく大量のポリゴンを含んでいますよね。地名などを重畳した場合、文字サイズが大きいと互いに重なってしまう、文字サイズが小さいと読めなくなってしまう、なんて問題が起こります。 前述のエントリでは、高解像度のpngに出力するという対処法を書いたんですけど、高解像度とはいえラスタである以上、拡大すれば文字はギザギザになるし、解像度を大きくするのにも限度がある。 使っている環境（メモリとか？ CPUのビット数とか？）によるんでしょうけど、あまり大きくすると、 > png("output.png", width= 10000 , height= 10000 )  以下にエラー png("output.png", width = 10000, height = 10000) :    デバイス png() を開始できませんでした  追加情報:  警告メッセージ: 1: In png("output.png", width = 10000, height = 10000) :    ビットマップを割り当てられません 2: In png("output.png", width = 10000, height = 10000) :   opening device failed みたいな感じにエラーになってしまいます。 で、結論としてはPDFに出力するといい感じじゃん、ってことになりました。 コードは↓こんな感じ。 # GIS関係のパッケージ library(maptools) # 座標系の情報を持つオブジェクトを作る # longlatは緯度経度指定、WGS84は世界測地系であることを表す pj <- CRS("+proj=longlat +datu...

まずはShapeファイルを入手するところから。 総務省統計局がe-Statというサイトを運営していて、そこに全国各地域のShapeファイルが公開されています。 地図で見る統計（統計ＧＩＳ) (1) 「平成２２年国勢調査（小地域） 2010/10/01」を選択 (2) 「男女別人口総数及び世帯総数」にチェック (3) 「統計表各種データダウンロードへ」ボタンをクリック (1) 「都道府県」を選択 (2) 「市区町村」を選択 (3) 「検索」ボタンをクリック (4) 境界データの欄に表示されている     「世界測地系緯度経度・Shape形式」のデータをクリック ダウンロードしたzipファイルを展開すると、 　xxx.shp 　xxx.dbf 　xxx.shx 　xxx.prj という、拡張子だけが違うファイルが4つできるんですが、これがshapeデータの1セットになります。 それぞれのファイルの中身は、 シェープファイルとは？ | 用語集とGISの使い方 | 株式会社パスコ によると、 　shp：    図形の座標が保存 　dbf：    属性の情報が保存 　shx：    shpの図形とdbfの属性の対応関係が保存 だそうです。あと、「.prj」は測地系の情報が入っています。（テキストとして開くとそれっぽい文言が入っている） では単純にshapeファイルをRで表示させてみましょう。 境界情報だけをプロット # GIS関係のパッケージ library ( maptools )   # 座標系の情報を持つオブジェクトを作る # longlatは緯度経度指定、WGS84は世界測地系であることを表す pj <- CRS ( "+proj=longlat +datum=WGS84" )   # シェープファイルの読み込み x <- readShapePoly ( "h22ka14137.shp" , proj4string=pj )   # 表示 plot ( x ) あら簡単。実質4行で表示できまし...