![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
//尾野久二 COLOR(red){SIZE(30){Rmapを使った地図表示}} [[Rmap:http://www.maths.lancs.ac.uk/Software/Rmap/]]はR言語で地図表示をおこなうパッケージで、CRANの外部のプロジェクトで開発が行われている [[Rmap:http://www.maths.lancs.ac.uk/Software/Rmap/]] というパッケージを使い方を説明していきます((Rmapは今後の開発の進展の見込みがないので、[[maptools:http://sal.agecon.uiuc.edu/csiss/Rgeo/index.html]]というパッケージの方がいいようだ。また、S言語のmapproject 関数をRに移植した[[mapproj :http://www.stat.cmu.edu/~minka/software/maps/]]ライブラリーもある（最近、Windowsでも利用できるようになった）。))。ここでは以下のサイトで試験的に公開されているRmapのWindowsバイナリ((Rの現行バージョンでは利用できません。現行バージョンでは、空間データのインポート機能（シェープファイル以外のMaoInfoなどの形式）は rgdal および これとリンケージのある sp パッケージを組み合わせることで、地理情報処理ことができます。))を使用します。 [[Rmap:http://www.maths.lancs.ac.uk/Software/Rmap/]]はR言語で地図表示をおこなうパッケージで、CRANの外部のプロジェクトで開発が行われている [[Rmap:http://www.maths.lancs.ac.uk/Software/Rmap/]] というパッケージを使い方を説明していきます((Rmapは今後の開発の進展の見込みがないようで、[[maptools:http://sal.agecon.uiuc.edu/csiss/Rgeo/index.html]]というパッケージの方がいいようだ。また、S言語のmapproject 関数をRに移植した[[mapproj :http://www.stat.cmu.edu/~minka/software/maps/]]ライブラリーもある（最近、Windowsでも利用できるようになった）。))。ここでは以下のサイトで試験的に公開されているRmapのWindowsバイナリ((Rの現行バージョンでは利用できません。現行バージョンでは、空間データのインポート機能（シェープファイル以外のMaoInfoなどの形式）は rgdal および これとリンケージのある sp パッケージを組み合わせることで、地理情報処理ことができます。))を使用します。 http://spatial.nhh.no/R/Devel/ Rmapは多数の地図フォーマットをサポートし((Windowsバイナリ版は今のところ、シェープファイルのみです))、地図投影変換もサポートしています。 ここでは、日本でも普及しているESRI シェープファイル((シェープファイルの日本語情報はhttp://www.esrij.com/gis_data/shape/shape.htmlにあります))を使った例を説明していくことにします((CRANにはその名もずばり [[shapefiles:http://www.odot.state.or.us/tddtpau/R.html]]というパッケージがあり、シェープファイルの読み込みばかりでなく、出力もサポートしています))。 #contents * Windows でのインストール * Windows でのインストール [#x4300e2a] バイナリ版は zipファイルをインストールして、shapelib.dll と proj.dll との２つのファイルをhttp://spatial.nhh.no/R/Devel/ より入手して、R のインストールされたディレクトリのbinディレクトリに入れる。 * パッケージのロード * パッケージのロード [#p3df5069] Rmapインストール後、パッケージをロードします。 library(rmap) * シェープファイルの所在 * シェープファイルの所在 [#o2387631] - 世界 -- [[DataServer: Free maps and GIS data :http://www.diva-gis.org/data/DataServer.htm]] 世界各国 -- [[Shapefiles for Epi Info:http://www.cdc.gov/epiinfo/shape.htm]] 各国の行政界 - アメリカ -- [[アメリカ地質調査所のナショナルアトラス（国勢地図帳）のデータセット:http://nationalatlas.gov/atlasftp.html#layerstable]] -- ESRI - Geography Network - 日本 -- [[統計GISプラザ:http://gisplaza.stat.go.jp]] シェープファイル形式の町丁字界とCSV形式の人口・世帯データ・事業所データ(無償)。 -- [[地球地図日本:http://www1.gsi.go.jp/geowww/globalmap-gsi/download/index.html]] 行政界・水系・鉄道・道路・空港・人口集中地区 地名等は英語 *変換ツール ** 数値地図データ変換ツール（試用版) *変換ツール [#sfef0a08] ** 数値地図データ変換ツール（試用版) [#c722dc36] 50mメッシュ標高・250mメッシュ標高・1kmメッシュ標高／平均標高・ 10,000総合・25,000行政界・海岸線をシェープファイルに変換 http://www.esrij.com/products/smap/index.shtml **数値地図2500変換ツール(sdf2ogr) **数値地図2500変換ツール(sdf2ogr) [#yf366326] [[国土地理院::数値地図(空間データ基盤)の閲覧（試験公開）:http://sdf.gsi.go.jp/]]((東京都は休止中))より無償でまたは日本地図センターや書店で有償で入手可能な数値地図2500を MapInfo の tab ファイルに変換。シェープファイルへの変換についても記述がある。 [[sdf2ogr:http://www.max.hi-ho.ne.jp/scream/sdf2ogr.html]] *属性データ *属性データ [#m26365b7] [[GISで利用可能なRパッケージ内データ]] * シェープファイルのロード * シェープファイルのロード [#o32131ab] シェープファイルの読み込みは以下のようになります。 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の命名規則に従って、"_"が"."に変更されます。 * 地図の表示 * 地図の表示 [#s571ab32] 表示の方法は以下のとおりです。 plot(jpn) 以下の図１に表示の例を示します。 #ref(jpnmap.jpg,center) CENTER: 図１ 地図の表示例 単一色で塗りつぶすには、以下のように（赤色）します。 plot(jpn, f="red") 以下の図２に表示の例を示します。 #ref(jpnmapInRed.jpg,center) CENTER: 図２ 地図の表示例（赤で塗りつぶした場合） * エクステント * エクステント [#t9945978] bnd<-defaultextent(jpn) xmin<-bnd[1,1] ymin<-bnd[1,2] xmax<-bnd[2,1] ymax<-bnd[2,2] * 属性の表示 * 属性の表示 [#z2f8be24] 属性のフィールドの表示は以下のようにします。 names(jpn) 上のコマンドの実行結果は以下のとおりです。 [1] "PREF" "CITY1" "CITY2" "TOWN1" "TOWN2" "JCODE" "P.NUM" "H.NUM" "FLAG1" [10] "FLAG2" 属性情報を一覧表示するには、以下のようにします。 as.data.frame(jpn) 属性情報の一覧表示は以下のようになります((今回使用した属性情報の日本語コードはシフトＪＩＳですが、EUC-JP などの他のコードに変換するには、[[dbf2txt:http://www.usf.uos.de/~fkoorman/software/dbftools.en.html]]などでテキストデータに変換し、nkfなどでコード変換をおこなった後に、[[txt2dbf:http://www.usf.uos.de/~fkoorman/software/dbftools.en.html]] などでDBaseファイルに戻します))。 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 == "北海道"] 表示は図３のようになります。 #ref(hokkaidomap.jpg,center) CENTER: 図３ 北海道の表示 南関東のみを扱いたければ次のようになります。 minamiKanto<-jpn[jpn$PREF == "埼玉県" | jpn$PREF == "千葉県" | jpn$PREF == "東京都" | jpn$PREF == "神奈川県" ] 表示は図４のようになります。 #ref(minamikantomap.jpg,center) CENTER: 図４ 南関東の表示 *DBase のみのロード *DBase のみのロード [#k728a05f] db <- dbfproxy("("ダウンロード先/japan_jdg.dbf") 属性値の地理扱いは前の節を参照のこと *Excelなどの他の属性データと地図の関係付け *Excelなどの他の属性データと地図の関係付け [#n4c43732] match 関数を利用 * 地図投影変換の方法 * 地図投影変換の方法 [#i6af3c6d] Rmap は100以上の地図投影法をサポートしています。 上の日本の行政界の地図は、座標は緯度経度で、測地系はJGD2000データム（測量成果2000）、回転体楕円体はGRS80です。 日本全図を表現するのに使用されている正距円錐図法（２標準緯線、中心経度：東経135度、中心緯度：北緯35度、第１標準緯線：北緯46度、第２標準緯線：北緯20度、距離の単位：メートル）で、行政界図を表示するには以下のようにします。 plot(jpn,proj="+proj=eqdc +ellps=GRS80 +units=m +lat_0=35 +lat_1=46 +lat_2=20 +lon_0=135") 投影変換結果をオブジェクトにするには 現状では投影変換の出力をシェープファイルなどの外部形式には変換できません((オープンソースであれば、Rmap で利用している shapelib の shpproj を利用できる))。 上の"proj"のパラメータの指定は以下のとおりです((パラメータの詳細は、[[projのホームページ:http://www.remotesensing.org/proj]] のドキュメントまたは http://www.remotesensing.org/geotiff/proj_list/ を参照してください))。 - proj=地図投影法 - ellps=回転楕円体 - units=距離の単位 - lon_0=中心経度 - lat_0=中心緯度 - lat_1=第１標準緯線 - lat_2=第２標準緯線 表示は図５のようになります。 #ref(jpnmap_eqdc.jpg,center) CENTER: 図５ 正距円錐図法による投影変換 *拡大表示 *拡大表示 [#df88eb72] splancs の zoom 関数を利用して、表示したい領域の対角線上にある任意の２点を指定後、plotを使わず地図を再描画 * コロプレス図（色分け主題図）の表示方法(([[RColoBrewer:http://cran.at.r-project.org/src/contrib/PACKAGES.html#RColorBrewer]]を使えば、論理的なカラー表現ができそう。日本語の地図学の教科書は、欧米と異なり、非常に数が少なく、過去２０年も、まったく新刊がない。これらの数少ない地図学のテキストは、カラーページも少ないか （カラーのページがあるのは翻訳書）、白黒のみ （日本人の著者のものはすべて）で、対象とするのはアナログの紙地図を扱ったもので、デジタル表現に対応していない。GISの色表示の原理も、日本の大学ではあまり教えていないのでは？)) * コロプレス図（色分け主題図）の表示方法(([[RColoBrewer:http://cran.at.r-project.org/src/contrib/PACKAGES.html#RColorBrewer]]を使えば、論理的なカラー表現ができそう。日本語の地図学の教科書は、欧米と異なり、非常に数が少なく、過去２０年も、まったく新刊がない。これらの数少ない地図学のテキストは、カラーページも少ないか （カラーのページがあるのは翻訳書）、白黒のみ （日本人の著者のものはすべて）で、対象とするのはアナログの紙地図を扱ったもので、デジタル表現に対応していない。GISの色表示の原理も、日本の大学ではあまり教えていないのでは？)) [#l19d17a9] コロプレス図で連続量を扱う場合、以下の表示方法があります。ここでは、人口数データ（P.NUM）を使って説明していきます。 ** ユーザ定義の階級区分 ** ユーザ定義の階級区分 [#tc8dd94c] 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") #アウトラインを描く #ref(userDefinedHokkaido1.png,center) CENTER: 図X ユーザ定義の階級区分(７階級) ** 等間隔 ** 等間隔 [#r639c697] 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") #アウトラインを描く #ref(equalsizehokkaido.png,center) CENTER: 図X 等間隔(６階級) ** 等サイズ ** 等サイズ [#q6b6438d] 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") #アウトラインを描く #ref(quantilehokkaido.png,center) CENTER: 図X 等サイズ(６階級) ** 標準化 ** 標準化 [#f5951fbc] 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") #アウトラインを描く #ref(standardhokkaido.png,center) CENTER: 図X 標準化 ** 自然階級分類(Jenk's Natural Break) ** 自然階級分類(Jenk's Natural Break) [#m209d952] ** 空間的コンテクストを考慮した階級区分 ** 空間的コンテクストを考慮した階級区分 [#j66aa5ff] ** 未分類 ** 未分類 [#lb02b6ba] rangePop=range(hokkaido$P.NUM) plot(hokkaido,f=gray(1 - (hokkaido$P.NUM - rangePop[1])/diff(rangePop))) #地図を描く lines(hokkaido, col="red") #アウトラインを描く #ref(unclassedhokkaido.png,center) CENTER: 図X 未分類 **凡例表示 ***legend関数 ***layout関数 **凡例表示 [#f7f3d8c3] ***legend関数 [#gc62fd70] ***layout関数 [#p5575b80] ** 階級分類の評価 ** 階級分類の評価 [#te021748] 階級内の差が小さく、階級間の差を大きくする。 評価関数 **二変量のコロプレス図 **二変量のコロプレス図 [#w777c8e8] **ドットマップ **ドットマップ [#e2b3d654] 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) #ref(saitama_dotmap.png,center) CENTER:図 埼玉県の人口のドットマップ（１点千人） *テキストラベルの表示 *テキストラベルの表示 [#h688337c] 属性のラベル表示は以下のようにおこないます（ポリゴンの重心に表示する場合）。 まず、ポリゴンの重心を求める関数(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) #都市名の表示 #ref(textlabelshokkaido.png,center) CENTER: 図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) #ref(hoklabellegend.png,center) CENTER: 図X 北海道の都市名を矩形で囲んで表示 *円ドット図 *円ドット図 [#aa0890ff] 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) #円ドット図表示 #ref(circledothokkaido.png,center) CENTER: 図X 円ドット図 円の大きさを決めるのに対数を利用する *等値線図作成 *等値線図作成 [#g8081553] **等間隔 **等間隔 [#x4238fb0] -- contour 表示 -- countourLines 座標値を返す **不等間隔 **不等間隔 [#p3a898b0] -- akima の intrerp で補間して、contour -- TIN(ドロネイ)を使ったコンター *フローデータ *フローデータ [#mcc45ada] *ポリゴンのオーバーレイ解析 *ポリゴンのオーバーレイ解析 [#kc84a31e] gpclib を使用 **Union（ポリゴンの合成） **Union（ポリゴンの合成） [#b9b85e55] library(gpclib) osaka<-jpn[jpn$PREF == "大阪府"] plot(osaka) #ref(osaka_original.png,center) CENTER:図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) #ref(osaka_union.png,center) CENTER:図XXX 大阪府のＵｎｉｏｎ結果(各市区町村ポリゴンをひとつのポリゴンにした) **Dissolve（任意の属性値ごとにポリゴンをグルーピングして各グループ内のポリゴンを合成） **Dissolve（任意の属性値ごとにポリゴンをグルーピングして各グループ内のポリゴンを合成） [#c65b7a41] 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) #ref(osaka_dissolve.png,center) CENTER:図XXX 大阪府の Dissolve 結果(大阪市(行政区),大阪市以外の市,郡部別) **ポリゴンクリッピング **ポリゴンクリッピング [#hd9a1913] 上記の union によって作成された大阪府のポリゴンより、任意の市のポリゴンをくり抜く matsubara <- as(polycoords(osaka[osaka$CITY1 == "松原市"])[[1]], "gpc.poly") plot(setdiff(osaka.union, matsubara)) #ref(osaka_clipping.png,center) CENTER:図XXX 大阪府のポリゴンクリッピング結果(松原市のポリゴンを大阪府ポリゴンよりクリップ) *バッファー **ポイントのバッファー *バッファー [#c374c493] **ポイントのバッファー [#odff9516] ポイントの座標を中心にして任意の半径の円を描画~ 作成された各円をユニオン~ *カルトグラム **非連続型面カルトグラム *カルトグラム [#m1396bf0] **非連続型面カルトグラム [#sa9d430c] 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") #ポリゴンの描画 } #ref(noncontiguous_area_cartogram.png,center) CENTER:図XXX 非連続型面カルトグラムの例(赤く塗りつぶされた部分) *maptools とのデータ交換 *maptools とのデータ交換 [#b6687ab3] maptools の新版はシェープファイルの書き込みをサポートしているため、Rmap でオーバーレイ解析などをおこなったデータを maptools へ送ることができれば、結果をシェープで保存できる。 **ポイント **ポイント [#d4c112bc] **ポリゴン **ポリゴン [#k177b8af] * 地図表示の問題点 * 地図表示の問題点 [#k6a60626] -道路データなどデータ量が多い場合、表示が遅くなる。 -やはり、R では日本語が一部しか表示されないこと。現状では、これが日本での普及の最大の障害の一つか？と書きましたが、[[Windows版日本語グラフパッチ]]のおかげで障害が一つなくなりました。 * Rmap 関連情報サイト * Rmap 関連情報サイト [#t44d1ee7] -[[なかまさんのサイト:http://www.nakama.ne.jp/memo/cran-R/Rmap/index.html]] Vine Linuxでのインストールと実行
