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ファイルを入手するところから。 総務省統計局が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行で表示できまし...