地図をプロットするパッケージはいろいろありますが、経緯線の描写には手間がかかる場合があります。そんな問題を解決する「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")実行コマンド
詳細はコマンド、パッケージヘルプを確認してください。近畿および名古屋圏と日本全体をプロットする例を紹介します。
#パッケージの読み込み
library("rgdal")
library("raster")
library("rworldmap")
library("graticule")
###近畿および名古屋圏のプロット#####
#地図データの読み込み
#rworldmapパッケージのデータを使用
data(countriesLow)
#rasterパッケージのprojectionコマンド
#座標を取得
llproj <- projection(countriesLow)
#地図データを抽出
map <- subset(countriesLow, SOVEREIGNT == "Japan")
#参考:地域のデータ一覧参照コマンド
countriesLow$SOVEREIGNT
#地図投影法の設定
#ランベルト等角円錐図法で指定
prj <- "+proj=lcc +lat_1=40 +lat_2=50 +lat_0=30 +lon_0=125 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
#データに適応:spTransformコマンド
pmap <- spTransform(map, CRS(prj))
###描写範囲の設定#####
#経度の設定
lons <- seq(135.02, 137.43, length = 5)
#緯度の設定
lats <- seq(33.34, 35.58, length = 6)
#####
###メッシュの遊びを指定#####
#経度の遊びを指定,カッコ内を変更
xl <- range(lons) + c(-0.04, 0.04)
#緯度の遊びを指定,カッコ内を変更
yl <- range(lats) + c(-0.04, 0.04)
#####
#メッシュプロットのデータを作成:graticuleコマンド
grat <- graticule(lons, lats, proj = prj, xlim = xl, ylim = yl)
#メッシュプロットラベルのデータを作成:graticule_labelsコマンド
labs <- graticule_labels(lons, lats, xline = min(xl), yline = max(yl), proj = prj)
#プロット領域の周囲余白調整
op <- par(mar = rep(0, 4))
#プロット領域の描写
#[c()]内の数字を大きくすると広域表示になります
plot(extent(grat) + c(0.5, 0.3) * 1e5, asp = 1, type = "n", axes = FALSE, xlab = "", ylab = "")
#地図の描写
plot(pmap, add = TRUE, col = "#00bfd4")
#メッシュの描写
plot(grat, add = TRUE, lty = 5, col = "red")
#メッシュテキストの描写
#経度側
text(subset(labs, labs$islon), lab = parse(text = labs$lab[labs$islon]), pos = 3)
#緯度側
text(subset(labs, !labs$islon), lab = parse(text = labs$lab[!labs$islon]), pos = 2)
#####
###日本全域をプロット#####
#llgridlinesコマンドとの組み合わせ
par(xpd = NA)
plot(pmap, col = "#00bfd4")
llgridlines(pmap, col = "red")
#####出力例
・近畿および名古屋圏のプロット
・日本全域をプロット
少しでも、あなたのウェブや実験の解析が楽になりますように！！
