Skip to main content

疫学のための R ハンドブック

28 GIS の基本

28.1 概要

データに空間情報があることで、アウトブレイクの状況について多くの洞察が得られ、次のような質問に答えることができます。

  • 現在の疾病のホットスポットはどこか？
  • ホットスポットは時系列でどのように変化しているか？
  • 医療施設へのアクセスはどうなっているか？改善が必要か？

この GIS の章では、疫学業務担当者がアウトブレイク対応において知りたいニーズに応えることを目的としています。ここでは、tmap と ggplot2 パッケージを使って、基本的な空間データの可視化方法を学びます。また、sf パッケージを使用して、基本的な空間データの管理とクエリの方法についても説明します。最後に、spdep パッケージを用いて、空間的関係、空間的自己相関、空間回帰など、空間統計の概念に簡単に触れます。

28.2 重要な用語

ここでは、いくつかの重要な用語を紹介します。GIS と空間分析について詳しく知りたい方は、「リソース」節に掲載されている長いチュートリアルやコースをご覧になることをお勧めします。

GIS（Geographic Information System, 地理情報システム） - GIS とは、空間データを収集、管理、分析、視覚化するためのフレームワークまたは環境のことです。

GIS ソフトウェア

一般的な GIS ソフトウェアの中には、マウス操作で地図の作成や空間分析ができるものがあります。これらのツールには、コードを覚える必要がない、アイコンや機能を手動で選択して地図上に配置するのが簡単であるなどの利点があります。ここでは、人気のある2つのソフトウェアを紹介します。

ArcGIS - ESRI 社が開発した商用の GIS ソフトウェアで、非常に人気がありますますが、かなり高価です。

QGIS - フリーのオープンソース GIS ソフトウェアで、ArcGIS でできることもおおむねできます。ここでQGIS をダウンロードすることができます。

R を GIS として使用することは、「マウス操作」ではなく「コマンドライン・インターフェース」（目的の結果を得るためにはコーディングが必要）であるため、最初は敷居が高く感じられるかもしれません。しかし、繰り返し地図を作成したり、再現性のある分析を行う必要がある場合には、これは大きな利点となります。

空間データ

GIS で使用される空間データには、主にベクトルデータとラスターデータの 2 種類があります。

ベクトルデータ - GIS で使用される空間データの最も一般的な形式です。ベクトルデータは、頂点とパスの幾何学的なフィーチャ (feature; 「地物」とも訳される) で構成されています。ベクトルの空間データは、さらに3つの広く使用されるタイプに分類されます。

  • 点 - 点は、座標系における特定の位置を表す座標ペア(x,y)で構成されます。点は空間データの最も基本的な形態であり、地図上の症例（例：患者の家）や場所（例：病院）を表すのに使用されます。

  • 線 - 線は、2つの接続された点で構成されます。線には長さがあり、道路や川などを表すのに使われます。

  • ポリゴン - ポリゴンは、点で接続された少なくとも3つの線分で構成されています。ポリゴンのフィーチャは、面積だけでなく、長さ（すなわち、面の外周の長さ）を持っています。ポリゴンは、地域（例：村）や構造（例：病院の実際の面積）を示すときに使用されます。

ラスターデータ - 空間データの代替フォーマットであるラスターデータは、セル（ピクセルなど）のマトリックスで、各セルには高さ、温度、傾斜、森林被覆などの情報が含まれています。これは、航空写真や衛星画像などでよく見られます。ラスターデータは、ベクトルデータの下の「基図」としても使用されます。

空間データの可視化

GIS ソフトウェアでは、空間データを地図上に視覚的に表現するために、それぞれのフィーチャがどこにあるのか、お互いの関係について十分な情報を提供する必要があります。ベクトルデータを使用している場合、ほとんどのユース症例では、この情報は通常シェープファイルに保存されます。

シェープファイル - シェープファイルは、線、点、ポリゴンで構成される「ベクトル」空間データを保存するための一般的なデータフォーマットです。1つのシェープファイルは、実際には少なくとも3つのファイル（.shp、.shx、.dbf）の集合体です。シェープファイルを読み取るためには、これらのサブコンポーネントファイルがすべて同一のディレクトリ（フォルダ）に存在する必要があります。これらの関連ファイルは、ZIPフォルダに圧縮してメールで送信したり、ウェブサイトからダウンロードすることができます。

シェープファイルには、フィーチャ自体の情報だけでなく、地球の表面のどこにあるかという情報も含まれています。空間データをどのように「平面化」するかは、地図の見た目や解釈に大きな影響を与えるため、重要なポイントとなります。

座標参照系 (Coordinate Reference Systems, CRS) - CRS とは、地球上の地理的フィーチャを特定するために使用される座標ベースのシステムです。CRSにはいくつかの重要な要素があります。

  • 座標系 (Coordinate System) - 非常にたくさんの座標系があるため、自分が使用する座標がどの座標系を用いているのかを必ず確認してください。緯度/経度の度数を用いたものは一般的ですが、UTM 座標が利用されていることもあります。

  • 単位 (Units) - 座標系の単位を知ることができます（例：小数点で表された度、メートル）。

  • 測地系 (Datum) - 地球上の特定のモデルの、特定のバージョンを指します。測地系は長年にわたって改訂されてきたので、地図として表示されるレイヤー（マップレイヤー）が同じ測地系を使用していることを確認してください。

  • 投影 (Projection) - 実際には丸い地球を平らな表面（地図）に投影するために使用された数式への参照。

以下に示すマッピングツールを使用しなくても、空間データを要約することができることも覚えておいてください。地域別（例：地区、国など）のシンプルな表があれば十分なこともあります。

訳注: 日本で主に利用されている CRS は次のような ものがあります。

  • 平面直角座標系（世界測地系）I〜XIII EPSG: 2443-2455
  • JGD2000 GRS80楕円体 EPSG:4612
  • JGD2011 GRS80楕円体 EPSG:6668 （東日本大震災以降の大規模地殻変動を反映）

28.3 はじめての GIS

地図を作るためには、以下の通り、いくつかの重要な要素について、考え、また入手しなければありません。

  • データセット - これは空間データ形式（上述のようにシェープファイルなど）の場合もあれば、空間形式ではない場合もあります（例えば、単なる csv など）。

  • データセットが空間フォーマットでない場合は、参照データも必要です。参照データは、データの空間表現と関連する属性で構成されており、特定のフィーチャの位置情報やアドレス情報を含む資料などがあります。

    • 事前に定義された地理的境界線（例えば、行政区域）を扱う場合、参照シェープファイルは、政府機関やデータ共有組織から自由にダウンロードできることが多いです。迷ったときは、「[地域名] shapefile」とグーグルで検索してみるとよいでしょう。

    • 住所情報はあっても緯度・経度がないような場合は、 ジオコーディング・エンジン（geocoding engine）を使ってデータの空間参照データを取得する必要があるかもしれません。

  • データセットの情報を対象者に、どのように見せたいかというアイデア。地図には様々な種類がありますが、どの種類の地図が自分のニーズに最も適しているかを考えることが重要です。

訳注：日本の主なGISデータ： * 国勢調査 * 国土数値地図 * 国土地理院基盤地図 * G空間情報センター

データを可視化するための地図のタイプ

色分け（Choropleth）地図 - 主題図の一種で、色、濃淡、またはパターンを使って、ある属性の値に関連して地理的な地域を表現するもの。例えば、大きな値は小さな値よりも濃い色で示されます。このタイプの地図は、ある変数と、それが定義された地域や地政学的エリアでどのように変化するかを視覚化する場合に特に便利です。

症例密度ヒートマップ - 主題図の一種で、値の強さを色で表しますが、データをグループ化するために定義された地域や地政学的な境界線は使用しません。このタイプの地図は、一般的に「ホットスポット」や、点が高密度または集中しているエリアを示すために使用されます。

点密度マップ - は、データの属性値をドットで表現する主題図タイプです。このタイプのマップは、データの散らばりを視覚化し、クラスターを視覚的にスキャンするのに適しています。

等級シンボルマップ（段階のあるシンボルマップ） - 色分け地図に似ている主題図ですが、属性の値を色で表すのではなく、その値に対応するシンボル（通常は円）を使用します。例えば、大きな値は小さな値よりも大きなシンボルで示されます。このタイプの地図は、地域ごとのデータの大きさや量を視覚化したい場合に最適です。

また、複数の異なるタイプの視覚化を組み合わせて、複雑な地理的パターンを示すこともできます。例えば、下の地図の症例（点）は、最も近い医療施設に応じて色分けされています（凡例参照）。大きな赤い円は、一定の半径を持つ医療施設のキャッチメントエリアを示し、赤色の症例点は、キャッチメントの範囲外にある症例を示しています。

注：この GIS の章の主な焦点は、現場でのアウトブレイク対応の状況に基づいています。そのため、この章の内容は、基本的な空間データの操作、視覚化、および分析をカバーしています。

28.4 準備

パッケージを読み込む

以下のコードを実行すると、分析に必要なパッケージが読み込まれます。このハンドブックでは、パッケージを読み込むために、pacman パッケージの p_load() を主に使用しています。p_load() は、必要に応じてパッケージをインストールし、現在の R セッションで使用するためにパッケージを読み込む関数です。また、すでにインストールされたパッケージは、R の基本パッケージである base （以下、base R）の library() を使用して読み込むこともできます。R のパッケージに関する詳細はR の基礎の章をご覧ください。

pacman::p_load(
  rio,           # データをインポート
  here,          # ファイルの位置を探す
  tidyverse,     # データを処理、プロット (ggplot2 パッケージを含む)
  sf,            # Simple Featureフォーマットによる空間データの管理
  tmap,          # シンプルな地図を作成、インタラクティブな地図と静的な地図の両方に対応
  janitor,       # 列名のクリーニング
  OpenStreetMap, # ggplot map に OSM 基図を追加
  spdep          # 空間統計
  )

CRAN “Spatial Task View”では、空間データを扱うすべての R パッケージの概要を見ることができます。

症例データサンプル

ここでは分かりやすくするため、シミュレーションされたエボラ出血熱の linelist データフレームから1000件のランダムなサンプルを使って作業します（計算上、少ない件数で作業すると、このハンドブックでの表示が容易になります）。続いて、「前処理された」ラインリスト（linelist） をクリックしてダウンロードできます（.rdsファイル）。

ここではランダムにサンプルを取っているので、実際にコードを実行してみると、ここに表示されている結果とは若干異なる結果になるかもしれません。

データを rio パッケージの import() を使ってインポートします（rio パッケージは、.xlsx、.csv、.rds など様々な種類のファイルを扱えるパッケージです。詳細はデータのインポート・エクスポートの章を参照してください）。

# import clean case linelist
linelist <- import("linelist_cleaned.rds")

次に、base R の sample() を使って 1000 行のランダムなサンプルを抽出します。

# ラインリストの行数から1000個のランダムな行番号を生成
sample_rows <- sample(nrow(linelist), 1000)

# サンプル行の全ての列だけを保持するラインリストにサブセットする
linelist <- linelist[sample_rows,]

そして、データフレーム型である linelist を “sf”（空間フィーチャ）型のオブジェクトに変換していきます。この ラインリストには、各症例の居住地の経度と緯度を表す “lon” と “lat” の2つの列があるので、これは簡単です。

ここでは、sf（空間的フィーチャ）パッケージとその関数、 st_as_sf() を使用して、linelist_sf という新しいオブジェクトを作成します。この新しいオブジェクトは、基本的には linelist と同じですが、列 lon と lat を座標列として指定し、点を表示する際の座標参照系 (CRS) を割り当てています。4326では、GPS座標の標準である World Geodetic System 1984 (WGS84) に基づいた座標が指定されています。

# sfオブジェクトを生成
linelist_sf <- linelist %>%
     sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

これがオリジナルの linelist データフレームの中身です。このデモでは、date_onset と geometry（上記の経度と緯度のフィールドから作成されたもので、データフレームの最後の列）の列のみを使用します。

DT::datatable(head(linelist_sf, 10), rownames = FALSE, options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )

行政境界シェープファイル

シエラレオネ (Sierra Leone): 行政境界シェープファイル

事前に、Humanitarian Data Exchange (HDX) ウェブサイトからシエラレオネのすべての行政（admin）境界データをダウンロードしておきました。あるいは、ハンドブックとデータのダウンロードの章で説明しているように、私たちの R パッケージを介して、このデータやその他のサンプルデータをダウンロードすることもできます。

では、いよいよ Admin Level 3 のシェープファイルを R に保存するために、以下の作業を行います。

  1. シェープファイルをインポートする
  2. 列名を加工する
  3. 興味のある地域だけを残すために行を抽出する

シェープファイルをインポートするには、sf パッケージの read_sf() を使います。ファイルのパスは here() で与えられます。- 今回の場合、ファイルはRプロジェクトの “data”、“gis”、“shp” というサブフォルダに “sle_adm3.shp” として格納されています（詳細は データのインポート・エクスポート と R プロジェクトの設定 の章を参照してください）。独自のファイルパスを設定してください。

次に、janitor パッケージの clean_names() を使って、シェープファイルの列名を標準化します。また、filter() を使って、admin2name が “Western Area Urban” または “Western Area Rural” の行だけを残します。

# ADM3 level clean
sle_adm3 <- sle_adm3_raw %>%
  janitor::clean_names() %>% # 列名の標準化
  dplyr::filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # 特定の地域を残してフィルター

下の表は、インポートし、クリーニングした後のシェープファイルの様子です。右側にスクロールして、admin レベル 0（国）、admin レベル 1、admin レベル 2、そして最後に admin レベル 3 の列があることがわかります。それぞれのレベルには、キャラクター名と固有の識別子 “pcode” があります。たとえば、SL（Sierra Leone） → 04（Western） → L0410（Western Area Rural）→ L0040101（Koya Rural）のように、pcode は admin レベルが上がるごとに拡張されます。

人口データ

シエラレオネ: ADM3 ごとの人口

これらのデータは、HDX からダウンロードすることもできますし（リンクはこちら）、この章ハンドブックとデータのダウンロードで説明しているように、epirhandbook パッケージを介してダウンロードすることもできます。.csv ファイルの読み込みには、import() を使用します。また、インポートしたファイルを clean_names() に渡して列名の構文を標準化します。

# ADM3 ごとの人口
sle_adm3_pop <- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
  janitor::clean_names()

人口ファイルはこのようになっています。右にスクロールすると、各行政区域に male、female、total の列があり、さらに年齢層別の列に人口の内訳が表示されています。

医療施設

シエラレオネ: OpenStreetMap からの医療施設

今回も HDX から医療施設の位置情報をダウンロードしました。こちら または ハンドブックとデータのダウンロードの章の指示に従ってください。

施設の位置を示す点のシェープファイルを read_sf() でインポートし、列名をクリーニングしてから、病院（“hospital”）、診療所（“clinic”）、医師（“doctor”） のいずれかのタグが付けられた点だけを残すように抽出しました。

# OSM 医療施設の shapefile
sle_hf <- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>% 
  janitor::clean_names() %>%
  dplyr::filter(amenity %in% c("hospital", "clinic", "doctors"))

右にスクロールすると、施設名と geometry 座標が表示されます。

28.5 座標のプロット

X-Y座標（経度・緯度で表される点）をプロットする最も簡単な方法は、今回の場合、準備編で作成した linelist_sf オブジェクトから直接、点として描画することです。

tmap パッケージ は、わずか数行のコードで、静的（“plot” モード）とインタラクティブ（“view” モード）の両方に対応したシンプルなマッピング機能を提供します。tmap パッケージ の構文は ggplot2 パッケージ と似ていて、コマンド同士を + で追加していきます。詳しくはこの vignette をご覧ください。

  1. tmap モードを設定。ここでは plot モードを使い、静的なアウトプットを作成します。
tmap_mode("plot") # "view" または "plot" を選択

以下では、点を単独でプロットしています。tm_shape() は、linelist_sf オブジェクトとともに提供されます。次に、tm_dots() でサイズと色を指定して点を追加します。linelist_sf は sf オブジェクトなので、緯度・経度座標と座標系（CRS）を格納する2つの列はすでに指定されています。

# 症例（点）のみ
tm_shape(linelist_sf) + tm_dots(size=0.08, col='blue')

点だけでは、あまり意味がありません。そこで、行政境界もマッピングする必要があります。

ここでも tm_shape() を使いますが（documentation 参照）、症例の点のシェープファイルを提供する代わりに、行政境界のシェープファイル（ポリゴン）を提供します。

bbox = の引数（bbox は “bounding box” の略）で、座標の境界を指定することができます。まず bbox なしの地図表示を行い、次に bbox ありの地図表示を行います。

# 行政境界（ポリゴン）のみ
tm_shape(sle_adm3) +               # 行政境界のシェープファイル
  tm_polygons(col = "#F7F7F7")+    # ポリゴンを薄い灰色で表示
  tm_borders(col = "#000000",      # 境界を色と線の太さで表示
             lwd = 2) +
  tm_text("admin3name")            # 各ポリゴンについて表示する列テキスト


# 上と同様、ただしバウンディングボックス引数から縮尺を指定
tm_shape(sle_adm3,
         bbox = c(-13.3, 8.43,    # 角
                  -13.2, 8.5)) +  # 角
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")

そして、今度は点とポリゴンの両方を一緒に表示します。

# すべてまとめる
tm_shape(sle_adm3, bbox = c(-13.3, 8.43, -13.2, 8.5)) +     #
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")+
tm_shape(linelist_sf) +
  tm_dots(size=0.08, col='blue', alpha = 0.5) +
  tm_layout(title = "Distribution of Ebola cases")   # 地図にタイトルを付ける

R でのマッピングオプションの比較については、こちらの ブログ記事 をご覧ください。

28.6 空間結合

データセットと別のデータセットを 結合 することについては知っているでしょう。いくつかの方法は、このハンドブックのデータの結合の章で説明しています。空間結合は、同様の目的を持ちますが、空間的な関係を利用します。観測値を正しく一致させるために列の共通の値に頼るのではなく、あるフィーチャが他のフィーチャの中に入っているとか、他のフィーチャに最近傍であるとか、他のフィーチャから一定の半径のバッファ内にあるなど、空間的な関係を利用することができます。

sf パッケージには、空間結合のための様々なメソッドが用意されています。st_join メソッドや空間結合の種類については、こちらの参考資料に詳しい説明があります。

ポリゴン内の点

症例に行政単位を空間割当て

ここで興味深い問題があります。症例のラインリストには、症例の admin 単位に関する情報が含まれていません。このような情報は、最初のデータ収集段階で収集するのが理想的ですが、空間的な関係（つまり、点がポリゴンと交差する）に基づいて、個々の症例に admin 単位を割り当てることもできます。

以下では、症例の位置（点）と ADM3 の境界（ポリゴン）を空間的に交差させます。

  1. linelist（点）から始める
  2. 結合のタイプを “st_intersects” に設定して、境界線に空間的に結合する
  3. select() を使用して、新しい行政境界列のうち特定のものだけを残す
linelist_adm <- linelist_sf %>%
  
  # 空間交差に基づいて行政境界を linelist に結合
  sf::st_join(sle_adm3, join = st_intersects)

sle_adms のすべての列がラインリストに追加されました！各症例には、該当する行政区域レベルの詳細を示す列が追加されました。この例では、新しい列のうち 2 つ（admin レベル3）だけを残したいので、古い列名を select() して、追加したい2つだけを選択します。

linelist_adm <- linelist_sf %>%
  
  # 空間交差に基づいて、行政境界ファイルをラインリストに結合
  sf::st_join(sle_adm3, join = st_intersects) %>% 
  
  # 古い列名はそのままに、新たに2つの admin の列名を追加
  select(names(linelist_sf), admin3name, admin3pcod)

下の図は、最初の10件と、ポリゴン図形と空間的に交差する点に基づいて付けられた admin レベル3（ADM3）の行政境界を表示しています。

# これで、各症例に ADM3 名がついたことが確認できます。
linelist_adm %>% select(case_id, admin3name, admin3pcod)
## Simple feature collection with 1000 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.27122 ymin: 8.448701 xmax: -13.20595 ymax: 8.491748
## Geodetic CRS:  WGS 84
## First 10 features:
##      case_id     admin3name admin3pcod                   geometry
## 2282  05067a Mountain Rural   SL040102 POINT (-13.21619 8.451049)
## 168   34f7ac Mountain Rural   SL040102  POINT (-13.20947 8.45053)
## 448   7174dc       West III   SL040208 POINT (-13.25744 8.453545)
## 4776  24f9d8       West III   SL040208 POINT (-13.26871 8.460778)
## 1306  83f0cf Mountain Rural   SL040102 POINT (-13.21768 8.472346)
## 2619  18ee39       West III   SL040208  POINT (-13.26451 8.45718)
## 3591  410db8      Central I   SL040201  POINT (-13.23145 8.47372)
## 3336  483a29      Central I   SL040201  POINT (-13.23354 8.47932)
## 4617  80663c         West I   SL040206 POINT (-13.25057 8.487206)
## 16    c97dd9      Central I   SL040201 POINT (-13.22434 8.471451)

空間結合の前にはできなかった、行政単位での症例の表現が可能になりました。

# 行政単位ごとに症例の数を含むデータフレームを新規作成
case_adm3 <- linelist_adm %>%          # 新しい admin 列を含むラインリストで始める
  as_tibble() %>%                      # 見やすいように tibble に変換
  group_by(admin3pcod, admin3name) %>% # 名前と pcode で行政単位をグループ化
  summarise(cases = n()) %>%           # 要約と行数のカウント
  arrange(desc(cases))                 # 下り順に並べ替え

case_adm3
## # A tibble: 10 × 3
## # Groups:   admin3pcod [10]
##    admin3pcod admin3name     cases
##    <chr>      <chr>          <int>
##  1 SL040102   Mountain Rural   292
##  2 SL040208   West III         208
##  3 SL040207   West II          166
##  4 SL040204   East II          114
##  5 SL040201   Central I         66
##  6 SL040203   East I            63
##  7 SL040206   West I            44
##  8 SL040205   East III          24
##  9 SL040202   Central II        21
## 10 <NA>       <NA>               2

また、行政単位ごとの症例数の棒グラフを作成することもできます。

この例では、ggplot() を linelist_adm で始めているので、頻度で棒グラフを並べる fct_infreq() などの因子型のための関数を適用することができます（ヒントは 因子（ファクタ）型データ の章を参照してください）。

ggplot(
    data = linelist_adm,                     # admin unit 情報を含むラインリストで始める
    mapping = aes(
      x = fct_rev(fct_infreq(admin3name))))+ # x軸は行政単位、件数で並べ替え
  geom_bar()+                                # 帽を作成、高さは行数
  coord_flip()+                              # adm を読みやすくするためにXとYを入れ替える
  theme_classic()+   # 背景を単純に
  labs(                                      # タイトルとラベル
    x = "Admin level 3",
    y = "Number of cases",
    title = "Number of cases, by adminstative unit",
    caption = "As determined by a spatial join, from 1000 randomly sampled cases from linelist"
  )

最近傍

最寄の医療施設／集客エリアを探す

病気のホットスポットに関連して、医療施設がどこにあるかを知っておくと便利かもしれません。

ここでは、st_join() (sf パッケージ)の st_nearest_feature 結合メソッドを使って、個々の症例に最も近い医療施設を可視化することができます。

  1. まず、シェープファイルのラインリスト linelist_sf を用意します。
  2. 保健施設や診療所の位置（点）である sle_hf と空間的に結合します。
# 各症例に最も近い医療施設
linelist_sf_hf <- linelist_sf %>%                  # linelist shapefile で始める
  st_join(sle_hf, join = st_nearest_feature) %>%   # 症例データから最も近い診療所からのデータ
  select(case_id, osm_id, name, amenity) %>%       # 残しておくべき列、例えば id, name, type, と医療施設の位置情報
  rename("nearest_clinic" = "name")                # 分かりやすいように名前を変更

以下（最初の50行）のように、それぞれの症例には最寄りの診療所や病院のデータがあることがわかります。

約30％の症例で “Den Clinic” が最も近い医療機関であることがわかります。

# 医療施設ごとに症例を数える
hf_catchment <- linelist_sf_hf %>%   # 最寄りの診療所データを含むラインリストで始める
  as.data.frame() %>%                # shapefile をデータフレームに変換
  count(nearest_clinic,              # （診療所の）"name" で行を数える
        name = "case_n") %>%         # 数えたデータの列に "case_n" と命名
  arrange(desc(case_n))              # 下がり順で並べ替え

hf_catchment                         # console に表示
##                          nearest_clinic case_n
## 1       Shriners Hospitals for Children    352
## 2                            Den Clinic    336
## 3         GINER HALL COMMUNITY HOSPITAL    177
## 4                             panasonic     53
## 5 Princess Christian Maternity Hospital     36
## 6                     ARAB EGYPT CLINIC     21
## 7                                  <NA>     15
## 8                  MABELL HEALTH CENTER     10

結果を可視化するために、tmap パッケージ を使用することができます。

tmap_mode("view")   # tmap モードをインタラクティブに設定

# 症例と診療所の点をプロット
tm_shape(linelist_sf_hf) +            # 症例をプロット
  tm_dots(size=0.08,                  # 最も近い診療所で症例を色分け
          col='nearest_clinic') +    
tm_shape(sle_hf) +                    # 診療所を大きい黒い点でプロット
  tm_dots(size=0.3, col='black', alpha = 0.4) +      
  tm_text("name") +                   # 施設名でオーバーレイ
tm_view(set.view = c(-13.2284, 8.4699, 13), # 縮尺を調整 (中心座標, zoom)
        set.zoom.limits = c(13,14))+
tm_layout(title = "Cases, colored by nearest clinic")

バッファ

また、最も近い医療施設から徒歩2.5km（約30分）以内にある症例がどれくらいあるかを調べることもできます。

注: より正確な距離を計算するためには、UTM（地球を平面に投影したもの）などの各地域の地図投影系に sf オブジェクトを再投影するのがよいでしょう。この例では、簡単にするために、世界測地系（WGS84）の地理的座標系（地球は球状/円形の表面で表現され、そのため単位は10進法の度数）にしておきます。ここでは、一般的な換算として、1度＝約111kmとします。

地図投影法と座標系については、こちらのesriの記事をご覧ください。こちらのブログ では、様々なタイプの地図投影について、また、興味のある分野や地図・分析の文脈に応じて、どのように適切な投影を選ぶことができるかについて述べられています。

まず、各医療施設の周囲に半径約2.5kmの円形バッファを作成します。これは tmap パッケージ の st_buffer() で行います。地図の単位は緯度・経度の10進法であるため、“0.02” は度数として解釈されます。地図の座標系がメートル単位の場合は、数値もメートル単位で指定する必要があります。

sle_hf_2k <- sle_hf %>%
  st_buffer(dist=0.02)       # 約2.5kmに度数を変換

以下のように、バッファゾーン自体を作ります：

tmap_mode("plot")
# 円バッファを作成
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) +                    # 診療所を大きい赤丸でプロット
  tm_dots(size=0.3, col='black')

次に、st_join() の st_intersects 結合タイプを使って、これらのバッファと症例（点）を交差させます。これによって、バッファからのデータが、交差している点に結合されます。

# バッファで症例を交差
linelist_sf_hf_2k <- linelist_sf_hf %>%
  st_join(sle_hf_2k, join = st_intersects, left = TRUE) %>%
  filter(osm_id.x==osm_id.y | is.na(osm_id.y)) %>%
  select(case_id, osm_id.x, nearest_clinic, amenity.x, osm_id.y)

これによって結果を数えることができます。1000件のうち1000件は、どのバッファとも交差していない（値が欠けている）ので、最寄りの医療施設から徒歩30分以上の場所に住んでいることになります。

# どの医療施設バッファとも交差しなかった症例
linelist_sf_hf_2k %>% 
  filter(is.na(osm_id.y)) %>%
  nrow()
## [1] 1000

その結果、どのバッファとも交わらなかった症例が赤で表示されるように可視化されます。

tmap_mode("view")

# まず症例を点で表示
tm_shape(linelist_sf_hf) +
  tm_dots(size=0.08, col='nearest_clinic') +

# 診療所を大きい黒点でプロット
tm_shape(sle_hf) +                    
  tm_dots(size=0.3, col='black')+   

# 医療施設バッファをポリラインで重ねる
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2) +

# どの医療施設バッファにもない症例を赤点で強調
tm_shape(linelist_sf_hf_2k %>%  filter(is.na(osm_id.y))) +
  tm_dots(size=0.1, col='red') +
tm_view(set.view = c(-13.2284,8.4699, 13), set.zoom.limits = c(13,14))+

# タイトルを追加
tm_layout(title = "Cases by clinic catchment area")