Cách tính diện tích vùng trong ArcGIS

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 tmap và ggplot2. 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.

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.

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.

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].

Để 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!

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.

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.

Đ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 và 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”.

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 SL0410 [Western Area Rural] -> SL040101 [Koya Rural].

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 % 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.

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 % 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.

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.

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.

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 % # 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 % # 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 % # 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 ## ## 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 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" ]

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 % # 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 % # 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 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"]

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 % 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 % 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"]

Các giá trị tùy chọn cho đối số join bao gồm [lấy từ tài liệu này]

• 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

Bản đồ Choropleth có thể hữu ích để trực quan hóa dữ liệu của bạn theo khu vực được xác định trước, thường là đơn vị hành chính hoặc khu vực y tế. Ví dụ, trong ứng phó với ổ dịch, điều này có thể giúp xác định mục tiêu phân bổ nguồn lực cho các khu vực cụ thể có tỷ lệ mắc bệnh cao.

Bây giờ chúng ta đã gán tên đơn vị hành chính cho tất cả các trường hợp [xem phần về phép nối không gian, ở trên], chúng ta có thể bắt đầu lập bản đồ số lượng trường hợp theo khu vực [bản đồ choropleth].

Vì chúng ta cũng có dữ liệu dân số theo cấp hành chính cấp 3 [ADM3], chúng ta có thể thêm thông tin này vào bảng case_adm3 đã được tạo trước đó.

Chúng ta bắt đầu với dataframe case_adm3 được tạo trong bước trước đó, là bảng tóm tắt của từng đơn vị hành chính và số lượng các trường hợp của nó.

1. Dữ liệu dân số sle_adm3_pop được nối sử dụng hàm left_join[] từ dplyr trên cơ sở các giá trị chung trên cột admin3pcod trong bộ dữ liệu case_adm3, và cột adm_pcode trong bộ dữ liệu sle_adm3_pop. Xem chương Nối dữ liệu].
   
2. select[] được áp dụng trên dataframe mới, để chỉ giữ những cột cần thiệt - total nghĩa là tổng dân số
   
3. Các trường hợp trên 10,000 dân được tính toán bằng cách tạo cột mới với hàm mutate[]

# Add population data and calculate cases per 10K population case_adm3 % left_join[sle_adm3_pop, # add columns from pop dataset by = c["admin3pcod" = "adm3_pcode"]] %>% # join based on common values across these two columns select[names[case_adm3], total] %>% # keep only important columns, including total population mutate[case_10kpop = round[cases/total * 10000, 3]] # make new column with case rate per 10000, rounded to 3 decimals case_adm3 # print to console for viewing

## # A tibble: 10 x 5 ## # Groups: admin3pcod [10] ## admin3pcod admin3name cases total case_10kpop ## ## 1 SL040102 Mountain Rural 307 33993 90.3 ## 2 SL040208 West III 244 210252 11.6 ## 3 SL040207 West II 163 145109 11.2 ## 4 SL040204 East II 110 99821 11.0 ## 5 SL040203 East I 49 68284 7.18 ## 6 SL040201 Central I 48 69683 6.89 ## 7 SL040206 West I 39 60186 6.48 ## 8 SL040205 East III 24 500134 0.48 ## 9 SL040202 Central II 14 23874 5.86 ## 10 2 NA NA

Nối bảng này với đa giác shapfile ADM3 để vẽ bản đồ

case_adm3_sf % # begin with cases & rate by admin unit left_join[sle_adm3, by="admin3pcod"] %>% # join to shapefile data by common column select[objectid, admin3pcod, # keep only certain columns of interest admin3name = admin3name.x, # clean name of one column admin2name, admin1name, cases, total, case_10kpop, geometry] %>% # keep geometry so polygons can be plotted drop_na[objectid] %>% # drop any empty rows st_as_sf[] # convert to shapefile

Vẽ bản đồ kết quả

# tmap mode tmap_mode["plot"] # view static map # plot polygons tm_shape[case_adm3_sf] + tm_polygons["cases"] + # color by number of cases column tm_text["admin3name"] # name display

Chúng ta cũng có thể lập bản đồ tỷ suất mới mắc

# Cases per 10K population tmap_mode["plot"] # static viewing mode # plot tm_shape[case_adm3_sf] + # plot polygons tm_polygons["case_10kpop", # color by column containing case rate breaks=c[0, 10, 50, 100], # define break points for colors palette = "Purples" # use a purple color palette ] + tm_text["admin3name"] # display text

Nếu bạn đã quen với việc sử dụng ggplot2, bạn có thể sử dụng package này để vẽ bản đồ tĩnh cho dữ liệu của bạn. Hàm geom_sf[] sẽ vẽ các đối tượng khác nhau dựa trên các đối tượng địa lý [điểm, đường thẳng hoặc đa giác] có trong dữ liệu của bạn. Ví dụ: bạn có thể sử dụng hàm geom_sf[] trong ggplot[] bằng cách sử dụng dữ liệu sf với dạng hình học đa giác để tạo bản đồ choropleth.

Để minh họa cách thức hoạt động của nó, Chúng ta có thể bắt đầu với shapefile đa giác ADM3 mà chúng ta đã sử dụng lúc trước. Xin nhớ lại rằng đây là các khu vực hành chính cấp 3 ở Sierra Leone:

## 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 ## # A tibble: 12 x 20 ## objectid admin3name admin3pcod admin3ref_n admin2name admin2pcod admin1name admin1pcod admin0name admin0pcod date valid_on valid_to ## * ## 1 155 Koya Rural SL040101 Koya Rural Western Are~ SL0401 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 2 156 Mountain Rural SL040102 Mountain Rural Western Are~ SL0401 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 3 157 Waterloo Rural SL040103 Waterloo Rural Western Are~ SL0401 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 4 158 York Rural SL040104 York Rural Western Are~ SL0401 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 5 159 Central I SL040201 Central I Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 6 160 East I SL040203 East I Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 7 161 East II SL040204 East II Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 8 162 Central II SL040202 Central II Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 9 163 West III SL040208 West III Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 10 164 West I SL040206 West I Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 11 165 West II SL040207 West II Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## 12 167 East III SL040205 East III Western Are~ SL0402 Western SL04 Sierra Le~ SL 2016-08-01 2016-10-17 NA ## # ... with 7 more variables: shape_leng , shape_area , rowcacode0 , rowcacode1 , rowcacode2 , rowcacode3 , ## # geometry

Chúng ta có thể sử dụng hàm left_join[] từ dplyr để thêm dữ liệu mà chúng ta muốn vẽ tới đối tượng shapefile. Trong trường hợp này, chúng ta sẽ sử dụng bộ dữ liệu case_adm3 mà chúng ta đã tạo trước đó để tóm tắt số lượng trường hợp theo khu vực hành chính; tuy nhiên, chúng ta có thể sử dụng phương pháp tương tự này để vẽ bất kỳ dữ liệu nào được lưu trữ trong data frame.

sle_adm3_dat % inner_join[case_adm3, by = "admin3pcod"] # inner join = retain only if in both data objects select[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 ## # A tibble: 9 x 3 ## admin3name.x cases geometry ## ## 1 Mountain Rural 307 [[[-13.21496 8.474341, -13.21479 8.474289, -13.21465 8.474296, -13.21455 8.474298, -... ## 2 Central I 48 [[[-13.22646 8.489716, -13.22648 8.48955, -13.22644 8.489513, -13.22663 8.489229, -1... ## 3 East I 49 [[[-13.2129 8.494033, -13.21076 8.494026, -13.21013 8.494041, -13.2096 8.494025, -13... ## 4 East II 110 [[[-13.22653 8.491883, -13.22647 8.491853, -13.22642 8.49186, -13.22633 8.491814, -1... ## 5 Central II 14 [[[-13.23154 8.491768, -13.23141 8.491566, -13.23144 8.49146, -13.23131 8.491294, -1... ## 6 West III 244 [[[-13.28529 8.497354, -13.28456 8.496497, -13.28403 8.49621, -13.28338 8.496086, -1... ## 7 West I 39 [[[-13.24677 8.493453, -13.24669 8.493285, -13.2464 8.493132, -13.24627 8.493131, -1... ## 8 West II 163 [[[-13.25698 8.485518, -13.25685 8.485501, -13.25668 8.485505, -13.25657 8.485504, -... ## 9 East III 24 [[[-13.20465 8.485758, -13.20461 8.485698, -13.20449 8.485757, -13.20431 8.485577, -...

Để tạo biểu đồ cột về số lượng trường hợp theo khu vực, sử dụng ggplot2, sau đó chúng ta có thể gọi geom_col[] như sau:

ggplot[data=sle_adm3_dat] + geom_col[aes[x=fct_reorder[admin3name.x, cases, .desc=T], # reorder x axis by descending 'cases' y=cases]] + # y axis is number of cases by region theme_bw[] + labs[ # set figure text title="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

Nếu chúng ta muốn sử dụng ggplot2 để tạo bản đồ choropleth về số lượng trường hợp, chúng ta có thể sử dụng cú pháp tương tự để gọi hàm geom_sf[]:

ggplot[data=sle_adm3_dat] + geom_sf[aes[fill=cases]] # set fill to vary by case count variable

Sau đó, chúng ta có thể tùy chỉnh hình thức bản đồ bằng cách sử dụng ngữ pháp nhất quán trên ggplot2, ví dụ:

ggplot[data=sle_adm3_dat] + geom_sf[aes[fill=cases]] + scale_fill_continuous[high="#54278f", low="#f2f0f7"] + # change color gradient theme_bw[] + labs[title = "Number of cases, by administrative unit", # set figure text subtitle = "Admin level 3" ]

Đối với người dùng R cảm thấy thoải mái khi làm việc với ggplot2, geom_sf[] cung cấp một cách làm đơn giản và trực tiếp, phù hợp với việc trực quan hóa bản đồ cơ bản. Để tìm hiểu thêm, hãy đọc hướng dẫn geom_sf[] này hoặc sách về ggplot2.

Dưới đây, chúng tôi mô tả cách lấy được bản đồ cơ sở cho bản đồ ggplot2 bằng cách sử dụng các tính năng của OpenStreetMap. Các phương pháp thay thế bao gồm sử dụng ggmap, yêu cầu bạn đăng ký miễn phí với Google [chi tiết].

OpenStreetMap là một dự án hợp tác nhằm tạo ra một bản đồ thế giới có thể chỉnh sửa miễn phí. Dữ liệu vị trí địa lý nền tảng ví dụ: vị trí của các thành phố, đường xá, đặc điểm tự nhiên, sân bay, trường học, bệnh viện, đường xá, v.v.] được coi là đầu ra chính của dự án.

Đầu tiên, chúng ta gọi package OpenStreetMap để lấy bản đồ cơ sở.

Sau đó chúng ta tạo đối tượng map, được xác định bằng cách sử dụng hàm openmap[] từ package OpenStreetMap [tài liệu]. Chúng ta cung cấp những thông tin sau:

• upperLeft và lowerRight Hai cặp tọa độ xác định giới hạn của ô bản đồ cơ sở

  • Trong trường hợp này, chúng tôi đã đưa giá trị tối đa và tối thiểu từ các hàng trong linelist, vì vậy bản đồ sẽ tương tác động với dữ liệu
    

• zoom = [nếu null nó sẽ được xác định tự động]
  

• type = loại bản đồ cơ sở nào - chúng tôi đã liệt kê một số khả năng ở đây và code hiện đang sử dụng cái đầu tiên [[1]] “osm”
  

• mergeTiles = chúng tôi đã chọn TRUE để tất cả các lớp nền được hợp nhất thành một

# load package pacman::p_load[OpenStreetMap] # Fit basemap by range of lat/long coordinates. Choose tile type map