Нанесение на географическую карту происхождений наших пациентов - PullRequest
Новая
       23

Нанесение на географическую карту происхождений наших пациентов

1 голос
/ 03 февраля 2020

Я пытаюсь нанести на итальянскую географическую карту точку, показывающую происхождение («провинция») наших пациентов. В идеале размер точек должен быть пропорционален количеству пациентов из этой «провинции». Пример списка, который я хотел бы построить, следующий. 

MI  8319
CO  537
MB  436
VA  338
BG  310
PV  254
CR  244
NO  210
RM  189
CS  179

В первом столбце указан код провинции: MI (Милан), CO (Комо), MB (Монца-Брианца) и др. c. Во втором столбце указано количество пациентов из этой «провинции». Таким образом, на выходе должна быть итальянская политическая карта, где самая большая точка находится вокруг города Милан (Мичиган), вторая самая большая точка находится около города Комо (Колорадо), третья - вокруг города Монца-Брианца (МБ). ), и др c. Есть ли какой-нибудь пакет, который мог бы сделать сюжет, который я ищу? Я нашел инструмент, который мог бы выполнить эту работу здесь, но, видимо, они ожидают, что я загружу географические координаты для построения графика.

https://www.littlemissdata.com/blog/maps

Спасибо заранее.

Ответы [ 2 ]

1 голос
/ 03 февраля 2020

Вот один из способов справиться с вашей задачей. У вас есть сокращения для итальянской провинции. Вы хотите использовать их для объединения ваших данных с данными многоугольника. Если вы загружаете полигоны Италии из GADM, вы можете получить данные, которые содержат сокращения. В частности, столбец HASC_2 является единственным. Вам необходимо объединить ваши данные с данными многоугольника. Затем вы хотите создать другой набор данных, который содержит центроид. Вы можете нарисовать карту с двумя наборами данных. 

library(tidyverse)
library(sf)
library(ggthemes)

# Get the sf file from https://gadm.org/download_country_v3.html
# and import it in R.

mysf <- readRDS("gadm36_ITA_2_sf.rds")

# This is your data, which is called mydata.
mydata <- structure(list(abbs = c("MI", "CO", "MB", "VA", "BG", "PV", "CR", 
"NO", "RM", "CS"), value = c(8319L, 537L, 436L, 338L, 310L, 254L, 
244L, 210L, 189L, 179L)), class = "data.frame", row.names = c(NA, 
-10L))

   abbs value
1    MI  8319
2    CO   537
3    MB   436
4    VA   338
5    BG   310
6    PV   254
7    CR   244
8    NO   210
9    RM   189
10   CS   179

# Abbreviations are in HASC_2 in mysf. Manipulate strings so that
# I can join mydata with mysf with the abbreviations. I also get
# longitude and latitude with st_centroid(). This data set is for
# geom_point().

mysf2 <- mutate(mysf, HASC_2 = sub(x = HASC_2, pattern = "^IT.", replacement = "")) %>% 
         left_join(mydata, by = c("HASC_2" = "abbs")) %>% 
         mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
                lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

# Draw a map

ggplot() +
geom_sf(data = mysf) +
geom_point(data = mysf2, aes(x = lon, y = lat, size = value)) +
theme_map()

enter image description here

0 голосов
/ 06 февраля 2020

ОБНОВЛЕНИЕ НА КАРТЕ ВСТАВКИ

Это обновление, сделанное после другого предложения по использованию карт-вставок, которое, я думаю, будет лучшим решением для ваших вопросов и комментариев: 

library(sf)
library(cartography)


EU = st_read("~/R/mapslib/EUROSTAT/NUTS_RG_03M_2016_3035_LEVL_3.geojson")
IT = subset(EU, CNTR_CODE == "IT")
mydata <-
  structure(list(
    abbs = c("MI", "CO", "MB", "VA", "BG", "PV", "CR",
             "NO", "RM", "CS"),
    value = c(8319L, 537L, 436L, 338L, 310L, 254L,
              244L, 210L, 189L, 179L),
    nuts = c("ITC4C","ITC42","ITC4D","ITC41",
             "ITC46", "ITC48","ITC4A","ITC15",
             "ITI43","ITF61")
  ),
  class = "data.frame",
  row.names = c(NA, -10L))

patients = merge(IT, mydata, by.x = "id", by.y = "nuts")

#Get breaks for map
br=getBreaks(patients$value)
#Delimit zone
#Based on NUTS1, Nortwest Italy
par(mar=c(0,0,0,0))
ghostLayer(IT[grep("ITC",IT$NUTS_ID),], bg="lightblue")
plot(st_geometry(EU), col="grey90", add=TRUE)
plot(st_geometry(IT), col = "#FEFEE9", border = "#646464", add=TRUE)
choroLayer(
  patients,
  var = "value",
  breaks = br,
  col = carto.pal(pal1 = "red.pal", n1 = length(br)-1),
  legend.pos = "topleft",
  legend.title.txt = "Total patients",
  add = TRUE,
  legend.frame = TRUE
)
labelLayer(patients,txt="abbs", halo=TRUE, overlap = FALSE)

#Inset
par(
  fig = c(0, 0.4, 0.01, 0.4),
  new = TRUE
)
inset=patients[patients$abbs %in% c("RM","CS"),]
ghostLayer(inset, bg="lightblue")
plot(st_geometry(EU), col="grey90", add=TRUE)
plot(st_geometry(IT), col = "#FEFEE9", border = "#646464", add=TRUE)
choroLayer(
  patients,
  var = "value",
  breaks = br,
  col = carto.pal(pal1 = "red.pal", n1 = length(br)-1),
  legend.pos = "n",
  add = TRUE
)
labelLayer(patients,txt="abbs", halo=TRUE, overlap = FALSE)
box(which = "figure", lwd = 1)

enter image description here

#RESTORE PLOT

par(fig=c(0,1,0,1))

СТАРЫЙ ОТВЕТ

После моего комментария к нанесению меток, возможно, с кружками не является лучший вариант для вашей карты, учитывая концентрацию. Я предлагаю вам использовать другой вид карты для этого, так как в качестве хориста я использовал { ссылка } для dataframe. 

    library(sf)
    library(cartography)

    EU = st_read("~/R/mapslib/EUROSTAT/NUTS_RG_03M_2016_3035_LEVL_3.geojson")
    IT = subset(EU, CNTR_CODE == "IT")
    mydata <-
      structure(list(
        abbs = c("MI", "CO", "MB", "VA", "BG", "PV", "CR",
                 "NO", "RM", "CS"),
        value = c(8319L, 537L, 436L, 338L, 310L, 254L,
                  244L, 210L, 189L, 179L),
        nuts = c("ITC4C","ITC42","ITC4D","ITC41",
                 "ITC46", "ITC48","ITC4A","ITC15",
                 "ITI43","ITF61")
      ),
      class = "data.frame",
      row.names = c(NA, -10L))

    patients = merge(IT, mydata, by.x = "id", by.y = "nuts")

    #Options1 - With circles
    par(mar = c(0, 0, 0, 0))
    plot(st_geometry(IT), col = "#FEFEE9", border = "#646464")
    propSymbolsLayer(
      x = patients,
      var = "value",
      col = carto.pal(pal1 = "red.pal", n1 = 6),
      legend.title.txt = "Total patients",
      add = TRUE
    )

    #Option 2 - Chorolayer with labels
    par(mar = c(0, 0, 0, 0))
    plot(st_geometry(IT), col = "#FEFEE9", border = "#646464")
    choroLayer(
      patients,
      var = "value",
      col = carto.pal(pal1 = "red.pal", n1 = 6),
      legend.title.txt = "Total patients",
      add = TRUE
    )
    #Create labels
    patients$label = paste(patients$abbs, patients$value, sep = " - ")
    labelLayer(
      patients,
      txt = "label",
      overlap = FALSE,
      halo = TRUE,
      show.lines = TRUE,
      )

Options1 Options2

Данные от https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/nuts-2016-files.html

...