地理热/等值线图的空间插值的最佳方法？

Question

我想使用类似ggplot2并ggmap以产生任意值的热图,如每平方米楼价在街道水平,以在地理区域(高分辨率).

不幸的是,这项任务似乎相当困难,因为虽然ggplot2可以产生一个很大的密度图,但似乎无法在没有事先插值的情况下将这样的空间数据可视化.

为此,我使用了库akima(不规则数据的网格双变量插值)和mgcv(带有积分平滑度估计的广义加性模型),但是我对插值方法的了解最多也是平庸的,而且我能够产生的结果不是足够令人满意

请考虑以下示例:

数据

library(ggplot2)
library(ggmap)

## data simulation
set.seed(1945)

df <- tibble(x = rnorm(500, -0.7406, 0.03),
             y = rnorm(500, 51.9976, 0.03),
             z = abs(rnorm(500, 2000, 1000)))

地图,散点图,密度图

## ggmap
map <- get_map("Bletchley Park, Bletchley, Milton Keynes", zoom = 13, source = "stamen", maptype = "toner-background")
q <- ggmap(map, extent = "device", darken = .5)

## scatterplot over map
q + geom_point(aes(x, y), data = df, colour = z)

## classic density heat map
q + 
  stat_density2d(aes(x=x, y=y, fill=..level..), data=df, geom="polygon", alpha = .2) + 
  geom_density_2d(aes(x=x, y=y), data=df, colour = "white", alpha = .4) +
  scale_fill_distiller(palette = "Spectral")

如您所见,所选区域的数据相当密集,密度热图在圆边和闭合曲线上看起来很棒(除了一些最外层).

使用akima进行插值和绘图

## akima interpolation
library(akima)

df_akima <-interp2xyz(interp(x=df$x, y=df$y, z=df$z, duplicate="mean", linear = T,
                             xo=seq(min(df$x), max(df$x), length=200),
                             yo=seq(min(df$y), max(df$y), length=200)), data.frame=TRUE)

## akima plot
q +
  geom_tile(aes(x = x, y = y, fill = z), data = df_akima, alpha = .4) +
  stat_contour(aes(x = x, y = y, z = z, fill = ..level..), data = df_akima, geom = 'polygon', alpha = .4) +
  geom_contour(aes(x = x, y = y, z = z), data = df_akima, colour = 'white', alpha = .4) +
  scale_fill_distiller(palette = "Spectral", na.value = NA)

这会产生一个密集的内插值网格(以确保足够的分辨率),虽然下面的瓷砖图是可以接受的,但是轮廓图太粗糙,许多曲线没有关闭.

非线性插值使用linear = F更平滑,但显然牺牲了分辨率并且与数字(负值z)相关.

使用mgcv进行插值和绘图

## mgcv interpolation
library(mgcv)

gam <- gam(z ~ s(x, y, bs = 'sos'), data = df)
df_mgcv <- data.frame(expand.grid(x = seq(min(df$x), max(df$x), length=200),
                                  y = seq(min(df$y), max(df$y), length=200)))
resp <- predict(gam, df_mgcv, type = "response")
df_mgcv$z <- resp

## mgcv plot
q +
  geom_tile(aes(x = x, y = y, fill = z), data = df_mgcv, alpha = .4) +
  stat_contour(aes(x = x, y = y, z = z, fill = ..level..), data = df_mgcv, geom = 'polygon', alpha = .4) +
  geom_contour(aes(x = x, y = y, z = z), data = df_mgcv, colour = 'white', alpha = .4) +
  scale_fill_distiller(palette = "Spectral", na.value = NA)

使用相同的过程可以mgcv得到漂亮而平滑的图,但分辨率要低得多,实际上所有曲线都没有关闭.

问题

1. 你能否建议一个更好的方法或修改我的尝试获得类似于第一个的图(清晰,连接和高分辨率的平滑线)？

2. 是否可以关闭曲线,例如在最后一个图中(阴影区域应该计算超出图像边界)？

感谢您的时间!

Answer 1

您的地图的问题不在于您使用的插值方法，而是 ggplot 显示密度线的方式。这是这个问题的答案：Remove gaps in a stat_desired2d ggplot Chart without altering XY Limits。

密度线超出了地图，因此任何超出绘图区域的多边形都会被不适当地渲染（ggplot 将使用相应级别的下一个点来闭合多边形）。这在您的第一张地图上不会显示太多，因为插值分辨率较低。

Andrew提出的技巧是首先扩大绘图区域，以便正确渲染密度线，然后切断显示区域以隐藏多余的空间。由于我用您的第一个示例测试了他的解决方案，因此代码如下：

q + 
  stat_density2d(
    aes(x = x, y = y, fill = ..level..),
    data = df,
    geom = "polygon",
    alpha = .2,
    color = "white",
    bins = 20
  ) + 
  scale_fill_distiller(
    palette = "Spectral"
  ) +
  xlim(
    min(df$x) - 10^-5,
    max(df$x) + 10^-5
  ) +
  ylim(
    min(df$y) - 10^-3,
    max(df$y) + 10^-3
  ) +
  coord_equal(
    expand = FALSE,
    xlim = c(-.778, -.688),
    ylim = c(51.965, 52.03)
  )

唯一的区别是我使用min()- / max() +而不是固定数字并coord_equal确保地图不变形。此外，我手动指定了更多的级别（使用bin），因为通过增加绘图面积，stat_密度会自动选择较低的分辨率。

至于最佳插值方法，这取决于您的目标和您拥有的数据类型。问题不是最适合您的地图的方法是什么，而是最适合您的数据的方法是什么。这是一个非常广泛的问题，超出了本领域的范围。但这是一个很好的指南：http ://www.rspatial.org/analysis/rst/4-interpolation.html

有关如何使用 ggplot 在 R 中制作良好地图的一般想法：http://spatial.ly/r/