Comments
Description
Transcript
牧山 文彦
[連載]フリーソフトによるデータ解析・マイニング 第52回 データキューブ㈱医療情報システム部主任 牧山 文彦(Makiyama Fumihiko) ■琉球大学医学部大学院保健学研究科修了。女子栄養大学大 学院保健学研究科修了、保健学博士(学術)。琉球大学医学 部保健管理学助手、秀和綜合病院企画・電算室学術情報主任、 ちばなクリニック健康管理センター事務チーフを経て、2007 年4月より現職。E-mail:[email protected] 1. maptoolsのクラス・オブジェクト DataFrame、SpatialLinesDataFrame、Spatial Rにおいて提供されている種々の空間統計 PointsDataFrameクラスも利用できるようにな 学的手法を理解する上での一番の早道は、各 パッケージ(ライブラリ)に用いられている 「空間統計用のクラス」の構造を知ることであ る。 った。 maptoolsはspのクラスを取り込むことで、 従来のS3クラスではできなかったより細かな 処理が可能となった。 現在Rには、S3クラスとS4クラスが用いら れており[1][2]、この2種類のクラスより作り 2. Mapクラス(S3) 出されるオブジェクトを多くの解析に利用す 表1にMapクラスの構造を示した。 る。 Mapオブジェクトは、read.shape( )関数に オブジェクトとは、色々な種類の変数の塊 よってShapeFileをRに取り込むことにより生 のことで、クラスとはオブジェクトを作り出 成される。 すための雛形のことを言う。 要がある。クラス構造は、生成されたオブジ >library(maptools) 要求されたパッケージforeignをロード中です 要求されたパッケージspをロード中です >xmap<-read.shape("c:/GISdata/okinawa.shp") Shapefile type: Polygon, (5), # of Shapes: 84 ェクトに対してstr( )関数を用いると調べるこ これで、Mapオブジェクトとしてxmapが生 通常、クラス構造はマニュアルには明示さ れていないものが多いため、自分で調べる必 とができる。 成された。 さて、maptoolsでは、従来よりS3のMap、 表1を見てもわかるように、Mapオブジェ polylistクラスが用いられてきたが、spパッケー クトは2階層の構造を持っており、第1階層 ジと連携することにより、S4のSpatialPolygons にはポリゴン関連のデータが保管された 62 ● ESTRELA 2007年11月 (No.164) 表1 S3クラス:Map class type class name S3 Map obj attribute1 attribute2 attribute3 * $Shapes [[n]] $Pstart $vert $shp.type $nVerts $nParts $bbox $shpID attr(,"nVert") attr(,"nParts") attr(,"shp.type") attr(,bbox") $att.data $PREF $CITY1 $CITY2 $TOWN1 $TOWN2 $JCODE $P_NUM $H_NUM $FLAG1 $FLAG2 $S_TKY2JGD attr(*,"data_types") attr(*,"class") atomic type comment1 list list ポリゴン開始点 int ポリゴン点列 num シェープタイプ int $verts行数 int $vertsの分割数 int 最小包含箱 num シェープID int ポリゴン点の数 int ポリゴン分割の数 int シェープタイプ int 最小包含箱 num data.frame 都道府県名 Factor 市町村名 Factor 政令指定都市の市名 Factor 郡名(町村部のみ) Factor 県市区町村名 Factor 市町村コード num 人口 num 世帯数 num 陸/水フラグ num 更新時のデータ変更フラグ num 自動作成データ num $*のデータタイプ chr クラスタイプ chr comment2 ポリゴン分割に利用(複数) 緯度・経度情報 5= ポリゴン(shapefile仕様書) ポリゴン点行数 ポリゴン分割数 バウンディングボックス DBFの並びと一致 vertsがいくつに分割されているか 5= ポリゴン(shapefile仕様書) バウンディングボックス DBFのデータ 政令指定都市の場合は区名 5桁のJISコード 陸部:1、 湖沼:2 変更無:0、ver3.0:1、ver5.0:3 "C" "N" Map $Shapesと、DBFデータが保管された$att.data 表2 Mapオブジェクトからポリゴンデータを取り出す が存在している。 > x$Shapes[[1]] $Pstart [1] 0 33 *:$vertsの開始点 次に、実際に生成されたMapオブジェクト から1番目の市町村のポリゴンデータを取り 出してみる(表2)。 細かく解説すると、$Pstartには「0 33」が 入っており、この行政界は2つに分離してい ることがわかる。これは、その下に続く $vertsのポリゴン点には開始点が2つあるこ とを示している(行番号−1で表現されている ので注意すること)。 このポリゴンは、沖縄県の久米島町のもの で、Google Earthなどで確認すると実際に久 米島と鳥島の2つに分かれている。更に表3 に示すようにDBFデータを確認すると、1番 目の市町村が久米島町になっていることが確 認できる。 3. polylistクラス(S3) 表4にpolylistクラスの構造を示した。 polylistオブジェクトは、Mapオブジェクト $verts [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] [11,] [12,] [13,] [14,] [15,] [16,] [17,] [18,] [19,] [20,] [21,] [22,] [23,] [24,] [25,] [26,] [27,] [28,] [29,] [30,] [31,] [32,] [33,] [34,] [35,] [36,] [37,] [38,] [39,] [40,] [41,] [,1] 126.7568 126.7709 126.7810 126.7977 126.8056 126.8108 126.8085 126.8117 126.8173 126.8209 126.8247 126.8157 126.8159 126.8053 126.8040 126.8101 126.8123 126.8012 126.7788 126.7760 126.7680 126.7702 126.7651 126.7554 126.7588 126.7532 126.7460 126.7249 126.7152 126.7091 126.7144 126.7332 126.7568 128.2345 128.2230 128.2176 128.2218 128.2267 128.2326 128.2365 128.2345 [,2] 26.38749 *0 26.39127 26.39023 26.36874 26.37015 26.36328 26.36122 26.35387 26.35060 26.35188 26.34519 26.34134 26.33567 26.32486 26.30762 26.30232 26.29171 26.29267 26.30636 26.31445 26.32224 26.32624 26.33291 26.33884 26.34143 26.34042 26.34605 26.34787 26.35532 26.36493 26.37341 26.37589 26.38749 27.86055 *33 27.86705 27.87755 27.88536 27.88356 27.87410 27.86551 27.86055 $shp.type [1] 5 $nVerts [1] 41 $nParts [1] 2 $bbox 27.88536 $shpID [1] 0 attr(,"nVerts") [1] 41 attr(,"nParts") [1] 2 attr(,"shp.type") [1] 5 attr(,"bbox") 27.88536 2007年11月 (No.164) ESTRELA ● 63 表3 MapオブジェクトからDBFデータの取り出し けしてみる。ここでは、沖縄本島を中心に ヒートカラーを16階調にして塗りわけして > x$att.data PREF CITY1 CITY2 TOWN1 TOWN2 JCODE P_NUM H_NUM FLAG1 FLAG2 S_TKY2JGD 1 沖縄県 久米島町 <NA> 島尻郡 久米島町 47361 9359 3177 1 1 33 2 沖縄県 伊平屋村 <NA> 島尻郡 伊平屋村 47359 1434 488 1 0 0 3 沖縄県 伊平屋村 <NA> 島尻郡 伊平屋村 47359 1434 488 1 0 0 4 沖縄県 伊是名村 <NA> 島尻郡 伊是名村 47360 1895 726 1 0 0 5 沖縄県 伊是名村 <NA> 島尻郡 伊是名村 47360 1895 726 1 0 0 6 沖縄県 国頭村 <NA> 国頭郡 国頭村 47301 6015 2090 1 0 32 7 沖縄県 伊江村 <NA> 国頭郡 伊江村 47315 5131 1766 1 0 0 8 沖縄県 大宜味村 <NA> 国頭郡 大宜味村 47302 3437 1223 1 0 0 みた(ただし、塗りわけた色に特に意味は ない)。 >plot(xpoly,xlim=c(127,128.5),ylim=c(26,27), col=heat.colors(16)) <中略> 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 75 沖縄県 与那国町 <NA> 八重山郡 与那国町 47382 1801 674 1 0 0 76 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 32 77 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 78 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 79 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 80 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 81 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 82 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 83 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 84 沖縄県 竹富町 <NA> 八重山郡 竹富町 47381 3508 1526 1 0 0 74 沖縄県 0 これで、沖縄本島の塗りわけ地図が表示 された(図1)。また、付加されたデータの 中に、attr(*,"ringDir") の項目があるが、こ れはポリゴン点の描画方向(時計回り、反 時計回り)を示している。 図1 沖縄県の塗りわけ地図(polylist) からMap2poly( )関数によって生成され、Map オブジェクトからポリゴン関連のデータを取 り出し、更に詳細なデータを付加したもので ある。ポリゴン内の色の塗りわけなどの細か い処理を行うためには、必ずMapオブジェク トからの変換が必要になる。 >xpoly<-Map2poly(xmap) これで、polylistオブジェクトとしてxpolyが 生成された。更に、ポリゴン内を色で塗りわ 表4 S3クラス:polylist class type class name obj S3 polylist * attribute1 attribute2 attribute3 [[n]] [[n,m]] attr(*, "nParts") attr(*, "pstart") $from $to attr(*, "after") attr(*, "plotOrder") attr(*, "bbox") attr(*, "nParts") attr(*, "ringDir") attr(*, "area") attr(*, "centroid") $x $y attr(*, "shpID") attr(*, "region.id") attr(*, "class") attr(*, "maplim") $x $y attr(*, "after") attr(*, "plotOrder") 64 ● ESTRELA 2007年11月 (No.164) atomic type list list num int list int int num num num int int num list num num int chr chr list num num int int comment1 ポリゴン点 ポリゴン分割の数 ポリゴン開始点 プロット順 最小包含箱 comment2 緯度・経度情報 バウンディングボックス リングの向き(ポリゴン点の方向)時計回り'clockwise' エリア 重心点 シェープID リージョンID クラスタイプ プロット順 polylist フリーソフトによるデータ解析・マイニング 4. SpatialPolygonsDataFrame クラス(S4) これで、SpatialPolygonsDataFrameオブジ ェクトとしてxspldfが生成された。 沖縄本島の塗りわけ地図を表示するには、 表5にSpatialPolygonsDataFrameクラスの 構造を示した。 >plot(xspldf,xlim=c(127,128.5),ylim=c(26,27), col=heat.colors(16)) 前述した3つのS4クラスは、spパッケージ とする(図2) 。 と連携することによりmaptoolsで利用可能に なったクラスである。 図2 沖縄県の塗りわけ地図(SpatialPolygonsDataFrame) 表1と比較してみると、階層が1層深く作 られており、ほぼMapクラスとpolylistクラス を合わせた構造になっている。これによって、 Mapやpolylistのようなポリゴン点を分割する 仕様ではなく、個々のポリゴンとして管理す る仕様になり、よりわかりやすくなった。 SpatialPolygonsDataFrameオブジェクトは、 readShapePoly( )関数によってShapeFileをR に取り込むことにより生成される。 >xspldf<-readShapePoly("c:/GISdata/okinawa. shp") 表5 S4クラス:SpatialPolygonsDataFrame class type S4 class name obj Spatial Polygons DataFrame * attribute1 attribute2 attribute3 atomic type comment1 comment2 sp package 'SpatialPolygonsDataFrame'class @data @polygons data.frame Factor 都道府県名 Factor 市町村名 Factor 政令指定都市の市名 Factor 郡名(町村部のみ) Factor 県市区町村名 num 市町村コード num 人口 num 世帯数 num 陸/水フラグ num 更新時のデータ変更フラグ num 自動作成データ chr $*のデータタイプ list $PREF $CITY1 $CITY2 $TOWN1 $TOWN2 $JCODE $P_NUM $H_NUM $FLAG1 $FLAG2 $S_TKY2JGD attr(*,"data_types") [[n]] @Polygons attr(*, "dimnames") $ $ list num num logi int num int num chr num int num list chr chr @projargs chr [[n]] @labpt @area @hole @ringDir @coords @plotOrder @labpt @ID @area @plotOrder @bbox @proj4string 政令指定都市の場合は区名 5桁のJISコード 陸部:1、 湖沼:2 変更無:0、ver3.0:1、ver5.0:3 "C" "N" sp package 'Polygons'class sp package 'Polygon'class ラベルポイント エリア ホール ホール(穴)チェック リングの向き(ポリゴン点の方向)時計回り'clockwise' ポリゴン点 緯度・経度情報 プロット順 ラベルポイント ID DBFの並びと一致 エリア プロット順 最小包含箱 バウンディングボックス bbox配列名 bbox配列名 "r1" "r2" bbox配列名 "min" "max" sp package 'CRS'class 地図投影法指定 2007年11月 (No.164) ESTRELA ● 65 とする(図3)。 このようにS4クラスでは直接塗りわけ地図 の作成が可能となっており、更に表5を細か 図3 沖縄県の線地図(SpatialLinesDataFrame) く見ると、SpatialPolygonsDataFrameクラス の特徴として、@hole、@proj4stringが実装さ れていることがわかる。 @holeは、ポリゴンホール(穴のポリゴン) の識別、@proj4stringは地図投影法の指定とな っており、この2つを用いることで、より柔 軟で詳細な地図の描画が可能となる。 5. SpatialLinesDataFrameクラス 表6にSpatialLinesDataFrameクラスの構造 を示した。このクラスはSpatialPolygonsData Frameクラスのライン(線)版であり、主に 線データの操作に使う。 6. SpatialPointsDataFrameクラス >xsldf<-readShapeLines("c:/GISdata/okinawa. shp") 表7にSpatialPointsDataFrameクラスの構造 これで、SpatialLinesDataFrameオブジェク を示した。SpatialPointsDataFrameクラスは トとしてxsldfが生成された。更に、沖縄本島 GISdata内のShapeFileでは生成できないため、 の線地図を青で表示するには、 maptoolsで提供されているデータを利用して >plot(xsldf, xlim=c(127,128.5),ylim=c(26,27), col="blue") 説明する。 表6 S4クラス:SpatialLinesDataFrame class type S4 class name obj SpatialLines DataFrame * attribute1 attribute2 attribute3 atomic type comment1 @data @lines data.frame Factor 都道府県名 Factor 市町村名 Factor 政令指定都市の市名 Factor 郡名(町村部のみ) Factor 県市区町村名 num 市町村コード num 人口 num 世帯数 num 陸/水フラグ num 更新時のデータ変更フラグ num 自動作成データ chr $* のデータタイプ list $PREF $CITY1 $CITY2 $TOWN1 $TOWN2 $JCODE $P_NUM $H_NUM $FLAG1 $FLAG2 $S_TKY2JGD attr(* ,"data_types") [[n]] @Lines [[n]] @coords attr(* , "dimnames") $ $ list num chr num list chr chr ポリゴン点 ID 最小包含箱 bbox配列名 bbox配列名 bbox配列名 @projargs chr 地図投影法指定 @ID @bbox @proj4string 66 ● ESTRELA comment2 sp package 'SpatialLinesDataFrame'class 2007年11月 (No.164) 政令指定都市の場合は区名 5桁のJISコード 陸部:1、 湖沼:2 変更無:0、ver3.0:1、ver5.0:3 "C" "N" sp package 'Lines'class sp package 'Line'class 緯度・経度情報 バウンディングボックス "r1" "r2" "min" "max" sp package 'CRS' class フリーソフトによるデータ解析・マイニング 表7 S4クラス:SpatialPointsDataFrame class type S4 class name obj attribute1 attribute2 attribute3 atomic type comment1 SpatialPoints * DataFrame comment2 sp package 'SpatiaPointsDataFrame'class @data data.frame num num num num num num num num num num num num num num num num 経度 num 緯度 chr $*のデータタイプ num $STATION $PRICE $NROOM $DWELL $NBATH $PATIO $FIREPL $AC $BMENT $NSTOR $GAR $AGE $CITCOU $LOTSZ $SQFT $X $Y attr(*,"data_types") @coods.nrs @coods attr(*, "dimnames") $ $ @bbox attr(*, "dimnames") $ $ list NULL chr num list chr chr "coords.x1" "coords.x2" バウンディングボックス 最小包含箱 bbox配列名 bbox配列名 bbox配列名 "coords.x1" "coords.x2" "min" "max" sp package 'CRS' class "NA" @proj4string @projargs chr "N" 地図投影法指定 注)GISdata.zipで提供しているShapeFileからSpatialPointsDataFrameは生成できないので、system.file(“shapes/baltim.shp”, package= “maptools”[1])を使った。 >xsptdf<- readShapePoints(system.file ("shapes/baltim.shp", package="maptools")[1]) 7. その他の機能 以上のように、maptoolsはsp、gpclib、 これで、SpatialPointsDataFrameオブジェク shapefiles、rgdal等、多くのパッケージと連携 トとしてxsptdfが生成された。更に、ポイン することで、高度なShapeFileの処理が可能と ト地図(赤いポイント)を表示するには、 なっている。 また、現在でも新しい機能の更新が続いて >plot(xsptdf,col="red") とする(図4)。 おり、最近ではGoogle Earth上にRで生成し たPNG画像をオーバーレイする機能 (kmlOverlay( )関数)も追加された。 図4 ポイントの描画(SpatialPointsDataFrame) *参考文献・URL [1] U . リ ゲ ス : R の 基 礎 と プ ロ グ ラ ミ ン グ 技 法 , pp.121 - 133. [2] RjpWiki:S4クラスとメソッド入門(http://www. okada.jp.org/RWiki/?S4%20%A5%AF%A5%E9%A5% B9%A4%C8%A5%E1%A5%BD%A5%C3%A5%C9%C6% FE%CC%E7). 2007年11月 (No.164) ESTRELA ● 67