5. 簡単な機能のプロット

Edzer Pebesma

For a better version of the sf vignettes see https://r-spatial.github.io/sf/articles/

この vignette は、シンプルフィーチャをプロットするのに役立つsfの関数について説明しています。また、他のパッケージ(mapview, tmap, ggplot2)でシンプルフィーチャを描画するための例やオプションへのポインタを提供します。

オブジェクト sf と sfc のプロット方法

ジオメトリのみ: sfc

ジオメトリのリストカラム（st_geometryメソッドで得られます sfc クラスのオブジェクト）には、ジオメトリのみが表示されます。

library(sf)
## Linking to GEOS 3.10.3, GDAL 3.5.1, PROJ 7.2.1; sf_use_s2() is TRUE
demo(nc, ask = FALSE, echo = FALSE)
plot(st_geometry(nc))

で、通常の基本プロットと同様に、色や記号などでさらに注釈を加えることができます。例えば、ポリゴンプロットにポイントを追加するには、次のようにします。

plot(st_geometry(nc), col = sf.colors(12, categorical = TRUE), border = 'grey', 
     axes = TRUE)
plot(st_geometry(st_centroid(nc)), pch = 3, col = 'red', add = TRUE)
## Warning in st_centroid.sf(nc): st_centroid assumes attributes are constant over
## geometries of x

で、凡例やタイトルなどは後から追加することができます。border=NA はポリゴンのボーダーを削除します。

見てのとおり、プロットされる軸はCRSに依存し、経度緯度座標の場合、axes = TRUE であれば度記号と方位が追加されます。

属性を持つジオメトリ sf

sf オブジェクトのデフォルトのプロットは、すべての属性のマルチプロットで、妥当な最大値までです。

plot(nc)
## Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to plot
## all

すべての属性が合理的にプロットできない場合は、警告が表示されます。プロットされるマップの最大数を増やすには

plot(nc, max.plot = 14)

行／列のレイアウトは、プロット領域が最大に満たされるように選択されます。max.plot のデフォルト値は、例えば、グローバルオプション sf_max.plot を設定することにより、制御することができます。

options(sf_max.plot=1)
plot(nc)

カラーキーの位置とサイズ

単一の属性を選択している場合、デフォルトでは、カラーキーは、プロットされているマップのためにできるだけ多くのスペースを残す側に与えられます; nc の場合、これは以下のようになります。

plot(nc["AREA"])

が、これは制御可能で、特定の側（1=下、2=左、3=上、4=右）に設定することができます。

plot(nc["AREA"], key.pos = 4)

カラーキーの大きさは、相対単位（0から1の間の数値）と絶対単位（2cmを表すlcm(2)のようなもの）のどちらかを使って制御できます。

plot(nc["AREA"], key.pos = 1, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)

因子変数のキーは少し異なり、通常、因子変数のテキストを回転させたくないからです。

nc$f = cut(nc$AREA, 10)
plot(nc["f"], axes = TRUE, key.pos = 4, pal = sf.colors(10), key.width = lcm(4.5))

クラス間隔

カラーブレーク（クラス間隔）は、プロットの引数 breaks と nbreaks によって制御することができます。nbreaks はブレークの数を指定します。 breaks はブレーク値を持つベクトルです。

plot(nc["AREA"], breaks = c(0,.05,.1,.15,.2,.25))

または breaks は、 classInt::classIntervals の style 引数として渡されるブレークファインディングメソッドを示すために使用されます。デフォルトの値である pretty は、クラスの区切りを丸くするものです。その他の方法としては、データ範囲を "nbreaks" 等しいクラスに分割する "equal" や、分位数をクラスの区切りとして利用する "quantile" 、そして他のソフトウェアで利用される "jenks" などがあります。

plot(nc["AREA"], breaks = "jenks")

sf はどのように地理座標を投影するのですか?

つまり、東経と北緯が x 軸と y 軸に直線的に写され、縦横比が 1 (東の 1 単位が北の 1 単位) になるように投影されます。地理データでは、座標が経度と緯度である場合、等緯経度投影 (equidistant circular とも呼ばれます)が選択され、プロットの中心 (または境界ボックス) では、北の1単位と東の1単位が等しいとします。

Proj.4では、この投影にデータを投影することもでき、その際のプロットは

plot(st_geometry(nc), axes = TRUE)

should, apart from the values along axes, be otherwise identical to

lat_ts = mean(st_bbox(nc)[c(2,4)]) # latitude of true scale
eqc = st_transform(nc, paste0("+proj=eqc +lat_ts=", lat_ts))
plot(st_geometry(eqc), axes = TRUE)

グラティキュール

等緯線とは、経度（子午線）や緯度（平行線）に沿いました格子線のことで、使用する投影法によっては、しばしば地図上に曲線で描かれ、経度や緯度の基準を与えることがあります。sf関数st_graticule` は、任意の地図に対してグラティキュールグリッドを作成しようとするものです。投影法は無限にあるので、うまくいかない場合も多いだろうが、そのような場合は sf issues として歓迎されます。

次のプロットは、それ自身に対するグラティキュール幾何のプロットです。

library(maps)
usa = st_as_sf(map('usa', plot = FALSE, fill = TRUE))
laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal area
usa <- st_transform(usa, laea)
g = st_graticule(usa)
plot(st_geometry(g), axes = TRUE)

ここで、このグラトリはプロットの境界には到達せず（しかし usa のバウンディングボックスで切り離されます）、軸は投影座標を表していることがわかります。

プロット関数の中で重心を計算するとき、プロット領域を知っているので、プロットの境界まで計算することができ、重心単位で軸を追加することができます。

plot(usa, graticule = TRUE, key.pos = NULL, axes = TRUE)

また、graticule に crs オブジェクトを渡すと、デフォルトのデータム (WGS84) とは異なるデータムのグラチュールを取得することができます。また、 st_graticule はパラメータを受け取るので、 plot の graticule パラメータにオブジェクトを渡して、より細かい制御を行うことができます。

g = st_graticule(usa, lon = seq(-130,-65,5))
plot(usa, graticule = g, key.pos = NULL, axes = TRUE,
     xlim = st_bbox(usa)[c(1,3)], ylim = st_bbox(usa)[c(2,4)],
     xaxs = "i", yaxs = "i")

ベースプロットのプロット領域を完全にコントロールすることは簡単ではありません。

他のパッケージで sf オブジェクトをプロットする

グリッド st_as_grob`

パッケージ sf は st_as_grob のためにいくつかのメソッドを提供しています。

methods(st_as_grob)
##  [1] st_as_grob.CIRCULARSTRING*      st_as_grob.COMPOUNDCURVE*      
##  [3] st_as_grob.CURVEPOLYGON*        st_as_grob.GEOMETRYCOLLECTION* 
##  [5] st_as_grob.LINESTRING*          st_as_grob.MULTILINESTRING*    
##  [7] st_as_grob.MULTIPOINT*          st_as_grob.MULTIPOLYGON*       
##  [9] st_as_grob.MULTISURFACE*        st_as_grob.POINT*              
## [11] st_as_grob.POLYGON*             st_as_grob.sfc*                
## [13] st_as_grob.sfc_CIRCULARSTRING*  st_as_grob.sfc_LINESTRING*     
## [15] st_as_grob.sfc_MULTILINESTRING* st_as_grob.sfc_MULTIPOINT*     
## [17] st_as_grob.sfc_MULTIPOLYGON*    st_as_grob.sfc_POINT*          
## [19] st_as_grob.sfc_POLYGON*        
## see '?methods' for accessing help and source code

これはシンプルフィーチャオブジェクトを grob (“graphics objects”) オブジェクトに変換するものです。grob は grid プロットパッケージのグラフィックプリミティブです。これらのメソッドは、 gridをベースにしているプロットパッケージ、例えば ggplot2 (これは geom_sf でこれらを使用しています) や tmap で使用することができます。さらに、 st_viewport を使用すると、 sf オブジェクトから base::plot.sf と同様の縦横比を持つグリッドビューポートをセットアップすることができます。

ggplot2

は、シンプルフィーチャオブジェクトのための特別な geom を含んでおり、 sf::st_graticule を使用してバックグラウンドでグラティキュール白線をサポートします。現在のところ、ポリゴンに対するサポートは良好です。線や点については、あなたのやり方次第です。

library(ggplot2)
ggplot() + geom_sf(data = usa)

多角形は aes を使って色をつけることができます。

ggplot() + 
  geom_sf(data = nc, aes(fill = BIR74)) + 
  scale_y_continuous(breaks = 34:36)

で、マップのセットは sf オブジェクトを再配置している後、ファセットプロットとしてプロットすることができます。

library(dplyr)
## 
##  次のパッケージを付け加えます: 'dplyr'
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     filter, lag
##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     intersect, setdiff, setequal, union
library(tidyr)
nc2 <- nc %>% select(SID74, SID79, geom) %>% gather(VAR, SID, -geom)
ggplot() + 
  geom_sf(data = nc2, aes(fill = SID)) + 
  facet_wrap(~VAR, ncol = 1) +
  scale_y_continuous(breaks = 34:36)

mapview

mapview パッケージは、 leaflet パッケージを利用して、html ページにインタラクティブな地図を作成します。豊富なサンプルが ここ にあります。

サンプルは以下のようにして得られます。

library(mapview)
mapviewOptions(fgb = FALSE) # needed when creating web pages
mapview(nc["BIR74"], col.regions = sf.colors(10), fgb = FALSE)

は、インタラクティブな地図を提供します。ズームやパンをしたり、クリックすることでフィーチャを問い合わせることができます。

tmap

パッケージ tmap は地図を描画するための別のパッケージであり、印刷できる地図に重点を置いています。

library(tmap)
qtm(nc)

また、tmapにはインタラクティブな leaflet の地図があります。

tmap_mode("view")
tm_shape(nc) + tm_fill("BIR74", palette = sf.colors(5))

非インタラクティブモードでの最後のマップの再描画は、次のように簡単です。

ttm()
tmap_last()

Martijn Tennekes と Jakub Nowosad による本 Elegant and informative maps with tmap のドラフト版は https://r-tmap.github.io/ で見つかります。