R语言绘制核密度估计的地图（欧洲作为底图）

核密度估计是研究地理空间上现象、事物集聚的重要方法，在QGIS和ArcGIS都有相应的工具进行操作。但是在R语言中，绘制核密度估计的地图并不那么友好。因为缺乏直观的图形操作界面，核密度估计的结果以及将之绘制在地图上都通过代码完成，这对使用者理解计算结果、绘制结果的基本过程有着更高的要求。

最近用R语言绘制地图有点多，因为一旦地图的样式确定之后，一段代码可以复用很多次。这样，就可以专心去分析、调整统计数据。

然而，现有计算核密度估计并绘制地图的中英文资源并不是非常完整，对于其中的过程也缺乏说明。也尝试过询问AI，给出的代码难以达到理想的效果，进一步反馈错误、询问细节，又往往“牛头不对马嘴”。

因此，遂自己尝试探索，将经验总结。

前言

在开始介绍R语言计算与绘制核密度估计的代码之前，首先介绍一下绘制核密度估计地图的基本原理。过程大致如下：

计算核密度估计结果，为一.im对象，可以理解为图像，主要通过spatstat包完成。
将核密度估计结果绘制在地图底图之上，地图为sf或sp对象，为矢量图形，主要通过ggplot2包完成。

代码主要是关于这两步，其中如果要得到理想的结果，还有很多的细节需要注意。

在正式开始写代码之前，需要用到的包如下：

sf，用于操作空间图形
sp，用于操作空间图形
spatstata，用于进行空间分析
ggplot2，用于绘图
dplyr，用于操作数据框
tidyr，用于操作数据框

需要预先安装一下，以免某个函数用不了。

library(sf)
library(sp)
library(spatstat)
library(ggplot2)
library(dplyr)
library(tidyr)

关于spatstat的安装
该包下面有许多依赖包，如spatstat.geom等。笔者使用的R语言版本是R 4.3，安装spatstat的时候显示依赖包安装错误，原因是获取到的依赖包版本低，不满足要求。
于是，前往https://cloud.r-project.org/一探究竟。发现收录的R 4.3的包版本确实不满足需求。那满足版本要求的依赖包在哪里呢？在R 4.4！！
因此，下载几个依赖包的.zip文件直接安装。
虽然加载的时候会出现版本警告，但使用起来毫无问题。

文末会给出完整的代码。

一、准备数据

没有数据怎么绘图呢？这里使用的数据是计算核密度估计并绘图成功的数据，因为是个人私用，不作分享。替换成实际用于绘图的数据即可。

geo <- st_read("【地理数据】data.geojson")
geo_country <- st_read("【地理数据】欧盟国家边界.geojson")

geo的地图是欧盟NUTS2区域。NUTS是欧盟空间规划单元，NUTS0等同国家一级，此下有NUTS1、NUTS2、NUTS3。NUTS2的底图可通过eurostat包得到（NUTS3在欧盟不同规划期内有调整，但NUTS2基本稳定）：

library(eurostat)

geo_nuts2 <- get_eurostat_geospatial(nuts_level = 2)

geo主要是在geo_nuts2以外另有一列number，实际含义为某产业的数量，用于对核密度估计的点进行赋值。

如果想复用接下来的代码，那么准备的数据至少包含一个空间对象，可以是点，也可以是面（geo是多边形，后面会有提取多边形中心点的代码），在基本的空间对象外，可以有一列用于加权，表示某一部分点更为重要。

geo_country实际通过eurostat包得到，后通过st_write保存为geojson文件，方便下次直接读取：

library(eurostat)

geo_country <- get_eurostat_geospatial(nuts_level = 0)

二、准备核密度估计使用的点

（一）转换原有的坐标系为UTM坐标系。因为spatstat对平面坐标系的点进行分析，如果原来的坐标系是WGS84或其他以经纬度表示的坐标，需要转换为UTM坐标系。欧洲绘图常用的UTM坐标系区号是3035。

geo_utm <- st_transform(geo, crs = 3035)

（二）填补number列的空值为0；从多边形提取中心点（如果本身是点空间对象，不需要这一步）。

points <- geo %>%
  mutate(number = replace_na(number, 0)) %>%
  st_centroid()

注意！！这里不能用filter(!is.na(number))等方法删除number列为空的行，否则会影响核密度估计的结果。

（三）将提取的点的坐标系也转换为UTM坐标系。

points_utm <- st_transform(points, crs = 3035)

（四）提取点的坐标数据，用于核密度估计。

coords <- st_coordinates(points_utm)
x <- coords[, "X"]
y <- coords[, "Y"]

这里的点即时参与计算核密度估计结果的点。

三、计算核密度估计结果

（一）准备ppp对象。ppp对象是spatstat包核密度估计方法的特定对象，可以把核密度估计的点的坐标、点权重封装到一起。

window <- as.owin(st_geometry(geo_utm))
ppp_points <- ppp(x, y, window = window, marks = points_utm$number)

参数说明：

x和y即为点的横坐标、纵坐标。

参数window，即核密度估计结果的呈现范围，需要是一个owin对象。很多资源这里用的是所有点的最大外接矩形，但是这和实际的地图情况相去甚远（底图怎么会是一个矩形！）。这里的呈现范围指定为geo空间对象，即为欧洲各国的地理边界（虽然是NUTS2一级，但会自动提取外围）：先通过st_geometry()提取空间对象，再通过as.owin()方法将空间对象转换为owin对象。

参数marks，用于接收点的权重，这里为number。

（二）核密度估计，使用density或density.ppp方法。

dens <- density(
  ppp_points,
  sigma = 100000,
  weights = marks(ppp_points),
  dimyx = 2048
)

参数说明：

ppp_points即为ppp对象。
sigma单位是米（m），为带宽或搜索半径，决定核密度估计的范围。由于这里是NUTS2区域的中心点，范围比较大，此处设定为100公里。具体实践过程中，可以多次调整，反复实验，找到最合适的。
weights为核密度估计时使用的点的权重。注意！！需要通过marks()方法从ppp对象提取，所以，在创建ppp对象时别忘了指定marks参数！
dimyx为点的数量，可以理解为分辨率，值越大，核密度估计结果的边缘越平滑。因为density()返回的结果是.im对象，为一个图形，图形的分辨率越高，越和矢量数据的底图更好匹配！同时，计算的时间也会越长。

这一步完成，可以用plot()方法查看结果（此处，因为欧洲一些国家的海外领土离欧洲大陆很远，图像把这一部分也包括了进来，所以，后续会对欧洲大陆做一个范围的裁剪）。如果出现核密度估计值为0，或呈现结果范围不对的情况，需要检查权重列的添加、底图指定为window是否操作正确。

四、绘制核密度估计结果

（一）创建用于ggplot绘图的核密度估计结果数据框。

dens_df <- as.data.frame(dens, long = TRUE)

dens_df会包括三列：x（横坐标）、y（纵坐标）、value（核密度估计值）。主要根据value列进行着色。

如果实际的点比较稀疏，使得核密度估计的结果数值比较小时（图例出现了1e-6等科学计数法表示），可以考虑进行标准化，让出图的图例更好看一些。具体如下：

dens_df$v_scaled <- scale(dens_df$value)

（二）使用ggplot绘制核密度估计结果。

ggplot() +
  geom_tile(
    data = dens_df,
    aes(x = x, y = y, fill = v_scaled)
  ) +
  scale_fill_distiller(
    palette = "Spectral",
    guide = guide_colorbar(
      barwidth = unit(0.3, "cm"),
      barheight = unit(1.5, "cm"),
    )
  ) +
  geom_sf(
    data = geo_utm,
    fill = NA,
    size = 1,
    color = "gray"
  ) +
  geom_sf(
    data = geo_country,
    fill = NA,
    size = 3,
    color = "black"
  ) +
  coord_sf(
    xlim = c(2500000, 6500000),
    ylim = c(1600000, 5200000),
    expand = FALSE
  ) +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "white", colour = NA),
    text = element_text(family = "serif"),
    legend.text = element_text(family = "serif", size = 12),
    legend.title = element_text(size = 12),
    legend.position = c(0.9, 0.9)
  )

这里分函数进行说明：

ggplot()可以理解为后续的结果都装在这里面。
geo_tile()用于绘制图像。如果用geo_raster的话，图像会和实际位置发生一些偏移。因此推荐使用这个函数。
scale_fill_distiller()指定填充颜色，可以多次试验，找到可读性最好的。
第一个geo_sf()绘制欧盟NUTS2区域的边界线。
第二个geo_sf()绘制欧洲各国的边界线。
coord_sf()限制地图绘制的范围为欧洲大陆。注意！！ggplot地图绘制会在限定的范围之外再扩展一圈，但是这一圈是没有核密度估计结果，颜色是白的！将expand参数设置为FALSE可以避免这一情况。
theme_void()指定绘图主题，去掉坐标轴、网格线，适合地图绘制。
theme()指定其他相关的参数。plot.background设置背景为白色，这样在一些图片查看器里会好看一点（不然透明的部分会干扰读图）；text设置字体为serif，即Times New Roman；legend.text设置图例字体，也是Times New Roman，注意！！该参数不会影响图例标题，图例标题字体由text参数决定；legend.title设置图例标题字体大小；legend.position设置图例位置，为右上角。

（三）保存图片。

ggsave(
  filename = paste0("【图片】核密度估计-100km.png"),
  device = "png",
  width = 4.5,
  height = 4.5,
  dpi = 300
)

结果展示：

欧洲地图绘制参考：4 Drawing maps of Europe | Using Eurostat with R

library(sf)
library(sp)
library(spatstat)
library(ggplot2)
library(dplyr)
library(tidyr)


# 读取需要的地理数据
geo <- st_read("【地理数据】data.geojson")
geo_country <- st_read("【地理数据】欧盟国家边界.geojson")


# 转变utm坐标系，为spatstat空间分析的需要
geo_utm <- st_transform(geo, crs = 3035)
# 提取底图的空间特征，以地图的空间边界限制核密度估计的范围
window <- as.owin(st_geometry(geo_utm))


# 准备核密度估计使用的点
points <- geo %>%
  # 填补缺失值为0，此处不建议删除缺失值，会影响核密度估计结果
  # number用于对点进行加权
  mutate(number = replace_na(number, 0)) %>%
  # 计算多边形区域的中心点
  st_centroid()
# 将点也转换为utm坐标系
points_utm <- st_transform(points, crs = 3035)


# 提取点的坐标数据，用于核密度估计
coords <- st_coordinates(points_utm)
x <- coords[, "X"]
y <- coords[, "Y"]


# 核密度估计需要ppp对象，在此创建
# marks用于对点进行加权
ppp_points <- ppp(x, y, window = window, marks = points_utm$number)


# 开始核密度估计
dens <- density(
  ppp_points,
  # simga指定搜索半径，单位：米
  sigma = 100000,
  # 对核密度估计的点进行加权
  weights = marks(ppp_points),
  # dens结果为一im对象，即图形，此参数决定图像分辨率
  # 分辨率越高，边界越接近真实地理边界，计算时间越长
  dimyx = 2048
)


# 创建用于ggplot绘图的数据框
dens_df <- as.data.frame(dens, long = TRUE)
# value一列即点的核密度估计结果，此处标准化使图例会好看一点
dens_df$v_scaled <- scale(dens_df$value)


# 使用ggplot绘制核密度估计结果
ggplot() +
  # 绘制核密度估计结果
  geom_tile(
    data = dens_df,
    aes(x = x, y = y, fill = v_scaled)
  ) +
  scale_fill_distiller(
    palette = "Spectral",
    guide = guide_colorbar(
      barwidth = unit(0.3, "cm"),
      barheight = unit(1.5, "cm"),
    )
  ) +
  # 叠加地理边界
  geom_sf(
    data = geo_utm,
    fill = NA,
    size = 1,
    color = "gray"
  ) +
  geom_sf(
    data = geo_country,
    fill = NA,
    size = 3,
    color = "black"
  ) +
  coord_sf(
    xlim = c(2500000, 6500000),
    ylim = c(1600000, 5200000),
    expand = FALSE
  ) +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "white", colour = NA),
    text = element_text(family = "serif"),
    legend.text = element_text(family = "serif", size = 12),
    legend.title = element_text(size = 12),
    legend.position = c(0.9, 0.9)
  )
ggsave(
  filename = paste0("【图片】核密度估计-100km.png"),
  device = "png",
  width = 4.5,
  height = 4.5,
  dpi = 300
)