29.1 冈比亚儿童疟疾

冈比亚地形

sp_path <- "data/" # 存储临时地形文件
if (!dir.exists(sp_path)) dir.create(sp_path, recursive = TRUE)
# Gambia 海拔数据
gambia_alt <- raster::getData(name = "alt", country = "GMB", mask = TRUE, path = sp_path)
# Gambia 市级行政边界数据
gambia_map <- raster::getData("GADM", country = "GMB", level = 2, path = sp_path)
# 绘制冈比亚地形
rasterVis::levelplot(gambia_alt,
  margin = FALSE,
  main = "Elevation",
  colorkey = list(
    space = "top",
    labels = list(at = seq(from = -5, to = 65, by = 10)),
    axis.line = list(col = "black")
  ),
  par.settings = list(
    axis.line = list(col = "transparent")
  ),
  scales = list(draw = FALSE),
  col.regions = hcl.colors,
  at = seq(-5, 65, len = 101)
) +
  latticeExtra::layer(sp::sp.polygons(gambia_map, lwd = 1.5))
冈比亚地形海拔数据

图 29.1: 冈比亚地形海拔数据

rgdal 包可以实现坐标变换

# 加载数据
data(gambia, package = "geoR")
# 坐标变换
library(rgdal)
sps <- SpatialPoints(gambia[, c("x", "y")],
  proj4string = CRS("+proj=utm +zone=28")
)
spst <- spTransform(sps, CRS("+proj=longlat +datum=WGS84"))
gambia[, c("x", "y")] <- coordinates(spst)
# 聚合数据
gambia_agg <- aggregate(
  formula = cbind(pos, netuse, treated) ~ x + y + green + phc,
  data = gambia, FUN = function(x) sum(x) / length(x)
)
# 抽取指定位置的海拔数据
# raster::extract(gambia_alt, gambia[, c("x", "y")])

\(Y \sim b(1,p)\) 每个人检验结果,就是感染 1 或是没有感染 0,感染率 \(p\) 的建模分析,个体水平

library(highcharter)
hchart(gambia_agg, "bubble", hcaes(x = x, y = y, fill = pos, size = pos),
  maxSize = "5%", name = "Gambia", showInLegend = FALSE
) %>%
  hc_yAxis(title = list(text = "Latitude")) %>%
  hc_xAxis(title = list(text = "Longitude"), labels = list(align = "center")) %>%
  hc_colorAxis(
    stops = color_stops(colors = hcl.colors(palette = "Plasma", n = 10))
  ) %>%
  hc_tooltip(
    pointFormat = "({point.x:.2f}, {point.y:.2f}) <br/> Size: {point.z:.2f}"
  )

图 29.2: 各个村庄疟疾流行度

# gm_data <- download_map_data("https://code.highcharts.com/mapdata/countries/gm/gm-all.js")
# get_data_from_map(gm_data)

hcmap("countries/gm/gm-all.js") %>%
  hc_title(text = "Gambia")
data("USArrests", package = "datasets")
data("usgeojson") # 加载地图数据 地图数据的结构

USArrests <- transform(USArrests, state = rownames(USArrests))

highchart() %>%
  hc_title(text = "Violent Crime Rates by US State") %>%
  hc_subtitle(text = "Source: USArrests data") %>%
  hc_add_series_map(usgeojson, USArrests,
    name = "Murder arrests (per 100,000)",
    value = "Murder", joinBy = c("woename", "state"),
    dataLabels = list(
      enabled = TRUE,
      format = "{point.properties.postalcode}"
    )
  ) %>%
  hc_colorAxis(stops = color_stops()) %>%
  hc_legend(valueDecimals = 0, valueSuffix = "%") %>%
  hc_mapNavigation(enabled = TRUE)

highcharter 包含三个数据集分别是: worldgeojson 世界地图(国家级)、 usgeojson 美国地图(州级)、 uscountygeojson 美国地图(城镇级)。其它地图数据见 https://code.highcharts.com/mapdata/

# 添加地图数据
hcmap(map = "countries/cn/custom/cn-all-sar-taiwan.js") %>%
  hc_title(text = "中国地图")
library(mapdeck)
# 多边形
mapdeck() %>%
  add_polygon(
    data = spatialwidget::widget_melbourne, 
    fill_colour = "SA2_NAME",
    palette = "spectral"
  )

mapdeck( location = c(145, -37.8), zoom = 10) %>%
  add_geojson(
    data = mapdeck::geojson
  )

googlewayggmapRgoogleMaps

library(RgoogleMaps)
lat <- c(40.702147, 40.718217, 40.711614)
lon <- c(-74.012318, -74.015794, -73.998284)
center <- c(mean(lat), mean(lon))
zoom <- min(MaxZoom(range(lat), range(lon)))
bb <- qbbox(lat, lon)
par(pty = "s")
# OSM
myMap <- GetMap(center, zoom = 15)
PlotOnStaticMap(myMap,
  lat = lat, lon = lon, pch = 20,
  col = c("red", "blue", "green"), cex = 2
)
library(ggmap)

us <- c(left = -125, bottom = 25.75, right = -67, top = 49)
get_stamenmap(us, zoom = 5, maptype = "toner-lite") %>%
  ggmap()