ggplot2のsizeが意味するもの

付け焼き刃な知識故に、間違いございましたら指摘頂けますと幸いです。

ggplot2を使っていると、数種類の“size”を設定するタイミングがある。例えば、geom_point()で散布図を描画する際のポイントの大きさを指定するためのsize、または連続した線を描画するのに用いられるgeom_line()でのsize、あるいは文字を重ねるためのgeom_text()やannotate()の中でのsizeである。また、少し異なる引数名であるがテーマ指定の関数theme_*()にも“base_size”という名の引数が用意されている。ここで肝心なことは、テキストとその他では、同じsizeであっても単位が異なるということだ。このことは次のStackOverflowの投稿にまとめられている。

stackoverflow.com

まとめて書いてしまうとこの通り。

• レイヤーとして重ねるgeom_point()やgeom_line()などの要素のサイズはmm単位（ポイントの大きさ、線の太さなど）
• 軸やタイトルなど、フォントを表示するテキスト系のサイズはpoint単位(1pointsはおよそ0.352778mm)
• geom_point()のポイントの輪郭の太さを制御するstrokeの単位もpoint（ggplot2::.stroke）

以下、この記事ではそのことを確かめてみたいと思う。

適当なggオブジェクトを用意し、それに手を加えて比較していくことにする。

library(ggplot2)
# 日本語の表示のためにフォントを指定
theme_set(theme_gray(base_family = "IPAexGothic"))

p <- 
  ggplot(mtcars, aes(wt, mpg))

まずはgeom_point()やgeom_line()でのsizeである。ここではsize = 10を指定した(既定では1.5)。

p_out <- 
  p + 
  geom_point(size = 10)

この結果をggsave()で出力した図を表示する。ポイントのサイズが10mmにしたので、検証のために出力時のサイズはその10倍となる10cmとした。

ggsave(filename = "out.png", 
       plot = p_out, 
       width = 10, 
       height = 10, 
       units = "cm")

確かにmmらしい。もっとも、これはソースコードにも示されている。

Size examples Should be specified with a numerical value (in millimetres), or from a variable source

github.com

次に、軸やタイトルテキストのsizeである。軸のラベルを除いて、size = 10を指定した。

p2 <- 
  p_out +
  # 適当な場所に文字「A」を表示
  annotate("text", x = 2, y = 15, size = 10, label = "A", color = "red") +
  ggtitle("文字の大きさはpointで指定")

p2 + 
  theme(title = element_text(size = 10),
         axis.title = element_text(size = 10),
         # グラフの両軸に用いられる数値の大きさ
         axis.text = element_text(size = 5))

フォントを扱っていて、なおかつ同じsize = 10を指定していてもannotate()が表示するフォントが大きい。これはannotate()がレイヤーの要素として評価されるためで、geom_point()のようにsizeはmm単位での指定となるためだ。

もし、タイトルや軸のサイズをmm単位にしたい時はpointをmmへ変換する必要がある。ggplot2::.ptはggplot2の内部で用いられている変数で、ドキュメントもある

p2 +
  theme(axis.title = element_text(size = ggplot2::.pt * 10))

最後にtheme_*()の引数“base_size”だが、こちらも単位はpointである。すなわち次の2つのコードによって出力される図は等しい。

ggplot2が生成するような描画オブジェクトの差分を確認するにはvdiffrが便利である。ここではwidget_diff()を用いて差分の検証を行うが、GitHubのように2つの画像を表示し、スライドさせながら差分を確認するwidget_slide()やテスト時に役立つ関数も備わっている。

qiita.com

# plotとaxisのラベルの大きさは、テーマのbase_sizeに対してI()で指定した値をかけた大きさとなる
p_default <- 
  p2 + 
  theme_gray(base_size = 11,
             base_family = "IPAexGothic")

p_custom <- 
  p2 +
  theme_replace(
    plot.title = element_text(size = 11 * I(1.2)),
    axis.title = element_text(size = 11 * I(0.8)))

# 2つのオブジェクトで差分がない場合は白い画面が表示される
vdiffr::widget_diff(
  p_default,
  p_custom
)

# 変更してみる
# フォントサイズが変更されたことで全体の配置も微妙にずれる
vdiffr::widget_diff(
  p_default,
  p_custom +
    theme(plot.title = element_text(size = 20))
)

学術誌に投稿する図では、軸のテキストサイズやらについて細かく指定しされているものもある。怒られないように、少しずつggplot2に慣れていこうと思う。 Enjoy!

地物範囲のタイル座標を得る

id:yutannihilationさんが、ggplot2::geom_sf()にOpen Street Map (OSM)のタイルを重ねるという面白い試みをしています。

yutani.rbind.io

OSMでは、画像データとして地図画像を配信していて、次のようなURLで参照されます。

http://[abc].tile.openstreetmap.org/zoom/x/y.png

ドメイン名の他、zoom、xのディレクトリがあり、ファイルはyという名前で置かれています。ここでzoomは地図のズームレベル、xとyはタイルの座標を示す、ピクセル座標と呼ばれる値です。OSMをはじめ、地図タイルの配信を行なっているサービスは同様にこの「XYZ地図タイル」形式を採用しています。

OSMに限らず、目的の地物を範囲に収めるタイルを用意するには、このピクセル座標を取得する必要があるわけですが、馴染みのある緯度経度からピクセル座標の変換は少し頭を使わなくてはいけません（私だけかもしれませんが）。今回の話題は、すでにyutannihilationさんの記事にも書かれていますが、緯度経度からピクセル座標を取得する処理をおさらいをするというものです（二番煎じです。本当にありがとうございました）。

基礎

緯度経度およびズームレベルからピクセル座標を求める処理は、OSMのWikiに書かれている通りです。計算式の詳細は参考にあげるページに詳しいのでそちらをご覧ください。

library(magrittr)

Wikiに書かれているコードをちょっと修正して、緯度経度とズームレベルを与えた時にピクセル座標とズームレベルをリストで返す関数を定義します。

deg2num <- function(longitude, latitude, zoom){
  lat_rad <- latitude * pi /180
  n <- 2.0 ^ zoom
  xtile <- floor((longitude + 180.0) / 360.0 * n)
  ytile = floor((1.0 - log(tan(lat_rad) + (1 / cos(lat_rad))) / pi) / 2.0 * n)
  return(list(xtile = xtile, 
              ytile = ytile,
              zoom = zoom))
}

実行結果です。

deg2num(140.113889, 35.613056, 10)

## $xtile
## [1] 910
## 
## $ytile
## [1] 403
## 
## $zoom
## [1] 10

deg2num(140.113889, 35.613056, 17)

## $xtile
## [1] 116549
## 
## $ytile
## [1] 51643
## 
## $zoom
## [1] 17

タイル座標が判明したので、次は一枚のタイルがカバーする領域の座標を求めることにします。これには、タイルの一角の座標を取得したあとで、周辺に隣接するタイルの座標を求めるという方法を採用します。まずは先ほどのタイル座標の領域に含まれる緯度経度座標の北東となる座標を取得する処理です。

num2deg <- function(xtile, ytile, zoom) {
  n = 2.0 ^ zoom
  longitude = xtile / n * 360.0 - 180.0
  lat_rad = atan(sinh(pi * (1 - 2 * ytile / n)))
  latitude = lat_rad / 0.01745329
  return(list(longitude = longitude, 
              latitude = latitude,
              zoom = zoom))
}

deg2num(140.113889, 35.613056, 10) %>% 
  purrr::pmap(
    ~ num2deg(..1, ..2, ..3)
  )

## [[1]]
## [[1]]$longitude
## [1] 139.921875
## 
## [[1]]$latitude
## [1] 35.7465174211
## 
## [[1]]$zoom
## [1] 10

同様の処理を隣接する全てのタイル座標で実行することで、元のタイル座標の領域を取得できます。これはポリゴンとして扱えるようにしておきます。

tile_bbox <- function(xtile, ytile, zoom) {
  
  ytile <- c(ytile, ytile + 1, ytile + 1, ytile, ytile)
  xtile <- c(xtile, xtile, xtile + 1, xtile + 1, xtile)
  
  purrr::map2_dfr(
  xtile,
  ytile,
  ~ num2deg(.x, .y, zoom)) %>% 
    dplyr::select(-zoom) %>% 
    as.matrix() %>% 
    list() %>% 
    sf::st_polygon() %>% 
    sf::st_sfc(crs = 4326)
}

deg2num(140.113889, 35.613056, 10) %>% 
  purrr::pmap(
    ~ tile_bbox(..1, ..2, ..3)
  )

## [[1]]
## Geometry set for 1 feature 
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 139.921875 ymin: 35.4606750714 xmax: 140.2734375 ymax: 35.7465174211
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

確認のめに一度描画してみましょう。北緯 35.613056、東経 140.113889をズームレベル10のタイル座標で示すと次のようになります。

library(mapview)

mapview(deg2num(140.113889, 35.613056, 10) %>% 
  purrr::pmap(
    ~ tile_bbox(..1, ..2, 10)
  ) %>% .[[1]])

本当かどうか怪しいので、さらに国土地理院のタイル座標確認ページと比較します。作成したポリゴンは、次のURLを表示した時に描画される領域と一致しています。

http://maps.gsi.go.jp/development/tileCoordCheck.html#10/35.6131/140.1139

作成したポリゴンは、次のURLを表示した時に描画される領域と一致しているようで一安心です。

地物を覆う領域のタイル座標

さてここからが本番です。一点の座標ではなく、地物を覆う領域のタイル座標にはどうすれば良いでしょう。答えは単に、対象の地物の矩形からタイルを作成することです。例として、「岡山県」を覆うタイルを描画してみましょう。

library(dplyr)
library(sf)
library(jpndistrict)

df_33 <- 
  jpn_pref(33)

b <- 
  df_33 %>% 
  # 「岡山県」の矩形を作成し、その領域を求めます
  st_bbox()

b

##          xmin          ymin          xmax          ymax 
## 133.266694605  34.298390466 134.413218589  35.352861842

st_bbox()により、地物の境界となる座標がわかったので、領域を覆うタイル座標を生成します。これにはexpand.grid()を使っても良いですが、ここではtidyr::crossing()を用いました。また、元のズームレベルも変数として残すようにしています。

df_tiles <- 
  tidyr::crossing(
    x = seq(
      deg2num(b["xmin"], b["ymin"], z = 10)[["xtile"]],
      deg2num(b["xmax"], b["ymax"], z = 10)[["xtile"]]
    ),
    y = seq(
      deg2num(b["xmin"], b["ymin"], z = 10)[["ytile"]],
      deg2num(b["xmax"], b["ymax"], z = 10)[["ytile"]]
    )
  ) %>% 
  dplyr::mutate(zoom = 10)

dplyr::glimpse(df_tiles)

## Observations: 20
## Variables: 3
## $ x    <int> 891, 891, 891, 891, 891, 892, 892, 892, 892, 892, 893, 89...
## $ y    <int> 404, 405, 406, 407, 408, 404, 405, 406, 407, 408, 404, 40...
## $ zoom <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1...

「岡山県」を覆うタイル座標を得ることができました。このタイル座標を再びポリゴンに変換します。

df_tiles <- 
  df_tiles %>% 
  dplyr::mutate(geometry = purrr::pmap(., ~ tile_bbox(..1, ..2, ..3) %>% 
                                         st_sf(crs = 4326))) %>% 
  tidyr::unnest() %>% 
  st_sf(crs = 4326)

おまけですが、各タイルの重心も変数として保存しておきましょう。

df_tiles <- 
  df_tiles %>% 
  dplyr::mutate(longitude = purrr::pmap_dbl(., ~ st_centroid(..4)[1]),
                latitude = purrr::pmap_dbl(., ~ st_centroid(..4)[2]))

library(ggplot2)

ggplot() +
  geom_sf(data = df_33) +
  geom_sf(data = df_tiles, alpha = 0.5) + 
  geom_text(data = df_tiles, 
            aes(longitude, latitude, label = paste(zoom, x, y, sep = "/")), 
            color = "red")

やりました！岡山県をカバーする範囲のタイル座標を重ねることができました。

ですが、右下のタイルなど、一部のタイルには「岡山県」が含まれていません。最後に、地物が含まれないタイルを除外する処理を加えておきましょう。これまでの処理を関数としてまとめ、一度に実行できるようにしておきます。

st_tiles <- function(x, zoom) {
  
  b <- st_bbox(x)
  zoom <- rlang::enquo(zoom)
 
  tidyr::crossing(
    x = seq(
      deg2num(b["xmin"], b["ymin"], z = !!zoom)[["xtile"]],
      deg2num(b["xmax"], b["ymax"], z = !!zoom)[["xtile"]]
    ),
    y = seq(
      deg2num(b["xmin"], b["ymin"], z = !!zoom)[["ytile"]],
      deg2num(b["xmax"], b["ymax"], z = !!zoom)[["ytile"]]
    )
  ) %>% 
    dplyr::mutate(zoom = !! zoom) %>% 
    dplyr::mutate(geometry = purrr::pmap(., ~ tile_bbox(..1, ..2, ..3) %>% 
                                           st_sf(crs = 4326))) %>% 
    tidyr::unnest() %>% 
    st_sf(crs = 4326) %>% 
    dplyr::mutate(longitude = purrr::pmap_dbl(., ~ st_centroid(..4)[1]),
                  latitude = purrr::pmap_dbl(., ~ st_centroid(..4)[2])) %>% 
    dplyr::mutate(check = purrr::pmap(., ~ st_intersects(..6,
                                                         !!rlang::eval_tidy(x),
                                                         sparse = FALSE) %>%
                                        rowSums()) > 0) %>%
    dplyr::filter(check == TRUE)
}

ズームレベルも11に変更します。

df_tiles_intersect <- 
  st_tiles(df_33, 11)

ggplot() +
  geom_sf(data = df_33) +
  geom_sf(data = df_tiles_intersect, alpha = 0.5) + 
  geom_text(data = df_tiles_intersect, 
            aes(longitude, latitude, label = paste(zoom, x, y, sep = "/")),
            color = "red",
            size = 1.8)

ちょっと長いコードになってしまったので、もっと簡単な書き方があれば教えてください。

Enjoy!

sfオブジェクトを描画する方法まとめ: ggplot2::geom_sf()もあるよ

Rの作図パッケージとして人気なggplot2の時期バージョンの2.3.0が間も無くリリースされるそうです。ggplot2は前回の更新が2016年12月末のバージョン2.2.1なので、久しぶりのバージョンアップとなります。

バージョンアップに伴う変更点はこのページを見て欲しいですが、私としてはsfオブジェクトを描画するための新たな関数geom_sf()に触れずには要られません。これまでずっと開発版でのみ利用可能でしたが、やっとCRAN版でも利用できるようになります。（人に紹介するときはGitHubの説明からやったり面倒でもあった。）

... と思っていたらid: yutannihilationさんに先を越されてしまいました。

ggplot2 2.3.0（RC版）を使ってみた - Technically, technophobic.

notchained.hatenablog.com

同じ話題をしても仕方がないので、ここでは、ggplot2に実装されるgeom_sf()の紹介を兼ねて、sfパッケージにより作成されるsfオブジェクトを描画する方法を整理しようと思います。また、sfオブジェクトの可視化に祭して、私の知っているいくつかのtipsを紹介します。(私の持ち合わせる知識の範囲です。他にあればぜひ教えてください。)

公式のvignettesにもsfオブジェクトを描画する方法について整理されています。

https://cran.r-project.org/web/packages/sf/vignettes/sf5.html

sfオブジェクトを描画する際はまず、静的あるいは動的に描画するという目的によって使い分けます。静的な描画はRの標準グラフィックスに利用されるplot()、冒頭のggplot2などです。一方、動的な描画は、興味のある箇所を拡大したり、表示・非表示をユーザが切り替え可能な描画の方法です。具体的にはleaflet、plotlyがこの方法での選択肢となります。

私は、静的な出力の場合はgeom_sf()、動的にズームしたり動かしたい場合はmapview()というような使いわけをしています（少し前はmapviewではなくてleafletでしたが鞍替えしました）。

用意

作図用のデータとして、国土数値情報の提供する行政区域のデータと市区町村役場データを用います。このデータはjpndistrictによりRから取得可能です。

library(sf)
library(jpndistrict)

## Loading required package: jpmesh
## This package provide map data is based on the Digital Map 25000(Map Image) published by
## Geospatial Information Authorityof Japan (Approval No.603FY2017 information usage
## <http://www.gsi.go.jp>)

# 岡山県の行政区域データ
# この記事の中の岡山県のポリゴンデータを使った図は、国土地理院長の承認を得て、
# 同院発行の数値地図（国土基本情報）電子国土基本図（地図情報）を使用したものです (承認番号　平28情使、第603号)
pref_33 <- 
  jpn_pref(33)
head(pref_33)

## Simple feature collection with 6 features and 4 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 133.6025 ymin: 34.4173 xmax: 134.169 ymax: 35.30721
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 5
##   pref_code prefecture city_code city                             geometry
##   <chr>     <chr>      <chr>     <chr>                      <GEOMETRY [°]>
## 1 33        岡山県     33101     岡山市 北区 POLYGON ((133.9098 34.948, 1…
## 2 33        岡山県     33102     岡山市 中区 MULTIPOLYGON (((133.9747 34.…
## 3 33        岡山県     33103     岡山市 東区 MULTIPOLYGON (((133.9958 34.…
## 4 33        岡山県     33104     岡山市 南区 MULTIPOLYGON (((133.9109 34.…
## 5 33        岡山県     33202     倉敷市      MULTIPOLYGON (((133.6414 34.…
## 6 33        岡山県     33203     津山市      POLYGON ((134.0871 35.30721,…

# 岡山県内の市役所等のポイントデータ
df_office33 <- 
  jpn_admins(33)

## options:        ENCODING=cp932 
## Reading layer `P34-14_33' from data source `/private/var/folders/0x/mb63hycs4k30_7httqyxh2rr0000gn/T/RtmpVmQLFw/P34-14_33_GML/P34-14_33.shp' using driver `ESRI Shapefile'
## Simple feature collection with 105 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 133.3315 ymin: 34.35497 xmax: 134.3596 ymax: 35.28274
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs

head(df_office33)

## Simple feature collection with 6 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 133.8795 ymin: 34.65511 xmax: 134.039 ymax: 34.90654
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs
##   jis_code type           name                    address
## 1    33100    1     岡山市役所        岡山市北区大供1-1-1
## 2    33101    2       御津支所     岡山市北区御津金川1020
## 3    33101    2       建部支所    岡山市北区建部町福渡489
## 4    33103    2       瀬戸支所     岡山市東区瀬戸町瀬戸45
## 5    33101    2     鶴田連絡所 岡山市北区建部町角石谷2063
## 6    33101    1 岡山市北区役所        岡山市北区大供1-1-1
##                    geometry
## 1 POINT (133.9196 34.65511)
## 2  POINT (133.9354 34.7992)
## 3 POINT (133.9038 34.86932)
## 4  POINT (134.039 34.73388)
## 5 POINT (133.8795 34.90654)
## 6 POINT (133.9196 34.65511)

plot

sfパッケージには、作図関数のplot()の総称関数が用意されています。そのため、sfが生成する、sfg (simple feature geometry)、sfc (simple feature geometry list-column)、sf (simple feature)のいずれかのクラスに属するsfオブジェクトであれば、plot()を適用することで、難なく地物を描画することが可能です。

par(mfrow = c(1, 2))
# sfg
st_point(0:1) %>% 
  plot()
# sfc
st_sfc(st_point(0:1), 
       st_point(2:1)) %>% 
  plot(col = "red")

sfcを含んだsfオブジェクト（データフレームの変数 + geometryのようなもの)を描画すると全ての変数を対象に描画が行われます。

# 既定では全ての変数について描画を行う
plot(pref_33)

これは一旦の確認には良いですが、最終的な出力では一変数だけを出力することが多いかと思います。その際はst_geometry()をかまして地物情報だけを抽出し、sfcとして描画するか、特定列を参照した状態で描画を行うと良いでしょう。

# 地物情報だけを描画
st_geometry(pref_33) %>% 
  plot()

# 特定列だけの描画
# 凡例の表示のため、余白を調節します
par(oma = c(3, 3, 3, 6.5), mar = c(2, 1, 1, 2))
pref_33["city"] %>% 
  plot(family = "IPAPGothic")

ggplot2

次に今回の目玉となるggplot2でのgeom_sf()を使う方法です。基本的なことは冒頭にあげたyutannihilationさんの記事に書かれています。

library(ggplot2)

ggplot(pref_33) +
  geom_sf()

ggplot2に詳しい方ならご存知かと存じますが、ggplot2 では描画したいデータをggplot()の第一引数dataに指定します（ここでは岡山県の行政区域データ）。それを+で連結させ、どのように描画するかという際にgeom_*()を利用します。今回実装されるgeom_sf()はsfオブジェクトを描画する関数で、これまで散布図ならgeom_point()、棒グラフならgeom_bar()のように使い分けてきたように、sfオブジェクトであればgeom_sf()という塩梅となります。

ggplot(data = )のgeom_sf(data = )に渡す方法もあります。これは、描画したいオブジェクトが複数（例えばポリゴンとポイント）ある際に、それぞれのgeom_sf()で対象のデータを指定する時に、どのデータをどのにレイヤーに渡すかを区別させるために利用されます。

ggplot() +
  geom_sf(data = pref_33, aes(fill = city_code)) +
  geom_sf(data = df_office33, color = "white")

次のコードでも同じ図が得られますが、geom_sf()にどのデータが渡っているかがわかりにくいので、上記の方法が好みです。

ggplot(pref_33) +
  geom_sf(aes(fill = city_code)) +
  geom_sf(data = df_office33, color = "white")

より地図っぽく描画するには nozmaさんの記事にあるようにcoords_sf()やtheme_void()を使うと良いでしょう。

ggplot() +
  geom_sf(data = pref_33, aes(fill = city_code)) +
  # 緯度経度の線、ラベルを描画しない
  coord_sf(datum = NA) +
  # 背景色を白にする
  theme_void()

tmap

tmapは主題図を描くために適したパッケージですが、こちらもsfオブジェクトに対応しています。描画のための関数はqtm()およびtm_shape()です。

library(tmap)

qtm(pref_33)

tm_shape()ではmodeと呼ばれる機能があり、それを切り替えることで静的・動的な描画の切り替えが可能です。

# 結果は上の図と同じなので省略
tmap_mode("plot")

tm_shape(pref_33, projection = "longlat") + 
    tm_polygons()

インタラクティブに動かせる図を描画するにはモードをviewに切り替えます。動作の様子をgif画像で表示します。

tmap_mode("view")

tm_shape(pref_33) + 
  tm_fill("city_code", palette = sf.colors(20))

leaflet

leafletはJavaScript製の地図描画ライブラリですが、Rに実装されたleafletパッケージはsfオブジェクトを描画し、ユーザが任意の箇所を拡大したり表示箇所を移動することが可能です。tmapやmapviewの基礎として利用されているライブラリでもあります。

leafletでsfを表示するには、対象とするsfのgeometryの種類を把握した上で適した関数に渡す必要があります。すなわちポイントの場合は、addCircles()、ラインはaddPolylines()、ポリゴンはaddPolygons()といった具合です。

library(leaflet)

leaflet() %>% 
  addTiles() %>% 
  # ポイントの表示なのでaddCircles()
  # 日本のへそ公園の位置を表示
  addCircles(data = st_point(c(135, 35.0)) %>% 
               st_sfc() %>% 
               st_sf(crs = 4326))

leaflet() %>% 
  addTiles() %>% 
  # ポリゴンの表示
  addPolygons(data = pref_33)

mapview

ggplot2のようにオブジェクトを重ねていくことができますが、種類に応じて関数を分ける必要がある点や、ベースタイルの表示や表示の切り替えに別の関数を利用しなければならないという点でちょっと面倒です。そこで、もっとサクッと表示したいという時にmapviewが便利です。

ggplot2のように複数のオブジェクトを+演算子で重ねていくことができ、地物の種類も関係なくmapview()で済むので楽チンです（spオブジェクトにも対応します）。

library(mapview)

mapview(pref_33) 

地図描画のためのオプションも豊富で一押しです。

mapview(pref_33, 
        zcol = "city_code",
        legend = FALSE) +
  mapview(st_point(c(133.9196, 34.65511)),
          color = "red", 
          col.regions = "red")

plotly

最後に、leafletライブラリとは別のアプローチでsfを動的に可視化する方法としてplotlyの例です。こちらはまだ開発版の機能となりますが、plotlyの開発ブログ記事に書書かれているように、plotlyの特徴を活かした3Dでの表示や他の種類のグラフと組み合わせが可能なようです。

こちらは開発版の機能ですので、将来実装が変更される可能性があります。

remotes::install_github("ropensci/plotly")
library(plotly)

nc <- 
  sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
plot_geo(nc)

以上です。Enjoy!