空间数据可视化地图绘制R语言可复现
空間數據可視化&地圖繪制&R語言可復現
繪制地理空間數據是一項常見的可視化任務，需要專門的工具。通常，問題可以分解為兩個問題：
有許多可以用于繪制地圖的工具，一些流行的 GIS 軟件允許點擊交互進行地圖開發和空間分析。這些工具具有不需要學習代碼以及在地圖上手動選擇和放置圖標和功能的便利性等優點，例如：
ArcGIS - 由 ESRI 公司開發的商業 GIS 軟件，非常流行但相當昂貴
QGIS - 一個免費的開源 GIS 軟件，幾乎可以做 ArcGIS 可以做的任何事情
使用 R 作為 GIS 一開始似乎更令人生畏，因為它不是“點擊”，而是有一個“命令行界面”（您必須編寫代碼才能獲得所需的結果）。但是，如果您需要重復生成地圖或創建可重現的分析，這是一個主要優勢。
本文主要解釋使用R語言繪制地圖，主要用到的包有:
- rio : 導入數據
- tidyverse : 清洗、處理、繪圖（包含ggplot2包）
- sf : 使用簡單要素格式管理空間數據
- tmap : 生成地圖，適用于交互式和靜態地圖
- OpenStreetMap : 在ggplot中免費添加 OSM 國內外地圖作為底圖
本目文錄：
- 空間數據類型 Spatial data
- 矢量數據 Vector data
- 柵格數據 Raster data
- 空間數據基本格式
- 基本地圖類型
- 必備環節
- 所需R包
- 實例案例數據
- 行政邊界數據
- 人口數據
- 衛生設施數據上
- 點數據與多邊形區域連接
- 最近的鄰域連接
- 緩沖區連接
- 其他空間連接類型
- 免費獲取國內外地圖OpenStreetMap
- 等高線密度熱圖 Contoured density heatmap
- 時間序列熱圖 Time series heatmap
1.綜述
您的數據的空間方面可以提供很多關于爆發情況的見解，并回答以下問題：
- 當前的疾病熱點在哪里？
- 隨著時間的推移，熱點發生了怎樣的變化？
- 衛生設施的使用情況如何？是否需要任何改進？
此 GIS 頁面的當前重點是解決應用流行病學家在疫情應對中的需求。我們將探索使用 tmap 和 ggplot2 包的基本空間數據可視化方法。我們還將通過 sf 包介紹一些基本的空間數據管理和查詢方法。最后，我們將使用 spdep 包簡要介紹空間統計的概念，例如空間關系、空間自相關和空間回歸。
1.1空間數據類型 Spatial data
**地理信息系統 (GIS) **- GIS 是用于收集、管理、分析和可視化空間數據的框架或環境，GIS 中使用的空間數據的兩種主要形式是矢量Vector和柵格Raster數據：
矢量數據 Vector data
GIS 中使用的最常見的空間數據格式，矢量數據由頂點vertices和路徑paths的幾何特征組成。矢量空間數據可以進一步分為三種廣泛使用的類型：
- 點 Points - 點由表示坐標系中特定位置的坐標對 (x,y) 組成。點是最基本的空間數據形式，可用于在地圖上表示一個病例（即患者家）或位置（即醫院）。
- 線 Lines - 一條線由兩個相連的點組成。線有長度，可用于表示道路或河流等事物。
- 多邊形 Polygons - 多邊形由至少三個由點連接的線段組成。多邊形特征具有長度（即區域的周長）以及面積測量值。多邊形可用于標注區域（即村莊）或結構（即醫院的實際區域）。
柵格數據 Raster data
一種替代格式 空間數據，柵格數據是一個單元格矩陣（例如像素），每個單元格包含高度、溫度、坡度、森林覆蓋等信息。這些通常是航空照片、衛星圖像等。柵格也可以用作“底圖base maps”放在矢量數據下方。
1.2空間數據基本格式
為了在地圖上直觀地表示空間數據，GIS 軟件要求您提供足夠的信息，說明不同要素應在何處相互關聯。如果您使用的是矢量數據，這對于大多數用例都是正確的，則此信息通常會存儲在 shapefile 中：
Shapefiles - shapefile 是一種常見的數據格式，用于存儲由線、點或多邊形組成的“矢量”空間數據.單個 shapefile 實際上是至少三個文件的集合 - .shp、.shx 和 .dbf。所有這些子組件文件必須存在于給定目錄（文件夾）中，以便 shapefile 可讀。這些相關文件可以壓縮到 ZIP 文件夾中，以便通過電子郵件發送或從網站下載。 shapefile 將包含有關要素本身的信息，以及它們在地球表面的位置。這很重要，因為雖然地球是一個球體，但地圖通常是二維的。關于如何的選擇 “展平flatten”空間數據會對最終地圖的外觀和解釋產生重大影響。
坐標參考系統 (CRS) - CRS 是一種基于坐標的系統，用于定位地球表面的地理特征。它有幾個關鍵組件：
-
坐標系 Coordinate System - 有許多不同的坐標系，因此請確保您知道您的坐標來自哪個系統。緯度/經度的度數很常見，但您也可以看到 UTM 坐標。
-
單位 Units - 了解坐標系的單位（例如十進制度 decimal degrees、米）
-
基準 Datum - 地球的特定模型版本。多年來，這些已被修訂，因此請確保您的地圖圖層使用相同的基準。
-
投影 Projection - 對用于將真正的圓形地球投影到平坦表面（地圖）的數學方程式的參考。
請記住，您可以在不使用下面顯示的映射工具的情況下匯總空間數據。有時只需要一張按地理（例如地區、國家等）的簡單表格！
1.3用于空間數據可視化的基本地圖類型
分區統計地圖 Choropleth map - 一種專題地圖，其中顏色、陰影或圖案用于表示與其屬性值相關的地理區域。例如，較大的值可以用比較小值更深的顏色來表示。這種類型的地圖在可視化變量以及它如何在定義的區域或地緣政治區域之間變化時特別有用。
案例密度熱圖 Case density heatmap - 一種主題圖，其中顏色用于表示值的強度，但是，它不使用定義的區域或地緣政治邊界對數據進行分組。這種類型的地圖通常用于顯示“熱點”或具有高密度或點集中的區域。
點密度圖 Dot density map - 一種專題地圖類型，使用點來表示數據中的屬性值。這種類型的地圖最適合用于可視化數據的分散情況并直觀地掃描集群。
比例符號地圖（分級符號地圖）Proportional symbols map (graduated symbols map) - 類似于等值線地圖的專題地圖，但不是使用顏色來指示屬性的值，而是使用與值相關的符號（通常是圓形）。例如，較大的值可以用比較小值更大的符號來表示。當您想要跨地理區域可視化數據的大小或數量時，最好使用這種類型的地圖。
您還可以組合幾種不同類型的可視化來顯示復雜的地理模式。例如，下圖中的病例（點）根據其最近的醫療機構（參見圖例）進行著色。大紅色圓圈表示一定半徑的醫療機構集水區，而鮮紅色的案例點表示在任何范圍之外的區域：
2.繪制地圖準備工作
2.1制作地圖時的幾個關鍵項目，其中包括：
- 數據集——可以是空間數據格式（例如 shapefile，如上所述），也可以不是空間格式（例如，只是 csv）。
- 如果您的數據集不是空間格式，您還需要一個參考數據集 reference dataset。參考數據由數據的空間表示和相關屬性組成，其中包括包含特定要素的位置和地址信息的材料。
- 如果您使用預定義的地理邊界（例如，行政區域），通常可以從政府機構或數據共享組織免費下載參考 shapefile。如有疑問，最好從谷歌“[regions] shapefile”開始。
- 如果您有地址信息，但沒有緯度和經度，您可能需要使用地理編碼引擎來獲取空間參考數據以供記錄。
關于如何將數據集中的信息呈現。有許多不同類型的地圖，重要的是要考慮哪種類型的地圖最適合您的需求。
2.2導入所需要的包
此代碼塊顯示分析所需的包的加載。在本文，我們強調來自 pacman 的 p_load()，它會在必要時安裝包并加載它以供使用。您還可以使用基本 R 中的 library() 加載已安裝的包。
if(!require(pacman))install.packages("pacman")pacman::p_load(rio, # to import datahere, # to locate filestidyverse, # to clean, handle, and plot the data (includes ggplot2 package)sf, # to manage spatial data using a Simple Feature formattmap, # to produce simple maps, works for both interactive and static mapsjanitor, # to clean column namesOpenStreetMap, # to add OSM basemap in ggplot mapspdep, # spatial statisticsshowtext # 支持ggplot顯示中文)2.3示例案例數據
出于演示目的，我們將使用來自模擬的埃博拉流行病線列表數據框中的 1000 個病例的隨機樣本（在計算上，使用較少的病例更容易在本手冊中顯示）。linelist_cleaned.rds,數據集獲取方式見文末。
由于我們對案例進行隨機抽樣，因此當您自行運行代碼時，您的結果可能與此處演示的結果略有不同。
使用 rio 包中的 import() 函數導入數據（它處理許多文件類型，如 .xlsx、.csv、.rds ）。
# 使用import函數導入數據 linelist = import(here("data", "Epidemiologist_R","linelist_cleaned.rds"))接著使用sample()函數隨機抽取1000個樣本
# 隨機抽取1000個行標簽 sample_rows = sample(nrow(linelist), 1000) # 子集只包含樣本行和所有的列 linelist = linelist[sample_rows,]現在我們想將linelist(dataframe格式)轉換為“sf”類（空間特征）的對象，其中linelist 有兩列“lon”和“lat”代表每個案例居住地的經度和緯度。
我們使用包 sf（空間特征）及其函數 st_as_sf() 來創建我們稱為 linelist_sf 的新對象。這個新對象看起來與 linelist 基本相同，但 lon 和 lat 列已指定為坐標列，并且已為顯示點時分配了坐標參考系 (CRS)。 4326 表示我們的坐標標識為基于 1984 年世界大地測量系統 (WGS84) - 這是 GPS 坐標的標準。
# 創建一個空間數據集 linelist_sf <- linelist %>%sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)這就是原始 linelist 數據框的樣子。在此演示中，我們將僅使用列 date_onset 和geometry（由上面的經度和緯度字段構成，是數據框中的最后一列）。
head(linelist_sf) ## Simple feature collection with 6 features and 28 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -13.26978 ymin: 8.453176 xmax: -13.21946 ymax: 8.484195 ## Geodetic CRS: WGS 84 ## case_id generation date_infection date_onset date_hospitalisation date_outcome outcome gender age age_unit age_years age_cat age_cat5 hospital infector source wt_kg ht_cm ct_blood fever chills cough aches vomit temp time_admission bmi days_onset_hosp geometry ## 5680 78703f 23 2014-12-10 2014-12-13 2014-12-16 2014-12-20 Death m 16 years 16 15-19 15-19 Missing 5bf94a other 57 167 17 yes no no no yes 39.3 15:01 20.43817 3 POINT (-13.25821 8.453866) ## 3020 525a03 14 2014-10-21 2014-11-05 2014-11-07 2014-11-16 Death f 5 years 5 5-9 5-9 St. Mark's Maternity Hospital (SMMH) 5bc9be other 46 71 23 yes no yes no yes 38.9 11:25 91.25174 2 POINT (-13.22085 8.453176) ## 5154 19d49c 15 2014-09-06 2014-09-09 2014-09-16 <NA> Death f 10 years 10 10-14 10-14 Port Hospital bf6471 other 39 101 17 yes no yes no no 38.5 06:39 38.23155 7 POINT (-13.26978 8.480407) ## 3963 76a224 21 2015-01-05 2015-01-20 2015-01-21 2015-01-27 Recover <NA> 15 years 15 15-19 15-19 Military Hospital 674201 other 59 152 22 yes no yes no no 40.0 09:16 25.53670 1 POINT (-13.22583 8.484195) ## 2268 5355c3 16 <NA> 2014-10-04 2014-10-05 2014-10-17 <NA> f 24 years 24 20-29 20-24 Port Hospital <NA> <NA> 63 153 23 yes no yes yes no 39.7 12:52 26.91273 1 POINT (-13.21946 8.481999) ## 3539 185ac2 20 2014-11-19 2014-12-10 2014-12-11 2014-12-20 Death m 5 years 5 5-9 5-9 Military Hospital 4d1c91 other 48 93 22 yes no yes no no 38.8 13:25 55.49775 1 POINT (-13.25619 8.483108) # 運行以下代碼，得到交互式表格 # DT::datatable(head(linelist_sf, 10), rownames = FALSE, # options = list(pageLength = 8, scrollX=T), class = 'white-space: nowrap' )2.4塞拉利昂行政邊界 Admin boundary shapefile
預先，我們已從或者，您可以通過我們的 The epidemiologist R下載所有示例數據或,者獲取方式見文末。
現在我們將執行以下操作以在 R 中保存 Admin Level 3 shapefile：
要導入 shapefile，我們使用 sf 中的 read_sf() 函數。它通過 here() 提供文件路徑。 - 在我們的例子中，該文件位于我們的 R 項目中的“data”、“gis”和“shp”子文件夾中，文件名為“sle_adm3.shp”（有關更多信息，請參見導入和導出和 R 項目頁面）。您需要提供自己的文件路徑。接下來我們使用 janitor 包中的 clean_names() 來標準化 shapefile 的列名。我們還使用 filter() 只保留 admin2name 為“Western Area Urban”或“Western Area Rural”的行。
# ADM3 level clean sle_adm3_raw = st_read(here("data","Epidemiologist_R", "sle_adm3.shp")) ## Reading layer `sle_adm3' from data source `/Users/cpf/Documents/paper/writting_blog/data/Epidemiologist_R/sle_adm3.shp' using driver `ESRI Shapefile' ## Simple feature collection with 167 features and 19 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -13.30901 ymin: 6.923379 xmax: -10.27056 ymax: 9.999253 ## Geodetic CRS: WGS 84 sle_adm3 <- sle_adm3_raw %>%clean_names() %>% # standardize column namesfilter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # filter to keep certain areas您可以在下面看到導入和清理后 shapefile 的外觀。
head(sle_adm3) ## Simple feature collection with 6 features and 19 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -13.29416 ymin: 8.094272 xmax: -12.91333 ymax: 8.494073 ## Geodetic CRS: WGS 84 ## objectid admin3name admin3pcod admin3ref_n admin2name admin2pcod admin1name admin1pcod admin0name admin0pcod date valid_on valid_to shape_leng shape_area rowcacode0 rowcacode1 rowcacode2 rowcacode3 geometry ## 1 155 Koya Rural SL040101 Koya Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.63822264 0.0136601326 SLE SLE004 SLE004001 SLE004001001 MULTIPOLYGON (((-13.02082 8... ## 2 156 Mountain Rural SL040102 Mountain Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.29268365 0.0031823435 SLE SLE004 SLE004001 SLE004001002 MULTIPOLYGON (((-13.21496 8... ## 3 157 Waterloo Rural SL040103 Waterloo Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.72265524 0.0136413976 SLE SLE004 SLE004001 SLE004001003 MULTIPOLYGON (((-13.12215 8... ## 4 158 York Rural SL040104 York Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 1.23916989 0.0198344753 SLE SLE004 SLE004001 SLE004001004 MULTIPOLYGON (((-13.24441 8... ## 5 159 Central I SL040201 Central I Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.06882179 0.0001882636 SLE SLE004 SLE004002 SLE004002001 MULTIPOLYGON (((-13.22646 8... ## 6 160 East I SL040203 East I Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.05749169 0.0001429506 SLE SLE004 SLE004002 SLE004002003 MULTIPOLYGON (((-13.2129 8.... # # DT::datatable(head(sle_adm3, 10), rownames = FALSE, # options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )2.5人口數據 Popluation Ddata
塞拉利昂：ADM3 的人口
這些數據可以再次通過 Epirhandbook R 包下載，使用 import() 來加載 .csv 文件。我們還將導入的文件傳遞給 clean_names() 以標準化列名語法。
# Population by ADM3 sle_adm3_pop <- import(here("data", "Epidemiologist_R", "sle_admpop_adm3_2020.csv")) %>%clean_names()這是人口文件的樣子,可以查看每個轄區的男性人口、女性人口、總人口以及按年齡組分列的人口列。
head(sle_adm3_pop) ## adm0_en adm0_pcode adm1_en adm1_pcode adm2_en adm2_pcode adm3_en adm3_pcode female male total t_00_04 t_05_09 t_10_14 t_15_19 t_20_24 t_25_29 t_30_34 t_35_39 t_40_44 t_45_49 t_50_54 t_55_59 t_60_64 t_65_69 t_70_74 t_75_79 t_80plus ## 1 Sierra Leone SL Southern SL03 Bo SL0301 Badjia SL030101 4440 4630 9070 1201 1418 1084 1117 847 777 555 539 382 310 240 142 144 95 84 51 86 ## 2 Sierra Leone SL Southern SL03 Bo SL0301 Bagbo SL030102 14275 14585 28859 3820 4513 3450 3556 2695 2471 1765 1713 1219 988 762 450 458 300 266 162 274 ## 3 Sierra Leone SL Southern SL03 Bo SL0301 Bagbwe(Bagbe) SL030103 11645 11687 23331 3088 3648 2788 2874 2179 1998 1427 1385 984 798 615 363 370 243 215 130 222 ## 4 Sierra Leone SL Southern SL03 Bo SL0301 Bo Town SL030191 93653 100759 194412 25732 30401 23239 23949 18154 16648 11891 11541 8208 6652 5127 3032 3087 2021 1796 1089 1847 ## 5 Sierra Leone SL Southern SL03 Bo SL0301 Boama SL030104 25066 26037 51104 6763 7991 6109 6295 4772 4376 3125 3034 2157 1748 1348 797 812 531 472 287 486 ## 6 Sierra Leone SL Southern SL03 Bo SL0301 Bumpe Ngao SL030105 24214 25154 49369 6535 7720 5901 6082 4610 4228 3019 2930 2084 1689 1302 770 784 513 456 277 469 # # DT::datatable(head(sle_adm3_pop, 10), rownames = FALSE, # options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )2.6衛生設施 health Facilities
塞拉利昂：來自 OpenStreetMap 的衛生設施數據
我們使用 read_sf() 導入設施點 shapefile，再次清理列名，然后過濾以僅保留標記為“醫院hospital”、“診所clinic”或“醫生doctors”的點。
sle_hf <- sf::read_sf(here("data", "Epidemiologist_R", "sle_hf.shp")) %>% clean_names() %>%filter(amenity %in% c("hospital", "clinic", "doctors")) head(sle_hf) ## Simple feature collection with 6 features and 11 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -13.264 ymin: 8.38842 xmax: -13.14276 ymax: 8.490478 ## Geodetic CRS: WGS 84 ## # A tibble: 6 × 12 ## osm_id source addrfull building healthcare operatorty addrcity name amenity healthca_1 capacitype geometry ## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <POINT [°]> ## 1 3197529881 UNMEER ; https://data.hdx.rwlabs.org/dataset/ebola-treatment-centers <NA> <NA> <NA> <NA> <NA> China-SL Friendship Hospital, Jui hospital <NA> <NA> (-13.14276 8.38842) ## 2 3197529888 UNMEER ; https://data.hdx.rwlabs.org/dataset/ebola-treatment-centers <NA> <NA> <NA> <NA> <NA> Lakka Hospital ETU hospital <NA> <NA> (-13.264 8.397) ## 3 3212073781 https://data.hdx.rwlabs.org/dataset/sierra-leone-ebola-care-facilities <NA> <NA> <NA> <NA> <NA> Princess Christian Maternity Hospital hospital <NA> <NA> (-13.21935 8.489861) ## 4 3339906977 MSF-CH <NA> <NA> <NA> <NA> <NA> COMMUNITY CLINIC clinic <NA> <NA> (-13.16954 8.441156) ## 5 3341831443 MSF-CH <NA> <NA> <NA> <NA> <NA> Den Clinic clinic <NA> <NA> (-13.24653 8.483513) ## 6 3341855623 <NA> <NA> <NA> <NA> <NA> <NA> MABELL HEALTH CENTER clinic <NA> <NA> (-13.22632 8.490479) # DT::datatable(head(sle_hf, 10), rownames = FALSE, # options = list(pageLength = 5, scrollX=T), class = 'white-space: nowrap' )3.畫坐標系 Plotting coordinates
這種情況下，繪制 X-Y 坐標（經度/緯度、點）的最簡單方法是直接從我們在準備部分創建的 linelist_sf 對象中將它們繪制為點。
tmap 包為靜態static（“繪圖”模式）和交互式interactive（“視圖”模式）提供了簡單的映射功能，只需幾行代碼。 tmap 的語法與 ggplot2 的語法類似，詳情點擊。
下面，這些點是單獨繪制的。tm_shape() 與 linelist_sf 對象一起提供。然后我們通過 tm_dots() 添加點，指定大小和顏色。因為 linelist_sf 是一個 sf 對象，我們已經指定了包含緯度/經度坐標和坐標參考系（CRS）的兩列：
# 繪制點 tm_shape(linelist_sf) + tm_dots(size = 0.08, col = 'blue')單獨來看，這些點并不能告訴我們太多。所以我們還應該繪制行政邊界：
我們再次使用 tm_shape()，但我們不提供案例點 shapefile，而是提供管理邊界 shapefile（多邊形）。
使用 bbox = argument（bbox 代表“邊界框”），我們可以指定坐標邊界。首先我們顯示沒有bbox的地圖顯示，然后再使用它。
# Just the administrative boundaries (polygons) # # 只是行政邊界（多邊形） tm_shape(sle_adm3) + # admin boundaries shapefiletm_polygons(col = "#F7F7F7")+ # show polygons in light greytm_borders(col = "#000000", # show borders with color and line weightlwd = 2) +tm_text("admin3name") # column text to display for each polygon # Same as above, but with zoom from bounding box # # 同上，但從邊界框縮放 tm_shape(sle_adm3,bbox = c(-13.3, 8.43, # corner 左上角-13.2, 8.5)) + # cornertm_polygons(col = "#F7F7F7") +tm_borders(col = "#000000", lwd = 2) +tm_text("admin3name")現在把點和多邊形邊界放在一起
# All together 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") # give title to map3.空間連接 Spatial joins
您可能熟悉將數據從一個數據集連接到另一個數據集。空間連接具有類似的目的，但利用了空間關系。您可以利用它們的空間關系，而不是依賴列中的公共值來正確匹配觀測值，例如一個要素在另一個要素內，或與另一個要素最近的鄰居，或在距另一個要素一定半徑的緩沖區內，等等。
sf 包提供了多種空間連接方法。請參閱此參考中有關 st_join 方法和空間連接類型的更多文檔。
3.1多邊形中的點 Points in polygon
空間分配行政單位到案例 Spatial assign administrative units to cases
這是一個有趣的難題：案例行列表不包含有關案例行政單位的任何信息。盡管在初始數據收集階段收集此類信息是理想的，但我們也可以根據其空間關系（即點與多邊形相交）將行政單元分配給各個案例。
下面，我們將在空間上將案例位置（points）與 ADM3 邊界（多邊形polygons）相交：
sle_adms 中的所有列都已添加到行列表中！現在，每個案例都有列詳細說明其所屬的管理級別。在此示例中，我們只想保留兩個新列（admin level 3），因此我們select(）舊列名和另外兩個感興趣的列：
linelist_adm <- linelist_sf %>%# join the administrative boundary file to the linelist, based on spatial intersectionsf::st_join(sle_adm3, join = st_intersects) %>% # Keep the old column names and two new admin ones of interestselect(names(linelist_sf), admin3name, admin3pcod)下面，僅出于顯示目的，您可以看到前十個案例以及已附加的管理級別 3 (ADM3) 管轄區，基于點在空間上與多邊形形狀相交的位置。
# Now you will see the ADM3 names attached to each case 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.27095 ymin: 8.447961 xmax: -13.20589 ymax: 8.490442 ## Geodetic CRS: WGS 84 ## First 10 features: ## case_id admin3name admin3pcod geometry ## 5680 78703f West III SL040208 POINT (-13.25821 8.453866) ## 3020 525a03 Mountain Rural SL040102 POINT (-13.22085 8.453176) ## 5154 19d49c West III SL040208 POINT (-13.26978 8.480407) ## 3963 76a224 East II SL040204 POINT (-13.22583 8.484195) ## 2268 5355c3 East II SL040204 POINT (-13.21946 8.481999) ## 3539 185ac2 West II SL040207 POINT (-13.25619 8.483108) ## 5418 798126 West III SL040208 POINT (-13.26513 8.47565) ## 3874 e3ba2c West II SL040207 POINT (-13.23237 8.460837) ## 2481 b345c5 West II SL040207 POINT (-13.22861 8.469415) ## 2482 3ca785 West III SL040208 POINT (-13.25258 8.457118)現在我們可以按行政單位來描述我們的案例——這是在空間連接之前我們無法做到的！
# Make new dataframe containing counts of cases by administrative unit # 創建包含行政單位案件計數的新數據框case_adm3 <- linelist_adm %>% # begin with linelist with new admin colsas_tibble() %>% # convert to tibble for better displaygroup_by(admin3pcod, admin3name) %>% # group by admin unit, both by name and pcode summarise(cases = n()) %>% # summarize and count rowsarrange(desc(cases)) # arrange in descending orderhead(case_adm3) ## # A tibble: 6 × 3 ## # Groups: admin3pcod [6] ## admin3pcod admin3name cases ## <chr> <chr> <int> ## 1 SL040102 Mountain Rural 281 ## 2 SL040208 West III 214 ## 3 SL040207 West II 204 ## 4 SL040204 East II 108 ## 5 SL040201 Central I 54 ## 6 SL040203 East I 54我們還可以按行政單位創建病例計數條形圖 bar plot。
在此示例中，我們以 linelist_adm 開始 ggplot()，以便我們可以應用 fct_infreq() 之類的因子函數，它按頻率對條形圖進行排序（有關提示，請參閱因子頁面）。
ggplot(data = linelist_adm, # begin with linelist containing admin unit infomapping = aes(x = fct_rev(fct_infreq(admin3name))))+ # x-axis is admin units, ordered by frequency (reversed)geom_bar()+ # create bars, height is number of rowscoord_flip()+ # flip X and Y axes for easier reading of adm unitstheme_classic()+ # simplify backgroundlabs( # titles and labelsx = "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")3.2最近的鄰居 Nearest neighboor
尋找最近的醫療機構/集水區
了解與疾病熱點相關的衛生設施的位置可能很有用。
我們可以使用 st_join() 函數（sf 包）中的 st_nearest_feature 連接方法來可視化離個別病例最近的醫療機構。
我們可以看到，約 30% 的病例中，“Den Clinic”是最近的醫療機構。
# Count cases by health facility hf_catchment <- linelist_sf_hf %>% # begin with linelist including nearest clinic dataas.data.frame() %>% # convert from shapefile to dataframecount(nearest_clinic, # count rows by "name" (of clinic)name = "case_n") %>% # assign new counts column as "case_n"arrange(desc(case_n)) # arrange in descending orderhf_catchment # print to console ## nearest_clinic case_n ## 1 Den Clinic 369 ## 2 Shriners Hospitals for Children 321 ## 3 GINER HALL COMMUNITY HOSPITAL 176 ## 4 panasonic 43 ## 5 Princess Christian Maternity Hospital 33 ## 6 <NA> 23 ## 7 MABELL HEALTH CENTER 22 ## 8 ARAB EGYPT CLINIC 13 # tmap_mode("view") # set tmap mode to interactive # plot the cases and clinic points tm_shape(linelist_sf_hf) + # plot casestm_dots(size=0.08, # cases colored by nearest cliniccol='nearest_clinic') + tm_shape(sle_hf) + # plot clinic facilities in large black dotstm_dots(size=0.3, col='black', alpha = 0.4) + tm_text("name") + # overlay with name of facility tm_view(set.view = c(-13.2284, 8.4699, 13), # adjust zoom (center coords, zoom)set.zoom.limits = c(13,14))+ tm_layout(title = "Cases, colored by nearest clinic")3.3緩沖區 Buffers
我們還可以探索有多少病例位于距離最近的醫療機構 2.5 公里（約 30 分鐘）的步行距離內。
注意：為了更準確地計算距離，最好將您的 sf 對象重新投影到相應的本地地圖投影系統，例如 UTM（地球投影到平面上）。在此示例中，為簡單起見，我們將堅持使用世界大地測量系統 (WGS84) 地理坐標系（地球以球形/圓形表面表示，因此單位為十進制度）。我們將使用一般轉換：1 十進制度 = ~111km。
在這篇 esri 文章中查看有關地圖投影和坐標系的更多信息。該博客討論了不同類型的地圖投影，以及如何根據感興趣的區域和地圖/分析的背景選擇合適的投影。
現在我們可以計算結果：nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y),0])在 1000 個案例中沒有與任何緩沖區相交（缺少該值），因此步行超過 30 分鐘最近的醫療機構。
# Cases which did not get intersected with any of the health facility buffers linelist_sf_hf_2k %>% filter(is.na(osm_id.y)) %>%nrow() ## [1] 1000我們可以將結果可視化，使得不與任何緩沖區相交的案例顯示為紅色。
# tmap_mode("view")# First display the cases in points 首先以點顯示案例 tm_shape(linelist_sf_hf) +tm_dots(size=0.08, col='nearest_clinic') +# plot clinic facilities in large black dots 用大黑點繪制診所設施 tm_shape(sle_hf) + tm_dots(size=0.3, col='black')+ # Then overlay the health facility buffers in polylines 然后在折線中覆蓋衛生設施緩沖區 tm_shape(sle_hf_2k) +tm_borders(col = "black", lwd = 2) +# Highlight cases that are not part of any health facility buffers 用紅點突出顯示不屬于任何衛生設施緩沖區的病例 # in red dots 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))+# add title tm_layout(title = "Cases by clinic catchment area")3.4其他空間連接類型
參數join的替代值包括文檔:
- st_contains_properly
- st_contains
- st_covered_by
- st_covers
- st_crosses
- st_disjoint
- st_equals_exact
- st_equals
- st_is_within_distance
- st_nearest_feature
- st_overlaps
- st_touches
- st_within
4.等值域圖 Choropleth 地圖
Choropleth 地圖可用于按預定義區域（通常是行政單位或衛生區域）可視化您的數據。例如，在疫情應對中，這有助于針對高發病率的特定地區分配資源。
現在我們已經為所有案例分配了行政單位名稱（請參閱上面的空間連接部分），我們可以開始按區域映射案例計數（等值域圖）。由于我們還有 ADM3 的人口數據，我們可以將此信息添加到之前創建的 case_adm3 表中。
我們從上一步中創建的數據框 case_adm3 開始，它是每個行政單位及其案例數量的匯總表。
將此表與 ADM3 多邊形 shapefile 連接以進行映射
case_adm3_sf <- case_adm3 %>% # begin with cases & rate by admin unitleft_join(sle_adm3, by="admin3pcod") %>% # join to shapefile data by common columnselect(objectid, admin3pcod, # keep only certain columns of interestadmin3name = admin3name.x, # clean name of one columnadmin2name, admin1name,cases, total, case_10kpop,geometry) %>% # keep geometry so polygons can be plotteddrop_na(objectid) %>% # drop any empty rowsst_as_sf() # convert to shapefile # tmap mode tmap_mode("plot") # view static map# plot polygons tm_shape(case_adm3_sf) + tm_polygons("cases") + # color by number of cases columntm_text("admin3name") # name display我們還可以繪制發病率圖
# Cases per 10K population tmap_mode("plot") # static viewing mode# plot tm_shape(case_adm3_sf) + # plot polygonstm_polygons("case_10kpop", # color by column containing case ratebreaks=c(0, 10, 50, 100), # define break points for colorspalette = "Purples" # use a purple color palette) +tm_text("admin3name") # display text5.使用 ggplot2
進行映射如果您已經熟悉使用 ggplot2，則可以使用該包來創建數據的靜態映射。 geom_sf() 函數將根據數據中的特征（點、線或多邊形）繪制不同的對象。例如，您可以在 ggplot() 中使用 geom_sf()，使用 sf 數據和多邊形幾何來創建等值線圖。
為了說明這是如何工作的，我們可以從我們之前使用的 ADM3 多邊形 shapefile 開始。回想一下，這些是塞拉利昂的 3 級管理員區域：
sle_adm3 ## Simple feature collection with 12 features and 19 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -13.29894 ymin: 8.094272 xmax: -12.91333 ymax: 8.499809 ## Geodetic CRS: WGS 84 ## First 10 features: ## objectid admin3name admin3pcod admin3ref_n admin2name admin2pcod admin1name admin1pcod admin0name admin0pcod date valid_on valid_to shape_leng shape_area rowcacode0 rowcacode1 rowcacode2 rowcacode3 geometry ## 1 155 Koya Rural SL040101 Koya Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.63822264 1.366013e-02 SLE SLE004 SLE004001 SLE004001001 MULTIPOLYGON (((-13.02082 8... ## 2 156 Mountain Rural SL040102 Mountain Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.29268365 3.182343e-03 SLE SLE004 SLE004001 SLE004001002 MULTIPOLYGON (((-13.21496 8... ## 3 157 Waterloo Rural SL040103 Waterloo Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.72265524 1.364140e-02 SLE SLE004 SLE004001 SLE004001003 MULTIPOLYGON (((-13.12215 8... ## 4 158 York Rural SL040104 York Rural Western Area Rural SL0401 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 1.23916989 1.983448e-02 SLE SLE004 SLE004001 SLE004001004 MULTIPOLYGON (((-13.24441 8... ## 5 159 Central I SL040201 Central I Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.06882179 1.882636e-04 SLE SLE004 SLE004002 SLE004002001 MULTIPOLYGON (((-13.22646 8... ## 6 160 East I SL040203 East I Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.05749169 1.429506e-04 SLE SLE004 SLE004002 SLE004002003 MULTIPOLYGON (((-13.2129 8.... ## 7 161 East II SL040204 East II Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.08397463 1.494827e-04 SLE SLE004 SLE004002 SLE004002004 MULTIPOLYGON (((-13.22653 8... ## 8 162 Central II SL040202 Central II Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.04878431 6.513035e-05 SLE SLE004 SLE004002 SLE004002002 MULTIPOLYGON (((-13.23154 8... ## 9 163 West III SL040208 West III Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.30219417 1.698296e-03 SLE SLE004 SLE004002 SLE004002008 MULTIPOLYGON (((-13.28529 8... ## 10 164 West I SL040206 West I Western Area Urban SL0402 Western SL04 Sierra Leone SL 2016-08-01 2016-10-17 <NA> 0.06952030 1.822183e-04 SLE SLE004 SLE004002 SLE004002006 MULTIPOLYGON (((-13.24677 8...我們可以使用 dplyr 的 left_join() 函數來添加我們想要映射到 shapefile 對象的數據。在這種情況下，我們將使用之前創建的 case_adm3 數據框來匯總行政區域的病例數；但是，我們可以使用相同的方法來映射存儲在數據框中的任何數據。
sle_adm3_dat <- sle_adm3 %>% inner_join(case_adm3, by = "admin3pcod") # inner join = retain only if in both data objectsselect(sle_adm3_dat, admin3name.x, cases) # print selected variables to console ## Simple feature collection with 9 features and 2 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -13.29894 ymin: 8.384533 xmax: -13.12612 ymax: 8.499809 ## Geodetic CRS: WGS 84 ## admin3name.x cases geometry ## 1 Mountain Rural 281 MULTIPOLYGON (((-13.21496 8... ## 2 Central I 54 MULTIPOLYGON (((-13.22646 8... ## 3 East I 54 MULTIPOLYGON (((-13.2129 8.... ## 4 East II 108 MULTIPOLYGON (((-13.22653 8... ## 5 Central II 25 MULTIPOLYGON (((-13.23154 8... ## 6 West III 214 MULTIPOLYGON (((-13.28529 8... ## 7 West I 42 MULTIPOLYGON (((-13.24677 8... ## 8 West II 204 MULTIPOLYGON (((-13.25698 8... ## 9 East III 11 MULTIPOLYGON (((-13.20465 8...要按區域制作病例計數柱形圖，使用 ggplot2，我們可以調用 geom_col()，如下所示：
ggplot(data=sle_adm3_dat) +geom_col(aes(x=fct_reorder(admin3name.x, cases, .desc=T), # reorder x axis by descending 'cases' #通過降序“案例”對 x 軸重新排序y=cases)) + # y axis is number of cases by regiontheme_bw() +labs( # set figure texttitle="Number of cases, by administrative unit",x="Admin level 3",y="Number of cases") + guides(x=guide_axis(angle=45)) # angle x-axis labels 45 degrees to fit better如果我們想使用 ggplot2 來制作案例計數的等值線圖，我們可以使用類似的語法來調用 geom_sf() 函數：
ggplot(data=sle_adm3_dat) + geom_sf(aes(fill=cases)) # set fill to vary by case count variable然后，我們可以使用與 ggplot2 一致的語法來自定義地圖的外觀，例如：
ggplot(data=sle_adm3_dat) + geom_sf(aes(fill=cases)) + scale_fill_continuous(high="#54278f", low="#f2f0f7") + # change color gradienttheme_bw() +labs(title = "Number of cases, by administrative unit", # set figure textsubtitle = "Admin level 3")對于習慣使用 ggplot2 的 R 用戶，geom_sf() 提供了一個簡單直接的實現，適用于基本的地圖可視化。要了解更多信息，請閱讀 geom_sf() 小插圖或 ggplot2 書。
6.底圖 Basemaps
6.1 OpenStreetMap
下面我們將描述如何使用 OpenStreetMap 功能為 ggplot2 地圖實現底圖。替代方法包括使用 ggmap，它需要在 Google 上免費注冊（詳情）。
OpenStreetMap 是一個協作項目，旨在創建免費的可編輯世界地圖。基礎地理位置數據（例如城市位置、道路、自然特征、機場、學校、醫院、道路等）被認為是項目的主要輸出。
首先我們加載 OpenStreetMap 包，我們將從中獲取我們的底圖。
然后，我們創建對象映射，我們使用 OpenStreetMap 包（文檔）中的函數 openmap() 定義它。我們提供以下內容：
- upperLeft 和 lowerRight 兩個坐標對指定底圖切片的范圍 在這種情況下，我們從 linelist 行中放入了 max 和 min，因此地圖將動態響應數據
- zoom = (如果為 null 它自動確定）
- type = 底圖的類型 - 我們列出了幾個 這里的可能性和代碼當前使用第一個（[1]）“osm”
- mergeTiles = 我們選擇了TRUE，所以底圖都合并為一個
如果我們現在使用 OpenStreetMap 包中的 autoplot.OpenStreetMap() 繪制此底圖，您會看到軸上的單位不是緯度/經度坐標。它使用不同的坐標系。要正確顯示案例住所（以緯度/經度存儲），必須更改此設置。
autoplot.OpenStreetMap(map)因此，我們想使用 OpenStreetMap 包中的 openproj() 函數將地圖轉換為緯度/經度。我們提供底圖地圖，還提供我們想要的坐標參考系 (CRS)。我們通過為 WGS 1984 投影提供“proj.4”字符串來做到這一點，但您也可以通過其他方式提供 CRS。 （請參閱此頁面以更好地了解 proj.4 字符串是什么）
# Projection WGS84 map_latlon <- openproj(map, projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")現在，當我們創建繪圖時，我們看到沿軸是緯度和經度坐標。坐標系已轉換。現在，如果重疊，我們的案例將正確繪制！
# Plot map. Must use "autoplot" in order to work with ggplot autoplot.OpenStreetMap(map_latlon)7.等高線密度熱圖 Contoured density heatmaps
下面我們描述了如何在底圖上實現案例的等高線密度熱圖，從一個線列表（每個案例一行）開始。
當我們有一個具有緯度/經度坐標的底圖時，我們可以使用其居住地的經度/緯度坐標將我們的案例繪制在頂部。基于函數 autoplot。 OpenStreetMap() 創建底圖，ggplot2 函數將很容易添加到頂部，如下所示的 geom_point()：
# Plot map. Must be autoplotted to work with ggplot autoplot.OpenStreetMap(map_latlon)+ # begin with the basemapgeom_point( # add xy points from linelist lon and lat columns data = linelist, aes(x = lon, y = lat),size = 1, alpha = 0.5,show.legend = FALSE) + # drop legend entirelylabs(x = "Longitude", # titles & labelsy = "Latitude",title = "Cumulative cases")上面的地圖可能難以解釋，尤其是在點重疊的情況下。因此，您可以改為使用 ggplot2 函數 stat_density_2d() 繪制 2d 密度圖。您仍在使用 linelist 緯度/經度坐標，但會執行 2D 內核密度估計，并使用等高線顯示結果 - 就像地形圖一樣。在此處閱讀完整文檔。
# begin with the basemap autoplot.OpenStreetMap(map_latlon)+# add the density plotggplot2::stat_density_2d(data = linelist,aes(x = lon,y = lat,fill = ..level..,alpha = ..level..),bins = 10,geom = "polygon",contour_var = "count",show.legend = F) + # specify color scalescale_fill_gradient(low = "black", high = "red")+# labels labs(x = "Longitude",y = "Latitude",title = "Distribution of cumulative cases")7.2 時間序列熱圖 Time series heatmap
上面的密度熱圖顯示了累積案例cumulative cases。我們可以通過基于癥狀發作月份的熱圖來檢查隨著時間和空間的爆發，從線條列表中派生出來。
我們從linelist開始，創建一個帶有發病年份和月份的新列。 base R 中的 format() 函數改變了日期的顯示方式。在這種情況下，我們想要“YYYY-MM”。
linelist = linelist %>% mutate(date_onset_ym = format(date_onset, "%Y-%m")) # Examine the values table(linelist$date_onset_ym, useNA = "always") ## ## 2014-04 2014-05 2014-06 2014-07 2014-08 2014-09 2014-10 2014-11 2014-12 2015-01 2015-02 2015-03 2015-04 <NA> ## 3 7 13 38 77 181 204 134 100 73 54 42 32 42現在，我們簡單地通過 ggplot2 將 facetting 引入密度熱圖。應用 facet_wrap()，將新列用作行。為了清楚起見，我們將分面列的數量設置為 3。
# packages pacman::p_load(OpenStreetMap, tidyverse)# begin with the basemap autoplot.OpenStreetMap(map_latlon)+# add the density plotggplot2::stat_density_2d(data = linelist,aes(x = lon,y = lat,fill = ..level..,alpha = ..level..),bins = 10,geom = "polygon",contour_var = "count",show.legend = F) + # specify color scalescale_fill_gradient(low = "black", high = "red")+# labels labs(x = "Longitude",y = "Latitude",title = "Distribution of cumulative cases over time")+# facet the plot by month-year of onsetfacet_wrap(~ date_onset_ym, ncol = 4)Reference
the handbook team. (2021). R for applied epidemiology and public health | The Epidemiologist R Handbook. Epirhandbook.com. https://epirhandbook.com/en/index.html
tmap: get started! (2020). R-Project.org. https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html
Geometric binary predicates on pairs of simple feature geometry sets — geos_binary_pred. (2017). Github.io. https://r-spatial.github.io/sf/reference/geos_binary_pred.html
Visualise sf objects — CoordSf. (2022). Tidyverse.org. https://ggplot2.tidyverse.org/reference/ggsf.html
?
總結
以上是生活随笔為你收集整理的空间数据可视化地图绘制R语言可复现的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 强烈给大家推荐一款简单好用免费的甘特图项
- 下一篇： 《炬丰科技-半导体工艺》 玻璃薄化蚀刻