...

牧山 文彦

by user

on
Category: Documents
28

views

Report

Comments

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
Fly UP