地理热力/等高线地图的最佳空间插值方法是什么？

Question

地理热力/等高线地图的最佳空间插值方法是什么？

13

我想使用类似于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进行相同的过程会得到一个漂亮而平滑的图形，但分辨率要低得多，并且实际上所有曲线都不闭合。

问题

请问您能否提供更好的方法或修改我的尝试，以获得类似第一个图的绘图（干净、连贯、平滑的线条，高分辨率）？
是否有可能关闭曲线，例如在最后一张图中（阴影区域应计算超出图像边界）？

感谢您的时间！

- Harold Cavendish

1

如果您在https://gis.stackexchange.com/上发布帖子，可能会得到更好的回应。 - Tung

1

@Tung 谢谢您的建议！我会在这里再多等一些时间，如果没有找到解决方案，我会在那里发布。 - Harold Cavendish

考虑设置悬赏以吸引更多关注您的问题。 - Tung

2个回答

2

抱歉，我现在无法运行您的示例以提供详细信息。但是请尝试使用 automap 包中的 autoKrige() 函数。

Kriging 是一种很好的插值方法，只需确保您的数据符合要求即可。这里有一个很好的指南： https://gisgeography.com/kriging-interpolation-prediction/

- Carlos Eduardo Lagosta

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Carlos Eduardo Lagosta · Accepted Answer

你的地图问题不在于你使用的插值方法，而是ggplot显示密度线的方式。这里有一个解决方案：去除stat_density2d ggplot图表中的空白，而不修改XY限制。

密度线超出了地图范围，所以任何超出绘图区域的多边形都会被不当地渲染（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_density 会自动选择较低的分辨率。

至于最佳插值方法，这取决于您的目标和数据类型。问题不是什么是您的地图的最佳方法，而是什么是您的数据的最佳方法。这是一个非常广泛的问题，超出了本空间的范围。但这里有一个好的指南：http://www.rspatial.org/analysis/rst/4-interpolation.html 关于如何在 R 中使用 ggplot 制作出好的地图的一般思路：http://spatial.ly/r/