28 GIS cơ bản

28.1 Tổng quan

Các khía cạnh không gian trong dữ liệu của bạn có thể cung cấp nhiều thông tin chi tiết về tình hình đợt bùng phát dịch và để trả lời các câu hỏi như:

  • Các điểm nóng về dịch bệnh hiện nay ở đâu?
  • Các điểm nóng đã thay đổi như thế nào theo thời gian?
  • Việc tiếp cận các cơ sở y tế như thế nào? Có cần thêm sự tăng cường nào không?

Trọng tâm hiện tại trong chương GIS này nhằm giải quyết nhu cầu của các nhà dịch tễ học ứng dụng trong ứng phó với các đợt bùng phát dịch. Chúng ta sẽ khám phá các phương pháp trực quan hóa dữ liệu không gian cơ bản bằng cách sử dụng package tmapggplot2. Chúng ta cũng sẽ đi qua một số phương pháp quản lý và truy vấn dữ liệu không gian cơ bản với package sf. Cuối cùng, chúng ta sẽ đề cập ngắn gọn đến các khái niệm về thống kê không gian (spatial statistics) như mối quan hệ không gian, tự tương quan không gian và hồi quy không gian bằng cách sử dụng package spdep.

28.2 Các thuật ngữ chính

Dưới đây chúng tôi giới thiệu một số thuật ngữ chính. Để được giới thiệu kỹ lưỡng về GIS và phân tích không gian, chúng tôi khuyên bạn nên xem lại một trong các hướng dẫn hoặc khóa học dài hơn được liệt kê trong phần Tài nguyên học liệu.

Hệ thống thông tin địa lý (Geographic Information System - GIS) - GIS là một framework hoặc môi trường để thu thập, quản lý, phân tích và trực quan hóa dữ liệu không gian.

Phần mềm GIS

Một số phần mềm GIS phổ biến cho phép tương tác bằng các thao tác chuột để phát triển bản đồ và phân tích không gian. Những công cụ này đi kèm với những lợi thế như không cần phải học code và dễ dàng lựa chọn và đặt các biểu tượng và tính năng trên bản đồ theo cách thủ công. Dưới đây là hai phần mềm phổ biến:

ArcGIS - Một phần mềm GIS thương mại do công ty ESRI phát triển, rất phổ biến nhưng khá đắt

QGIS - Một phần mềm GIS mã nguồn mở nhưng làm được hầu hết mọi thứ mà ArcGIS có thể làm được. Bạn có thể tải QGIS tại đây

Sử dụng R để thao tác GIS thoạt đầu có vẻ đáng sợ hơn bởi vì thay bằng các “thao tác chuột”, nó có “giao diện dòng lệnh” (bạn phải viết code để có được kết quả mong muốn). Tuy nhiên, đây là một lợi thế lớn nếu bạn cần tạo bản đồ lặp đi lặp lại hoặc tạo các phân tích có thể tái lập được.

Dữ liệu không gian

Hai dạng dữ liệu không gian chính được sử dụng trong GIS là dữ liệu vectơ và raster:

Vector Data - Định dạng phổ biến nhất của dữ liệu không gian được sử dụng trong GIS, dữ liệu vectơ bao gồm các đặc điểm hình học của vertices and paths. Dữ liệu không gian dạng vectơ có thể được chia thành ba loại nhỏ được sử dụng rộng rãi:

  • Điểm (Points) - Một điểm bao gồm một cặp tọa độ (x, y) đại diện cho một vị trí cụ thể trong một hệ tọa độ. Điểm là dạng dữ liệu không gian cơ bản nhất và có thể được sử dụng để biểu thị một trường hợp (vd: nhà bệnh nhân) hoặc một vị trí (vd: bệnh viện) trên bản đồ.

  • Đường (Lines) - Một đường bao gồm hai điểm kết nối với nhau. Đường có độ dài và có thể được sử dụng để biểu thị những thứ như con đường hoặc sông.

  • Đa giác (Polygons) - A polygon is composed of at least three line segments connected by points. Polygon features have a length (i.e. the perimeter of the area) as well as an area measurement. Polygons may be used to note an area (i.e. a village) or a structure (i.e. the actual area of a hospital).

Raster Data - Một đa giác bao gồm ít nhất ba đoạn thẳng được nối với nhau bằng các điểm. Các đối tượng đa giác có chiều dài (vd: chu vi của khu vực) cũng như số đo diện tích. Đa giác có thể được sử dụng để biểu diễn một khu vực (vd: một ngôi làng) hoặc một cấu trúc (vd: diện tích thực tế của một bệnh viện).

Trực quan hóa dữ liệu không gian

Để thể hiện trực quan dữ liệu không gian trên bản đồ, phần mềm GIS yêu cầu bạn cung cấp đầy đủ thông tin về vị trí của các đối tượng địa lý khác nhau, trong mối quan hệ của chúng với nhau. Nếu bạn đang sử dụng dữ liệu vectơ, điều này sẽ đúng cho hầu hết các trường hợp sử dụng, thông tin này thường sẽ được lưu trữ trong một shapefile:

Shapefiles - Shapefile là một định dạng dữ liệu phổ biến để lưu trữ dữ liệu không gian “vectơ” bao gồm hoặc đường, điểm hoặc đa giác. Một shapefile thực chất là một tập hợp của ít nhất ba tệp - .shp, .shx và .dbf. Tất cả các tệp thành phần phụ này phải nằm trong cùng một thư mục để chúng có thể đọc được. Các tệp liên quan này có thể được nén vào một thư mục ZIP để gửi qua email hoặc tải xuống từ một trang web.

Shapefile sẽ chứa thông tin về bản thân các đối tượng địa lý cũng như vị trí định vị chúng trên bề mặt Trái đất. Điều này rất quan trọng bởi vì trong khi Trái đất là một quả địa cầu, các bản đồ thường là hai chiều; các lựa chọn về cách “làm phẳng” dữ liệu không gian có thể có tác động lớn đến giao diện và cách giải thích các kết quả..

Hệ trục tọa độ tham chiếu (Coordinate Reference Systems - CRS) - CRS là một hệ thống dựa trên tọa độ được sử dụng để xác định vị trí các đối tượng địa lý trên bề mặt Trái đất. Nó có một số thành phần chính:

  • Hệ tọa độ - Có nhiều hệ tọa độ khác nhau, vì vậy hãy đảm bảo rằng bạn biết hệ tọa độ của mình là từ hệ nào. Các kinh độ/vĩ độ là phổ biến nhẩt, nhưng bạn cũng có thể gặp hệ tọa độ UTM.

  • Đơn vị - Các đơn vị dành cho hệ tọa độ của bạn (ví dụ: độ thập phân, mét)

  • Dữ liệu - Một phiên bản được mô hình hóa cụ thể của Trái đất. Chúng đã được sửa đổi trong nhiều năm, vì vậy hãy đảm bảo rằng các lớp bản đồ của bạn đang sử dụng cùng một dữ liệu..

  • Phép chiếu - Tham chiếu đến phương trình toán học được sử dụng để chiếu trái đất hình tròn lên một bề mặt phẳng (bản đồ).

Hãy nhớ rằng bạn có thể tóm tắt dữ liệu không gian mà không cần sử dụng các công cụ lập bản đồ được hiển thị bên dưới. Đôi khi một bảng đơn giản phân chia theo địa lý (ví dụ: huyện, quốc gia, v.v.) là tất cả những gì cần thiết!

28.3 Bắt đầu với GIS

Có một số thứ quan trọng bạn sẽ cần phải có và suy nghĩ tới khi tạo bản đồ. Bao gồm:

  • Một tập dữ liệu – tập dữ liệu này có thể ở định dạng dữ liệu không gian (chẳng hạn như shapefiles, như đã lưu ý ở trên) hoặc nó có thể không ở định dạng không gian (chẳng hạn như csv).

  • Nếu tập dữ liệu của bạn không ở định dạng không gian, bạn cũng sẽ cần một tập dữ liệu tham chiếu. Dữ liệu tham chiếu bao gồm biểu diễn không gian của dữ liệu và các thuộc tính liên quan, sẽ bao gồm tài liệu chứa thông tin vị trí và địa chỉ của các đối tượng địa lý cụ thể.

    • Nếu bạn đang làm việc với các ranh giới địa lý được xác định trước (ví dụ: khu vực hành chính), các shapefiles tham chiếu thường có sẵn miễn phí để tải xuống từ cơ quan chính phủ hoặc tổ chức chia sẻ dữ liệu. Khi nghi ngờ, bạn có thể tìm kiếm trên Google cụm từ “[tên vùng] shapefile”

    • Nếu bạn có thông tin địa chỉ, nhưng không có vĩ độ và kinh độ, bạn có thể cần sử dụng một công cụ mã hóa địa lý để lấy dữ liệu không gian tham chiếu cho bản ghi của mình.

  • Ý tưởng về cách bạn muốn trình bày thông tin trong bộ dữ liệu của mình cho đối tượng mục tiêu. Có nhiều loại bản đồ khác nhau, và điều quan trọng là phải suy nghĩ xem loại bản đồ nào phù hợp nhất với nhu cầu của bạn.

Các loại bản đồ để trực quan hóa dữ liệu

Bản đồ Choropleth - một loại bản đồ trong đó màu sắc, sự đổ bóng, hoặc các họa tiết được sử dụng để thể hiện các vùng địa lý liên quan đến giá trị của một thuộc tính. Ví dụ: giá trị lớn hơn có thể được biểu thị bằng màu tối hơn giá trị nhỏ hơn. Loại bản đồ này đặc biệt hữu ích khi trực quan hóa một biến số và xem cách nó thay đổi trên các vùng hoặc khu vực địa chính trị xác định.

Bản đồ nhiệt mật độ trường hợp - một loại bản đồ trong đó màu sắc được sử dụng để thể hiện cường độ của một giá trị, tuy nhiên, nó không sử dụng các vùng hoặc ranh giới địa chính trị xác định để nhóm dữ liệu. Loại bản đồ này thường được sử dụng để hiển thị ‘điểm nóng’ hoặc các khu vực có mật độ hoặc tập trung điểm cao.

Bản đồ mật độ điểm - một loại bản đồ sử dụng các điểm để biểu thị các giá trị thuộc tính trong dữ liệu của bạn. Loại bản đồ này được sử dụng tốt nhất để trực quan hóa phân tán dữ liệu của bạn và scan một cách trực quan các cụm.

Bản đồ ký hiệu tỷ lệ (bản đồ ký hiệu chia độ) - một loại bản đồ tương tự như bản đồ choropleth, nhưng thay vì sử dụng màu sắc để biểu thị giá trị của một thuộc tính, nó sử dụng một ký hiệu (thường là một vòng tròn) liên quan đến giá trị. Ví dụ, một giá trị lớn hơn có thể được biểu thị bằng một ký hiệu lớn hơn một giá trị nhỏ hơn. Loại bản đồ này được sử dụng tốt nhất khi bạn muốn trực quan hóa kích thước hoặc số lượng dữ liệu của mình trên các vùng địa lý.

Bạn cũng có thể kết hợp một số loại trực quan hóa khác nhau để hiển thị các trường hợp địa lý phức tạp. Ví dụ, các trường hợp (điểm) trong bản đồ dưới đây được tô màu theo cơ sở y tế gần nhất của họ (xem phần chú thích). Các vòng tròn lớn hiển thị cơ sở y tế gần nhất trong khu vực ở một bán kính xác định, và các điểm màu đỏ tươi là những cơ sở y tế không nằm trong bất kỳ bán kính nào:

Lưu ý: Trọng tâm chính trong chương GIS này được dựa trên tình huống đáp ứng các vụ dịch tại thực địa. Do đó, nội dung của chương sẽ bao gồm các thao tác, hình ảnh hóa và phân tích dữ liệu không gian cơ bản.

28.4 Chuẩn bị

Gọi package

Đoạn code này hiển thị việc gọi các package cần thiết cho các phân tích. Trong cuốn sách này, chúng tôi nhấn mạnh việc sử dụng hàm p_load() từ package pacman, giúp cài đặt các package cần thiết gọi chúng ra để sử dụng. Bạn cũng có thể gọi các packages đã cài đặt với hàm library() của base R. Xem thêm chương R cơ bản để có thêm thông tin về các packages trong R.

pacman::p_load(
  rio,           # to import data
  here,          # to locate files
  tidyverse,     # to clean, handle, and plot the data (includes ggplot2 package)
  sf,            # to manage spatial data using a Simple Feature format
  tmap,          # to produce simple maps, works for both interactive and static maps
  janitor,       # to clean column names
  OpenStreetMap, # to add OSM basemap in ggplot map
  spdep          # spatial statistics
  ) 

Bạn có thể xem tổng quan về tất cả các package trong R xử lý dữ liệu không gian tại CRAN “Spatial Task View”.

Dữ liệu trường hợp mẫu

Với mục đích minh họa, chúng ta sẽ làm việc với một mẫu ngẫu nhiên gồm 1000 trường hợp từ bộ dữ liệu một vụ dịch Ebola mô phỏng có tên linelist (về mặt tính toán, việc làm việc với ít trường hợp hơn sẽ dễ hiển thị hơn trong sổ tay này). Để tiện the dõi, bấm để tải dữ liệu linelist “đã được làm sạch” (dưới dạng tệp .rds).

Vì chúng ta đang lấy một mẫu ngẫu nhiên của các trường hợp, nên kết quả của bạn có thể hơi khác so với những gì được minh họa ở đây khi bạn tự chạy code của mình.

Nhập dữ liệu vào bằng hàm import() từ package rio (nó xử lý nhiều loại tệp như .xlsx, .csv, .rds - xem thêm chương Nhập xuất dữ liệu để biết thêm chi tiết).

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

Tiếp theo, chúng ta chọn một mẫu ngẫu nhiên gồm 1000 hàng bằng hàm sample() từ base R.

# generate 1000 random row numbers, from the number of rows in linelist
sample_rows <- sample(nrow(linelist), 1000)

# subset linelist to keep only the sample rows, and all columns
linelist <- linelist[sample_rows,]

Bây giờ chúng ta muốn chuyển đổi bộ dữ liệu linelist đang là một dataframe, thành một đối tượng “sf” (spatial features - các đặc tính không gian). Do bộ dữ liệu linelist có hai cột “lon” và “lat” đại diện cho kinh độ và vĩ độ của nơi cư trú của từng trường hợp, nên việc chuyển đổi này khá dễ dàng.

Chúng ta sử dụng package sf (spatial features - các đặc tính không gian) và hàm st_as_sf() của nó để tạo đối tượng mới có tên linelist_sf. Đối tượng mới này về cơ bản giống với linelist, nhưng các cột lonlat đã được chỉ định là cột tọa độ, và một hệ thống tham chiếu tọa độ (CRS) đã được gán khi các điểm được hiển thị. 4326 xác định tọa độ của chúng ta dựa trên Hệ thống trắc địa thế giới 1984 (WGS84) - là tiêu chuẩn cho tọa độ GPS.

# Create sf object
linelist_sf <- linelist %>%
     sf::st_as_sf(coords = c("lon", "lat"), crs = 4326)

Dưới đây là cách bộ dữ liệu linelist gốc được hiển thị. Trong phần minh họa này, chúng ta sẽ chỉ sử dụng cột date_onsetgeometry (được xây dựng từ các trường kinh độ và vĩ độ ở trên và là cột cuối cùng trong khung dữ liệu).

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

Shapefiles phân cấp hành chính

Sierra Leone: shapefiles phân cấp hành chính

Trước đó, chúng tôi đã tải xuống tất cả các phân cấp hành chính của Sierra Leone từ website của Humanitarian Data Exchange (HDX). Ngoài ra, bạn có thể tải xuống các dữ liệu này và tất cả các dữ liệu mẫu khác cho sổ tay này qua gói R của chúng tôi, như được giải thích trong trang Tải sách và dữ liệu.

Bây giờ chúng ta sẽ thực hiện những bước sau để lưu shapefile hành chính cấp 3 vào R:

  1. Nhập shapefile
  2. Làm sạch tên cột
  3. Lọc các hàng để chỉ giữ lại các khu vực quan tâm

Để nhập một shapefile, chúng ta sử dụng hàm read_sf() từ package sf. Nó được cung cấp đường dẫn tệp thông qua hàm here(). - trong trường hợp này tệp nằm bên trong Dự án R của bạn tại thư mục “data”, rồi tới “gis”, và thư mục con “shp”, với tệp có tên “sle_adm3.shp” (xem chương Nhập xuất dữ liệuDự án R để biết thêm chi tiết). Bạn sẽ cần cung cấp đường dẫn tệp của riêng mình.

Tiếp theo, chúng ta sử dụng hàm clean_names() từ package janitor để chuẩn hóa tên các cột của shapefile. Chúng ta cũng sử dụng hàm filter() để chỉ giữ lại các hàng có tên admin2 là “Western Area Urban” hoặc “Western Area Rural”.

# ADM3 level clean
sle_adm3 <- sle_adm3_raw %>% 
  janitor::clean_names() %>% # standardize column names
  filter(admin2name %in% c("Western Area Urban", "Western Area Rural")) # filter to keep certain areas

Dưới đây, bạn có thể thấy shapefile trông như thế nào sau khi nhập và làm sạch. Cuộn sang phải to see how there are columns hành chính cấp 0 (quốc gia), hành chính cấp 1, hành chính cấp 2, và cuối cùng là hành chính cấp 3. Mỗi cấp độ có một tên và một mã định danh duy nhất “pcode”. Pcode mở rộng với mỗi sự tăng lên của cấp hành chính, vd. SL (Sierra Leone) -> SL04 (Western) -> SL0410 (Western Area Rural) -> SL040101 (Koya Rural).

Dữ liệu dân số

Sierra Leone: Dân số theo ADM3

Những dữ liệu này có thể tải xuống từ HDX (link tại đây) hoặc thông qua R package epirhandbook như đã được giải thích trong chương này. Chúng ta cũng sử dụng hàm import() để nạp tệp .csv. Chúng tôi cũng chuyển file được nhập tới hàm clean_names() để chuẩn hóa cú pháp tên cột.

# Population by ADM3
sle_adm3_pop <- import(here("data", "gis", "population", "sle_admpop_adm3_2020.csv")) %>%
  janitor::clean_names()

Bộ dữ liệu dân số trông sẽ như bên dưới. Cuộn sang bên phải để xem các cột: dân số nam (male population), dân số nữ (female populaton), tổng dân số (total population), và dân số theo từng nhóm tuổi.

Các cơ sở y tế

Sierra Leone: Dữ liệu cơ sở y tế từ OpenStreetMap

Một lần nữa, bạn có thể tải xuống thông tin vị trí của các cơ sở y tế từ HDX tại đây hoặc thông qua hướng dẫn trong chương Tải sách và dữ liệu.

Chúng ta nhập shapefile tọa độ điểm các cơ sở với hàm read_sf(), làm sạch lại tên cột, sau đó lọc để chỉ giữ lại các điểm được gắn thẻ là “hospital”, “clinic”, hoặc “doctors”.

# OSM health facility shapefile
sle_hf <- sf::read_sf(here("data", "gis", "shp", "sle_hf.shp")) %>% 
  janitor::clean_names() %>%
  filter(amenity %in% c("hospital", "clinic", "doctors"))

Dưới đây là dataframe kết quả - Cuộn phải để xem tên cơ sở và tọa độ geometry.

28.5 Vẽ đồ thị tọa độ

Cách dễ nhất để vẽ đồ thị tọa độ X-Y (kinh độ / vĩ độ, điểm), trong trường hợp này, là vẽ chúng dưới dạng điểm trực tiếp từ đối tượng linelist_sf mà chúng ta đã tạo trong phần chuẩn bị.

Package tmap cung cấp khả năng lập bản đồ đơn giản ở cả dạng tĩnh (“plot” mode) và tương tác (“view” mode) chỉ với một vài dòng code. Cú pháp của tmap tương tự như cú pháp của ggplot2, chẳng hạn như các lệnh được thêm vào nhau bằng dấu +. Đọc thêm chi tiết trong hướng dẫn này.

  1. Thiết lập tmap mode. Trong trường hợp này chúng ta sẽ sử dụng “plot” mode để tạo ra các biểu đồ tĩnh.
tmap_mode("plot") # choose either "view" or "plot"

Dưới đây, các điểm được vẽ độc lập. Hàm tm_shape() được cung cấp bởi đối tượng linelist_sf. Chúng ta sau đó thêm các điểm thông qua hàm tm_dots(), và cụ thể kích thước và màu sắc. linelist_sf là một đối tượng sf đã được chúng ta chỉ định hai cột chứa tọa độ vĩ độ/kinh độ và hệ quy chiếu tọa độ (CRS):

# Just the cases (points)
tm_shape(linelist_sf) + tm_dots(size=0.08, col='blue')

Khi đứng một mình, các điểm không cho biết nhiều thông tin. Vì vậy, chúng ta nên lập bản đồ địa giới hành chính:

Một lần nữa chúng ta sử dụng hàm tm_shape() (xem tài liệu này) nhưng thay vì cung cấp shapefile tọa độ điểm các trường hợp, chúng ta cung cấp shapefile địa giới hành chính (đa giác).

Với đối số bbox = (bbox là viết tắt của “bounding box”), chúng ta có thể chỉ định các ranh giới tọa độ. Đầu tiên, chúng ta hiển thị bản đồ mà không có bbox, và sau đó thêm nó vào.

# Just the administrative boundaries (polygons)
tm_shape(sle_adm3) +               # admin boundaries shapefile
  tm_polygons(col = "#F7F7F7")+    # show polygons in light grey
  tm_borders(col = "#000000",      # show borders with color and line weight
             lwd = 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)) +  # corner
  tm_polygons(col = "#F7F7F7") +
  tm_borders(col = "#000000", lwd = 2) +
  tm_text("admin3name")

Và bây giờ hiển thị các điểm và đa giác cùng nhau:

# 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 map

Để đọc một bài so sánh hay về các tùy chọn vẽ bản đồ trong R, hãy xem bài đăng blog này.

28.6 Phép nối theo không gian

Bạn có thể đã quen với việc nối dữ liệu từ một tập dữ liệu này sang một tập dữ liệu khác. Một số phương pháp được thảo luận trong chương Nối dữ liệu của cuốn sổ tay này. Một phép nối theo không gian phục vụ một mục đích tương tự nhưng tận dụng các mối quan hệ không gian. Thay vì dựa vào các giá trị chung trong các cột để khớp chính xác các quan sát, bạn có thể sử dụng các mối quan hệ không gian của chúng, chẳng hạn như một đối tượng địa lý nằm trong đối tượng khác, hoặc hàng xóm gần nhất với đối tượng khác, hoặc ở trong vùng đệm của một bán kính nhất định từ đối tượng khác, v.v.

Package sf cung cấp nhiều phương thức khác nhau để nối theo không gian. Xem thêm tài liệu về phương pháp st_join và các kiểu nối theo không gian trong tài liệu tham khảo này.

Điểm trong đa giác

Gán đơn vị hành chính cho các trường hợp theo không gian

Sau đây là một câu hỏi hóc búa thú vị: linelist không chứa bất kỳ thông tin nào về các đơn vị hành chính của các trường hợp. Mặc dù lý tưởng nhất là thu thập thông tin như vậy trong giai đoạn thu thập dữ liệu ban đầu, chúng ta cũng có thể gán các đơn vị hành chính cho các trường hợp riêng lẻ dựa trên mối quan hệ không gian của chúng (tức là điểm giao với một đa giác)..

Dưới đây, chúng ta sẽ giao các vị trí (điểm) theo không gian với địa giới hành chính cấp 3 (đa giác):

  1. Bắt đầu với linelist (các điểm)
  2. Nối theo không gian tới các địa giới, thiết lập kiểu nối là “st_intersects”
  3. Sử dụng hàm select() để chỉ giữ lại một số cột địa giới hành chính mới
linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects)

Tất cả các cột từ sle_adms đã được thêm vào linelist! Mỗi trường hợp bây giờ có các cột thể hiện chi tiết cấp hành chính mà nó nằm trong đó. Trong ví dụ này, chúng ta chỉ muốn giữ lại hai trong số các cột mới (địa giới hành chính cấp 3), vì vậy chúng dùng hàm select() để chọn tên các cột cũ và chỉ hai cột bổ sung quan tâm:

linelist_adm <- linelist_sf %>%
  
  # join the administrative boundary file to the linelist, based on spatial intersection
  sf::st_join(sle_adm3, join = st_intersects) %>% 
  
  # Keep the old column names and two new admin ones of interest
  select(names(linelist_sf), admin3name, admin3pcod)

Dưới đây, chỉ với mục đích hiển thị, bạn có thể thấy mười trường hợp đầu tiên và địa giới hành chính cấp 3 của chúng đã được gắn với nhau, dựa trên vị trí giao nhau trong không gian giữa các điểm với đa giác.

# 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.27125 ymin: 8.447887 xmax: -13.20522 ymax: 8.491748
## Geodetic CRS:  WGS 84
## First 10 features:
##      case_id     admin3name admin3pcod                   geometry
## 4447  63607f Mountain Rural   SL040102 POINT (-13.21982 8.463532)
## 3160  28c791       West III   SL040208 POINT (-13.25923 8.453884)
## 3967  9c3a56 Mountain Rural   SL040102 POINT (-13.21903 8.479095)
## 5293  1cf7a4       East III   SL040205 POINT (-13.20613 8.464462)
## 3555  5b1ded         West I   SL040206 POINT (-13.24692 8.485495)
## 5316  c68c31       West III   SL040208  POINT (-13.2671 8.450587)
## 297   7f3916 Mountain Rural   SL040102 POINT (-13.21223 8.458014)
## 1609  2591c9 Mountain Rural   SL040102 POINT (-13.22016 8.463424)
## 4776  24f9d8       West III   SL040208 POINT (-13.26871 8.460778)
## 5485  e48c67       West III   SL040208 POINT (-13.26651 8.474174)

Bây giờ chúng ta có thể mô tả các trường hợp theo đơn vị hành chính - điều mà chúng ta đã không thể làm trước khi nối theo không gian!

# Make new dataframe containing counts of cases by administrative unit
case_adm3 <- linelist_adm %>%          # begin with linelist with new admin cols
  as_tibble() %>%                      # convert to tibble for better display
  group_by(admin3pcod, admin3name) %>% # group by admin unit, both by name and pcode 
  summarise(cases = n()) %>%           # summarize and count rows
  arrange(desc(cases))                 # arrange in descending order

case_adm3
## # A tibble: 10 x 3
## # Groups:   admin3pcod [10]
##    admin3pcod admin3name     cases
##    <chr>      <chr>          <int>
##  1 SL040102   Mountain Rural   307
##  2 SL040208   West III         244
##  3 SL040207   West II          163
##  4 SL040204   East II          110
##  5 SL040203   East I            49
##  6 SL040201   Central I         48
##  7 SL040206   West I            39
##  8 SL040205   East III          24
##  9 SL040202   Central II        14
## 10 <NA>       <NA>               2

Chúng ta cũng có thể tạo một biểu đồ cột về số lượng trường hợp theo đơn vị hành chính.

Trong ví dụ này, chúng ta bắt đầu hàm ggplot() với bộ dữ liệu linelist_adm, trong đó chúng ta có thể áp dụng các hàm làm việc với factor như fct_infreq() để sắp xếp các cột theo tần suất (xem chương Factors để biết thêm các mẹo).

ggplot(
    data = linelist_adm,                       # begin with linelist containing admin unit info
    mapping = aes(
      x = fct_rev(fct_infreq(admin3name))))+ # x-axis is admin units, ordered by frequency (reversed)
  geom_bar()+                                # create bars, height is number of rows
  coord_flip()+                              # flip X and Y axes for easier reading of adm units
  theme_classic()+                           # simplify background
  labs(                                      # titles and labels
    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"
  )

Hàng xóm gần nhất

Tìm cơ sở y tế gần nhất / khu vực cung cấp dịch vụ y tế

Có thể sẽ hữu ích nếu biết các cơ sở y tế nằm ở đâu trong mối liên quan đến các điểm nóng về dịch bệnh.

Chúng ta có thể sử dụng phương pháp nối st_nearest_feature trong hàm st_join() (sf package) để trực quan hóa cơ sở y tế gần nhất với các trường hợp bệnh nhân.

  1. Chúng ta bắt đầu với shapefile linelist có tên linelist_sf
  2. Chúng ta nối theo không gian với đối tượng sle_hf, chứa thông tin về vị trí của các cơ sở y tế và phòng khám (các điểm)
# Closest health facility to each case
linelist_sf_hf <- linelist_sf %>%                  # begin with linelist shapefile  
  st_join(sle_hf, join = st_nearest_feature) %>%   # data from nearest clinic joined to case data 
  select(case_id, osm_id, name, amenity) %>%       # keep columns of interest, including id, name, type, and geometry of healthcare facility
  rename("nearest_clinic" = "name")                # re-name for clarity

Chúng ta có thể thấy bên dưới (50 hàng đầu tiên) rằng mỗi trường hợp hiện đã có dữ liệu về phòng khám/bệnh viện gần nhất

Chúng ta có thể thấy rằng “Den Clinic” là cơ sở y tế gần nhất với khoảng ~30% các trường hợp.

# Count cases by health facility
hf_catchment <- linelist_sf_hf %>%   # begin with linelist including nearest clinic data
  as.data.frame() %>%                # convert from shapefile to dataframe
  count(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 order

hf_catchment                         # print to console
##                          nearest_clinic case_n
## 1                            Den Clinic    361
## 2       Shriners Hospitals for Children    325
## 3         GINER HALL COMMUNITY HOSPITAL    177
## 4                             panasonic     53
## 5                     ARAB EGYPT CLINIC     32
## 6 Princess Christian Maternity Hospital     28
## 7                  MABELL HEALTH CENTER     12
## 8                                  <NA>     12

Để trực quan hóa kết quả, chúng ta có thể sử dụng tmap - lúc này interactive mode sẽ xem dễ dàng hơn

tmap_mode("view")   # set tmap mode to interactive  

# plot the cases and clinic points 
tm_shape(linelist_sf_hf) +            # plot cases
  tm_dots(size=0.08,                  # cases colored by nearest clinic
          col='nearest_clinic') +    
tm_shape(sle_hf) +                    # plot clinic facilities in large black dots
  tm_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")

Vùng đệm

Chúng ta cũng có thể tìm hiểu xem có bao nhiêu trường hợp nằm trong khoảng cách đi bộ 2,5 km (~30 phút) từ cơ sở y tế gần nhất.

Chú ý: Để tính toán khoảng cách chính xác hơn, tốt hơn nên chiếu lại đối tượng sf của bạn lên hệ thống chiếu bản đồ địa phương tương ứng chẳng hạn như UTM (Trái đất được chiếu lên bề mặt phẳng). Trong ví dụ này, để đơn giản hơn, chúng ta sẽ dựa vào Hệ tọa độ địa lý của Hệ thống trắc địa thế giới (WGS84) (Trái đất được biểu diễn trong một bề mặt hình cầu/tròn, do đó các đơn vị được tính bằng độ thập phân). Chúng ta sẽ sử dụng quy đổi chung là: 1 độ thập phân = ~111km.

Xem thêm thông tin về phép chiếu bản đồ và hệ tọa độ tại bài viết này esri article. Blog này nói về các loại phép chiếu bản đồ khác nhau và cách người ta có thể chọn phép chiếu phù hợp tùy thuộc vào khu vực quan tâm và bối cảnh của bản đồ/phân tích của bạn.

Đầu tiên, tạo một vùng đệm hình tròn với bán kính ~2.5km xung quanh mỗi cơ sở y tế. Điều này được thực hiện với hàm st_buffer() trong package tmap. Bởi vì đơn vị của bản đồ là kinh/vĩ độ thập phân, đó là cách “0.02” được diễn giải. Nếu hệ tọa độ bản đồ của bạn tính bằng mét, thì số đó phải được cung cấp bằng mét.

sle_hf_2k <- sle_hf %>%
  st_buffer(dist=0.02)       # decimal degrees translating to approximately 2.5km 

Sau đây chúng ta vẽ biểu đồ của chính các vùng đệm, với:

tmap_mode("plot")
# Create circular buffers
tm_shape(sle_hf_2k) +
  tm_borders(col = "black", lwd = 2)+
tm_shape(sle_hf) +                    # plot clinic facilities in large red dots
  tm_dots(size=0.3, col='black')      

Thứ hai, chúng ta giao các vùng đệm này với các trường hợp (điểm) bằng cách sử dụng hàm st_join() và kiểu nối là st_intersects. Tức là, dữ liệu từ vùng đệm được nối với các điểm mà chúng giao với nhau.

# Intersect the cases with the buffers
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)

Bây giờ chúng ta có thể đếm kết quả: có nrow(linelist_sf_hf_2k[is.na(linelist_sf_hf_2k$osm_id.y),]) trường hợp trong số 1000 trường hợp không giao nhau với bất kỳ vùng đệm nào (giá trị đó bị thiếu), và do đó họ sống ở nơi cách nhiều hơn 30 phút đi bộ tới cơ sở y tế gần nhất.

# 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

Chúng ta có thể trực quan hóa kết quả sao cho các trường hợp không giao nhau với bất kỳ vùng đệm nào sẽ xuất hiện với màu đỏ.

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