Rで解析：地図に経緯線の描写が楽々！「graticule」パッケージ

地図をプロットするパッケージはいろいろありますが、経緯線の描写には手間がかかる場合があります。そんな問題を解決する「graticule」パッケージです。

本パッケージは「rgdal」「raster」「rworldmap」のパッケージを利用します。なお、世界地図のデータは「rworldmap」パッケージに収録されていますので、別途用意する必要はありません。

発想次第では、地図の描写以外にも利用が広がりそうなパッケージです。公式の参考URLに示されているよう、複雑なカラーパレット作成にも利用できそうです。

公式参考URL
https://cran.r-project.org/web/packages/graticule/vignettes/graticule.html

また、経度、緯度の取得は下記サイトがオススメです。
地理院地図
http://maps.gsi.go.jp/

パッケージのバージョンは1.3-1。R version 3.2.1でコマンドを確認しています。

必要パッケージのインストール

下記コマンドを実行してください。

#パッケージのインストール
install.packages("graticule")
install.packages("rgdal")
install.packages("raster")
install.packages("rworldmap")

実行コマンド

詳細はコマンド、パッケージヘルプを確認してください。近畿および名古屋圏と日本全体をプロットする例を紹介します。

#&#12497;&#12483;&#12465;&#12540;&#12472;&#12398;&#35501;&#12415;&#36796;&#12415;
library("rgdal")
library("raster")
library("rworldmap")
library("graticule")

###&#36817;&#30079;&#12362;&#12424;&#12403;&#21517;&#21476;&#23627;&#22287;&#12398;&#12503;&#12525;&#12483;&#12488;#####
#&#22320;&#22259;&#12487;&#12540;&#12479;&#12398;&#35501;&#12415;&#36796;&#12415;
#rworldmap&#12497;&#12483;&#12465;&#12540;&#12472;&#12398;&#12487;&#12540;&#12479;&#12434;&#20351;&#29992;
data(countriesLow)
#raster&#12497;&#12483;&#12465;&#12540;&#12472;&#12398;projection&#12467;&#12510;&#12531;&#12489;
#&#24231;&#27161;&#12434;&#21462;&#24471;
llproj <- projection(countriesLow)
#&#22320;&#22259;&#12487;&#12540;&#12479;&#12434;&#25277;&#20986;
map <- subset(countriesLow, SOVEREIGNT == "Japan")
#&#21442;&#32771;:&#22320;&#22495;&#12398;&#12487;&#12540;&#12479;&#19968;&#35239;&#21442;&#29031;&#12467;&#12510;&#12531;&#12489;
countriesLow$SOVEREIGNT

#&#22320;&#22259;&#25237;&#24433;&#27861;&#12398;&#35373;&#23450;
#&#12521;&#12531;&#12505;&#12523;&#12488;&#31561;&#35282;&#20870;&#37648;&#22259;&#27861;&#12391;&#25351;&#23450;
prj <- "+proj=lcc +lat_1=40 +lat_2=50 +lat_0=30 +lon_0=125 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

#&#12487;&#12540;&#12479;&#12395;&#36969;&#24540;:spTransform&#12467;&#12510;&#12531;&#12489;
pmap <- spTransform(map, CRS(prj))

###&#25551;&#20889;&#31684;&#22258;&#12398;&#35373;&#23450;#####
#&#32076;&#24230;&#12398;&#35373;&#23450;
lons <- seq(135.02, 137.43, length = 5)
#&#32239;&#24230;&#12398;&#35373;&#23450;
lats <- seq(33.34, 35.58, length = 6)
#####

###&#12513;&#12483;&#12471;&#12517;&#12398;&#36938;&#12403;&#12434;&#25351;&#23450;#####
#&#32076;&#24230;&#12398;&#36938;&#12403;&#12434;&#25351;&#23450;,&#12459;&#12483;&#12467;&#20869;&#12434;&#22793;&#26356;
xl <- range(lons) + c(-0.04, 0.04)
#&#32239;&#24230;&#12398;&#36938;&#12403;&#12434;&#25351;&#23450;,&#12459;&#12483;&#12467;&#20869;&#12434;&#22793;&#26356;
yl <- range(lats) + c(-0.04, 0.04)
#####

#&#12513;&#12483;&#12471;&#12517;&#12503;&#12525;&#12483;&#12488;&#12398;&#12487;&#12540;&#12479;&#12434;&#20316;&#25104;:graticule&#12467;&#12510;&#12531;&#12489;
grat <- graticule(lons, lats, proj = prj, xlim = xl, ylim = yl)

#&#12513;&#12483;&#12471;&#12517;&#12503;&#12525;&#12483;&#12488;&#12521;&#12505;&#12523;&#12398;&#12487;&#12540;&#12479;&#12434;&#20316;&#25104;:graticule_labels&#12467;&#12510;&#12531;&#12489;
labs <- graticule_labels(lons, lats, xline = min(xl), yline = max(yl), proj = prj)

#&#12503;&#12525;&#12483;&#12488;&#38936;&#22495;&#12398;&#21608;&#22258;&#20313;&#30333;&#35519;&#25972;
op <- par(mar = rep(0, 4))

#&#12503;&#12525;&#12483;&#12488;&#38936;&#22495;&#12398;&#25551;&#20889;
#[c()]&#20869;&#12398;&#25968;&#23383;&#12434;&#22823;&#12365;&#12367;&#12377;&#12427;&#12392;&#24195;&#22495;&#34920;&#31034;&#12395;&#12394;&#12426;&#12414;&#12377;
plot(extent(grat) + c(0.5, 0.3) * 1e5, asp = 1, type = "n", axes = FALSE, xlab = "", ylab = "")

#&#22320;&#22259;&#12398;&#25551;&#20889;
plot(pmap, add = TRUE, col = "#00bfd4")
#&#12513;&#12483;&#12471;&#12517;&#12398;&#25551;&#20889;
plot(grat, add = TRUE, lty = 5, col = "red")
#&#12513;&#12483;&#12471;&#12517;&#12486;&#12461;&#12473;&#12488;&#12398;&#25551;&#20889;
#&#32076;&#24230;&#20596;
text(subset(labs, labs$islon), lab = parse(text = labs$lab[labs$islon]), pos = 3)
#&#32239;&#24230;&#20596;
text(subset(labs, !labs$islon), lab = parse(text = labs$lab[!labs$islon]), pos = 2)
#####

###&#26085;&#26412;&#20840;&#22495;&#12434;&#12503;&#12525;&#12483;&#12488;#####
#llgridlines&#12467;&#12510;&#12531;&#12489;&#12392;&#12398;&#32068;&#12415;&#21512;&#12431;&#12379;
par(xpd = NA)
plot(pmap, col = "#00bfd4")
llgridlines(pmap, col = "red")
#####

出力例

・近畿および名古屋圏のプロット

・日本全域をプロット

少しでも、あなたのウェブや実験の解析が楽になりますように！！