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")

実行コマンド

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

#パッケージの読み込み
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") ##### [/code]


出力例

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

・日本全域をプロット
japanmat


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

スポンサードリンク

関連コンテンツ


スポンサードリンク