Rmapを使った地図表示

 RmapはR言語で地図表示をおこなうパッケージで、CRANの外部のプロジェクトで開発が行われている Rmap というパッケージを使い方を説明していきます*1。ここでは以下のサイトで試験的に公開されているRmapのWindowsバイナリ*2を使用します。

  http://spatial.nhh.no/R/Devel/

 Rmapは多数の地図フォーマットをサポートし*3、地図投影変換もサポートしています。

 ここでは、日本でも普及しているESRI シェープファイル*4を使った例を説明していくことにします*5

Windows でのインストール

 バイナリ版は zipファイルをインストールして、shapelib.dll と proj.dll との2つのファイルをhttp://spatial.nhh.no/R/Devel/ より入手して、R のインストールされたディレクトリのbinディレクトリに入れる。

パッケージのロード

 Rmapインストール後、パッケージをロードします。

library(rmap)

シェープファイルの所在

  • 日本
    • 統計GISプラザ  シェープファイル形式の町丁字界とCSV形式の人口・世帯データ・事業所データ(無償)。
    • 地球地図日本 行政界・水系・鉄道・道路・空港・人口集中地区 地名等は英語

変換ツール

数値地図データ変換ツール(試用版)

50mメッシュ標高・250mメッシュ標高・1kmメッシュ標高/平均標高・ 10,000総合・25,000行政界・海岸線をシェープファイルに変換

 http://www.esrij.com/products/smap/index.shtml

数値地図2500変換ツール(sdf2ogr)

国土地理院::数値地図(空間データ基盤)の閲覧(試験公開)*6より無償でまたは日本地図センターや書店で有償で入手可能な数値地図2500を MapInfo? の tab ファイルに変換。シェープファイルへの変換についても記述がある。

sdf2ogr

属性データ

 GISで利用可能なRパッケージ内データ

シェープファイルのロード

 シェープファイルの読み込みは以下のようになります。

 mySF <- shapefile("シェープファイル名(例:columbus.shp")

 ここでは、例として以下のサイトより無償で入手できる日本の行政界のシェープファイルを使用することにします。

  http://www.esrij.com/gis_data/japanshp/japanshp.html

 上のサイトより、ファイルをダウンロード・解凍後、以下のようにして、Rに読み込ます。

 jpn <- shapefile("ダウンロード先/japan_jdg.shp")

 シェープファイル読み込み時に以下のようなメッセージが出力されます。

Warning message: 
Invalid names in DBF file, new names are:
P_NUM is now P.NUM
H_NUM is now H.NUM 
in: dbfproxy(dbf) 

 上で、DBaseのフィールドのうち、フィールド名に"_"を含んでいるものは、Rの命名規則に従って、"_"が"."に変更されます。

地図の表示

 表示の方法は以下のとおりです。

 plot(jpn)

 以下の図1に表示の例を示します。

jpnmap.jpg
 図1 地図の表示例

 単一色で塗りつぶすには、以下のように(赤色)します。

 plot(jpn, f="red")

 以下の図2に表示の例を示します。

jpnmapInRed.jpg
 図2 地図の表示例(赤で塗りつぶした場合)

エクステント

bnd<-defaultextent(jpn)
xmin<-bnd[1,1]
ymin<-bnd[1,2]
xmax<-bnd[2,1]
ymax<-bnd[2,2]

属性の表示

属性のフィールドの表示は以下のようにします。
names(jpn)

 上のコマンドの実行結果は以下のとおりです。

[1] "PREF"  "CITY1" "CITY2" "TOWN1" "TOWN2" "JCODE" "P.NUM" "H.NUM" "FLAG1"
[10] "FLAG2"

 属性情報を一覧表示するには、以下のようにします。

as.data.frame(jpn)

 属性情報の一覧表示は以下のようになります*7

     PREF      CITY1 CITY2  TOWN1          TOWN2 JCODE P.NUM H.NUM FLAG1 FLAG2
1  北海道     稚内市              稚内市      1214 45754 17638     1     0
2  北海道     礼文町       礼文郡 礼文町      1517  4375  1707     1     0
3  北海道     豊富町       天塩郡 豊富町      1516  5504  1915     1     0
4  北海道 利尻富士町       利尻郡 利尻富士町    1519  4398  1526     1     0
5  北海道     利尻町       利尻郡 利尻町      1518  4104  1403     1     0
6  北海道     幌延町       天塩郡 幌延町      1488  3095  1141     1     0
7  北海道     枝幸町       枝幸郡 枝幸町      1514  8428  3021     1     0
8  北海道   中頓別町       枝幸郡 中頓別町     1513  2754   997     1     0
9  北海道     歌登町       枝幸郡 歌登町      1515  2716  1076     1     0
10 北海道     天塩町       天塩郡 天塩町      1487  4931  1910     1     0

 都道府県名のフィールド PREF を使って、北海道だけを抽出するには以下のようにします。

hokkaido<-jpn[jpn$PREF == "北海道"]

 表示は図3のようになります。

hokkaidomap.jpg
 図3 北海道の表示

 南関東のみを扱いたければ次のようになります。

 minamiKanto<-jpn[jpn$PREF == "埼玉県" | jpn$PREF == "千葉県" | jpn$PREF == "東京都"  | jpn$PREF == "神奈川県" ]

 表示は図4のようになります。

minamikantomap.jpg
 図4 南関東の表示

DBase のみのロード

db <- dbfproxy("("ダウンロード先/japan_jdg.dbf")

 属性値の地理扱いは前の節を参照のこと

Excelなどの他の属性データと地図の関係付け

match 関数を利用

地図投影変換の方法

 Rmap は100以上の地図投影法をサポートしています。  上の日本の行政界の地図は、座標は緯度経度で、測地系はJGD2000データム(測量成果2000)、回転体楕円体はGRS80です。

 日本全図を表現するのに使用されている正距円錐図法(2標準緯線、中心経度:東経135度、中心緯度:北緯35度、第1標準緯線:北緯46度、第2標準緯線:北緯20度、距離の単位:メートル)で、行政界図を表示するには以下のようにします。

plot(jpn,proj="+proj=eqdc +ellps=GRS80 +units=m +lat_0=35 +lat_1=46 +lat_2=20 +lon_0=135")

 投影変換結果をオブジェクトにするには

 現状では投影変換の出力をシェープファイルなどの外部形式には変換できません*8。 

 上の"proj"のパラメータの指定は以下のとおりです*9

  • proj=地図投影法
  • ellps=回転楕円体
  • units=距離の単位
  • lon_0=中心経度
  • lat_0=中心緯度
  • lat_1=第1標準緯線
  • lat_2=第2標準緯線

 表示は図5のようになります。

jpnmap_eqdc.jpg
 図5 正距円錐図法による投影変換

拡大表示

splancs の zoom 関数を利用して、表示したい領域の対角線上にある任意の2点を指定後、plotを使わず地図を再描画

コロプレス図(色分け主題図)の表示方法*10

 コロプレス図で連続量を扱う場合、以下の表示方法があります。ここでは、人口数データ(P.NUM)を使って説明していきます。

ユーザ定義の階級区分

hokkaidoPop <- hokkaido$P.NUM
classes <- c(5000,10000,50000,100000,500000,1000000)
breaks <- length(classes) + 1
classNo <- hokkaidoPop < classes[1]
for(i in 1:(breaks - 1))
  classNo <- classNo + (hokkaidoPop >= classes[i] & hokkaidoPop < classes[i + 1]) * i

classNo <- classNo + (hokkaidoPop >= max(classes)) * breaks
plot(hokkaido, f=gray(1 - classNo/breaks)) #地図を描く
lines(hokkaido, col="red") #アウトラインを描く
userDefinedHokkaido1.png
 図X ユーザ定義の階級区分(7階級)

等間隔

rangePop <- range(hokkaido$P.NUM)
breaks<-6 # 階級数
classNo <- trunc((hokkaido$P.NUM -rangePop[1])/(diff(rangePop)/breaks))
classNo[classNo == breaks] <- breaks - 1
hist(classNo)
plot(hokkaido,f=gray(1 - classNo/breaks)) #地図を描く
lines(hokkaido, col="red") #アウトラインを描く
equalsizehokkaido.png
 図X 等間隔(6階級)

等サイズ

breaks<-6  # 階級数
classNo <- trunc(rank(hokkaido$P.NUM)/(dim(hokkaido)[1] / breaks))
classNo[classNo == breaks] <- breaks - 1
hist(classNo)
plot(hokkaido,f=gray(1 - classNo/breaks)) #地図を描く
lines(hokkaido, col="red") #アウトラインを描く
quantilehokkaido.png
 図X 等サイズ(6階級)

標準化

sd <- 1
pal <- cm.colors(8 / sd)
classNo <- trunc(scale(hokkaido$P.NUM) + 4 / sd) + 1
classNo[classNo <= 0] <- 1
classNo[classNo > (8 / sd)] <- (8 / sd)
hist(classNo)
plot(hokkaido,f=pal[classNo]) #地図を描く
lines(hokkaido, col="red") #アウトラインを描く
standardhokkaido.png
 図X 標準化

自然階級分類(Jenk's Natural Break)

空間的コンテクストを考慮した階級区分

未分類

rangePop=range(hokkaido$P.NUM)
plot(hokkaido,f=gray(1 - (hokkaido$P.NUM - rangePop[1])/diff(rangePop))) #地図を描く
lines(hokkaido, col="red") #アウトラインを描く
unclassedhokkaido.png
 図X 未分類

凡例表示

legend関数

layout関数

階級分類の評価

 階級内の差が小さく、階級間の差を大きくする。  評価関数

二変量のコロプレス図

ドットマップ

library(splancs) #splancs パッケージのロード

saitama<-jpn[jpn$PREF == "埼玉県"]

saitama.coords <- polycoords(saitama) #座標取得
np<-length(saitama.coords) #ポリゴンの数

plot(saitama, col="red")

for(i in 1:np)
   points(csr(saitama.coords[[i]],as.integer(saitama$P_NUM/1000)),
          pch=19,cex=0.1)
saitama_dotmap.png
図 埼玉県の人口のドットマップ(1点千人)

テキストラベルの表示

 属性のラベル表示は以下のようにおこないます(ポリゴンの重心に表示する場合)。

 まず、ポリゴンの重心を求める関数(Rmapの作者 Barry Rowlingson 氏に提供していただきました)を定義します。

centroid <- function(poly) {
 np<-dim(poly)[1]
 poly<-rbind(poly,poly[1,])
# Initialize variables
 area <- 0.0
 xbara <- 0.0
 ybara <- 0.0
 for(i in  1:np){
  area <- area - ( (poly[i+1,2]-poly[i,2])*(poly[i+1,1]+poly[i,1]))/2.
  xbara <- xbara - ( ( (poly[i+1,2]-poly[i,2])/8.)*
                     ((poly[i+1,1]+poly[i,1])**2 + ((poly[i+1,1]-poly[i,1])**2)/3.))
  ybara <- ybara + ( ( (poly[i+1,1]-poly[i,1])/8.)*
                     ( (poly[i+1,2]+poly[i,2])**2 + ( (poly[i+1,2]-poly[i,2])**2)/3.))
 }
 c( xbara/area, ybara/area)
}

 次にテキスト・ラベルを表示します。都市名を表示するので、Windows版日本語グラフパッチを使用しました。

hokkaidocity <- hokkaido[hokkaido$JCODE < 1300] #市のみ選択
citycoords <- polycoords(hokkaidocity)          #各ポリゴンの座標を抽出
xy<-sapply(citycoords, centroid)                #各ポリゴンの重心計算
plot(hokkaido, col="grey")                      #行政界はグレーで表示
text(xy[1,],xy[2,], labels=hokkaidocity$CITY1, col="red", cex=0.7) #都市名の表示
textlabelshokkaido.png
 図X 北海道の都市名の表示

 テキスト・ラベルの周囲を矩形で囲むには legend を使用します。

for(i in 1:length(hokkaidocity$CITY1)) 
   legend(xy[1,i],xy[2,i], legend=hokkaidocity$CITY1[i], bg="yellow",cex=0.7) 
hoklabellegend.png
 図X 北海道の都市名を矩形で囲んで表示

円ドット図

coords <- polycoords(hokkaido)
xy<-sapply(coords, centroid)
plot(hokkaido, col="grey")
points(xy[1,],xy[2,], col="blue",cex=2*(hokkaido$P.NUM/max(hokkaido$P.NUM)),pch=19) #円ドット図表示
circledothokkaido.png
 図X 円ドット図

円の大きさを決めるのに対数を利用する

等値線図作成

等間隔

  • contour 表示
  • countourLines 座標値を返す

不等間隔

  • akima の intrerp で補間して、contour
  • TIN(ドロネイ)を使ったコンター

フローデータ

ポリゴンのオーバーレイ解析

 gpclib を使用

Union(ポリゴンの合成)

library(gpclib)
osaka<-jpn[jpn$PREF == "大阪府"]
plot(osaka)
osaka_original.png
図XXX 大阪府の各市区町村ポリゴン
osaka.coords <- polycoords(osaka) #座標取得
np<-length(osaka.coords) #ポリゴンの数

#gpclib の座標形式に変換
osaka.union = as(osaka.coords[[1]],"gpc.poly") #gpclib のオブジェクトに変換

for(i in 2:np){
  poly <- as(osaka.coords[[i]],"gpc.poly") #gpclib のオブジェクトに変換
  osaka.union<-union( osaka.union, poly) #Union
}

plot(osaka.union)
osaka_union.png
図XXX 大阪府のUnion結果(各市区町村ポリゴンをひとつのポリゴンにした)

Dissolve(任意の属性値ごとにポリゴンをグルーピングして各グループ内のポリゴンを合成)

osaka.hist <- table(as.integer(osaka$JCODE/100)) #大阪市(行政区),大阪市以外の市,郡部

osaka.dissolve <- NULL

for(i in 1:length(osaka.hist)){
  coords <- polycoords(osaka[as.integer(osaka$JCODE/100) == as.integer(names(osaka.hist)[i])]) #座標取得
  np<-length(coords) #ポリゴンの数

  osaka.part = as(coords[[1]],"gpc.poly") #gpclib のオブジェクトに変換

  for(j in 2:np){
    poly <- as(coords[[j]],"gpc.poly") #gpclib のオブジェクトに変換
    osaka.part<-union( osaka.part, poly) #Union
  }

  if(i == 1)
    osaka.dissolve <- osaka.part
  else
    osaka.dissolve <- append.poly(osaka.dissolve,osaka.part) #ポリゴンの append
}

plot(osaka.dissolve)
osaka_dissolve.png
図XXX 大阪府の Dissolve 結果(大阪市(行政区),大阪市以外の市,郡部別)

ポリゴンクリッピング

 上記の union によって作成された大阪府のポリゴンより、任意の市のポリゴンをくり抜く

matsubara  <- as(polycoords(osaka[osaka$CITY1 == "松原市"])[[1]], "gpc.poly")
plot(setdiff(osaka.union, matsubara))
osaka_clipping.png
図XXX 大阪府のポリゴンクリッピング結果(松原市のポリゴンを大阪府ポリゴンよりクリップ)

バッファー

ポイントのバッファー

ポイントの座標を中心にして任意の半径の円を描画
作成された各円をユニオン

カルトグラム

非連続型面カルトグラム

columbus<-shapefile("columbus.shp") # シェープファイルの読み込み

#
# 等サイズのコロプレス図
#
breaks<-6  # 階級数
classNo <- trunc(rank(columbus$CRIME)/(dim(columbus)[1] / breaks))
classNo[classNo == breaks] <- breaks - 1
plot(columbus,f=gray(1 - classNo/breaks)) #地図を描く
lines(columbus, col="red") #アウトラインを描く

#
# カルトグラム計算
#
coords<-polycoords(columbus) #ポリゴンの座標抽出 
densities<- sqrt(columbus$CRIME/columbus$AREA1) #密度(任意の属性値/ポリゴンの面積) 
maxDensity<-max(densities) #密度の最大値 
k <- 1/maxDensity
L<-k * densities 
centroidXY<-sapply(coords, centroid) #重心

for(i in 1:length(coords)){
  newX<-L[i] * (coords[[i]][,1] -  centroidXY[1,i]) + centroidXY[1,i] #新X座標
  newY<-L[i] * (coords[[i]][,2] -  centroidXY[2,i]) + centroidXY[2,i] #新Y座標 
  polygon(newX,newY,col="red") #ポリゴンの描画
}
noncontiguous_area_cartogram.png
図XXX 非連続型面カルトグラムの例(赤く塗りつぶされた部分)

maptools とのデータ交換

 maptools の新版はシェープファイルの書き込みをサポートしているため、Rmap でオーバーレイ解析などをおこなったデータを maptools へ送ることができれば、結果をシェープで保存できる。

ポイント

ポリゴン

地図表示の問題点

  • 道路データなどデータ量が多い場合、表示が遅くなる。
  • やはり、R では日本語が一部しか表示されないこと。現状では、これが日本での普及の最大の障害の一つか?と書きましたが、Windows版日本語グラフパッチのおかげで障害が一つなくなりました。

Rmap 関連情報サイト

 Vine Linuxでのインストールと実行

 


*1 Rmapは今後の開発の進展の見込みがないので、maptoolsというパッケージの方がいいようだ。また、S言語のmapproject 関数をRに移植したmapproj ライブラリーもある(最近、Windowsでも利用できるようになった)。
*2 Rの現行バージョンでは利用できません。現行バージョンでは、空間データのインポート機能(シェープファイル以外のMaoInfo?などの形式)は rgdal および これとリンケージのある sp パッケージを組み合わせることで、地理情報処理ことができます。
*3 Windowsバイナリ版は今のところ、シェープファイルのみです
*4 シェープファイルの日本語情報はhttp://www.esrij.com/gis_data/shape/shape.htmlにあります
*5 CRANにはその名もずばり  shapefilesというパッケージがあり、シェープファイルの読み込みばかりでなく、出力もサポートしています
*6 東京都は休止中
*7 今回使用した属性情報の日本語コードはシフトJISですが、EUC-JP などの他のコードに変換するには、dbf2txtなどでテキストデータに変換し、nkfなどでコード変換をおこなった後に、txt2dbf などでDBaseファイルに戻します
*8 オープンソースであれば、Rmap で利用している shapelib の shpproj を利用できる
*9 パラメータの詳細は、projのホームページ のドキュメントまたは http://www.remotesensing.org/geotiff/proj_list/ を参照してください
*10 RColoBrewerを使えば、論理的なカラー表現ができそう。日本語の地図学の教科書は、欧米と異なり、非常に数が少なく、過去20年も、まったく新刊がない。これらの数少ない地図学のテキストは、カラーページも少ないか (カラーのページがあるのは翻訳書)、白黒のみ (日本人の著者のものはすべて)で、対象とするのはアナログの紙地図を扱ったもので、デジタル表現に対応していない。GISの色表示の原理も、日本の大学ではあまり教えていないのでは?

添付ファイル: fileunclassedhokkaido.png 2679件 [詳細] filehoklabellegend.png 2673件 [詳細] fileequalsizehokkaido.png 2394件 [詳細] fileosaka_union.png 2582件 [詳細] fileosaka_clipping.png 2523件 [詳細] filenoncontiguous_area_cartogram.png 2661件 [詳細] fileminamikantomap.jpg 2775件 [詳細] filesaitama_dotmap.png 3143件 [詳細] fileosaka_original.png 2982件 [詳細] filejpnmap_eqdc.jpg 2850件 [詳細] fileosaka_dissolve.png 2527件 [詳細] fileuserDefinedHokkaido1.png 2875件 [詳細] filehokkaidomap.jpg 3152件 [詳細] filecircledothokkaido.png 2958件 [詳細] filestandardhokkaido.png 2733件 [詳細] filejpnmapInRed.jpg 3198件 [詳細] filequantilehokkaido.png 2780件 [詳細] filetextlabelshokkaido.png 2702件 [詳細] filejpnmap.jpg 3775件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1381d)