【问题标题】：Smoothing out ggplot2 map平滑 ggplot2 地图
【发布时间】：2015-06-08 09:05:05
【问题描述】：

以前的帖子

Cleaning up a map using geom_tile

Get boundaries to come through on states

问题/疑问

我正在尝试平滑一些数据以使用 ggplot2 进行映射。感谢@MrFlick 和@hrbrmstr，我取得了很大进展，但是在我需要列出的状态上获得“渐变”效果时遇到了问题。

这是一个示例，可以让您了解我在寻找什么：

**** 这正是我想要实现的目标。

http://nrelscience.org/2013/05/30/this-is-how-i-did-it-mapping-in-r-with-ggplot2/

(1) 如何充分利用 ggplot2 的数据？

(2) 有没有更好的方法来实现渐变效果？

目标

我想从这个赏金中实现的目标是：

(1) 对数据进行插值构建栅格对象，然后用ggplot2绘图

（或者，如果可以对当前绘图进行更多操作并且光栅对象不是一个好的策略）

(2) 使用 ggplot2 构建更好的地图

当前结果

我一直在玩很多这些不同的情节，但仍然对结果不满意，原因有两个：（1）渐变没有我想要的那么多； (2) 演示文稿可以改进，虽然我不知道该怎么做。

正如@hrbrmstr 所指出的，如果我对数据进行一些插值以产生更多数据，然后将它们放入栅格对象并使用 ggplot2 绘图，它可能会提供更好的结果。我认为这是我现在应该追求的，但鉴于我拥有的数据，我不确定如何做到这一点。

我在下面列出了到目前为止我已经完成的代码和结果。我非常感谢在这件事上的任何帮助。谢谢。

数据集

这里有两个数据集：

(1) 完整数据集 (175 mb)：PRISM_1895_db_all.csv（不可用）

https://www.dropbox.com/s/uglvwufcr6e9oo6/PRISM_1895_db_all.csv?dl=0

(2) 部分数据集 (14 mb)：PRISM_1895_db.csv（不可用）

https://www.dropbox.com/s/0evuvrlm49ab9up/PRISM_1895_db.csv?dl=0

*** 编辑：对于那些感兴趣的人，数据集不可用，但我在我的网站上发表了一篇文章，将这段代码与加利福尼亚数据的子集联系起来 http://johnwoodill.com/pages/r-code.html

情节 1

PRISM_1895_db <- read.csv("/.../PRISM_1895_db.csv")

regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
  geom_point(data = PRISM_1895_db, aes(x = longitude, y = latitude, color = APPT), alpha = .5, size = 5) +
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA) +
  coord_equal()

情节 2

PRISM_1895_db

regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
    geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
    geom_point(data = PRISM_1895_db, aes(x = longitude, y = latitude, color = APPT), alpha = .5, size = 5, shape = 15) +
    geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA) +
    coord_equal()

情节 3

   PRISM_1895_db <- read.csv("/.../PRISM_1895_db.csv")

    regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas","minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")

ggplot() + 
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group)) +
  stat_summary2d(data=PRISM_1895_db, aes(x = longitude, y = latitude, z = APPT)) +
  geom_polygon(data=subset(map_data("state"), region %in% regions), aes(x=long, y=lat, group=group), color="white", fill=NA)

【问题讨论】：

感谢您发布这个精彩的问题。我也对它很感兴趣，但数据的链接不可用。你能提供一个新的链接吗？非常感谢。
johnwoodill.github.io/code.html
@YangYang 抱歉回复晚了，感谢您发布指向我网站的链接。请参阅链接以获取使用加利福尼亚的示例，而不是此处描述的示例。

标签： r ggplot2 geospatial interpolation

【解决方案1】：

CRAN spatial view 让我开始学习“克里金法”。下面的代码在我的笔记本电脑上运行大约需要 7 分钟。您可以尝试更简单的插值（例如，某种样条）。您还可以从高密度区域中删除一些位置。您不需要所有这些点来获得相同的热图。据我所知，使用ggplot2 创建真正的渐变并没有简单的方法（gridSVG 有几个选项，但没有像您在精美的 SVG 编辑器中找到的“网格渐变”那样）。

根据要求，这里是使用样条线的插值（快得多）。很多代码取自Plotting contours on an irregular grid。

克里金法代码：

library(data.table)
library(ggplot2)
library(automap)

# Data munging
states=c("AR","IL","MO")
regions=c("arkansas","illinois","missouri")
PRISM_1895_db = as.data.frame(fread("./Downloads/PRISM_1895_db.csv"))
sub_data = PRISM_1895_db[PRISM_1895_db$state %in% states,c("latitude","longitude","APPT")]
coord_vars = c("latitude","longitude")
data_vars = setdiff(colnames(sub_data), coord_vars)
sp_points = SpatialPoints(sub_data[,coord_vars])
sp_df = SpatialPointsDataFrame(sp_points, sub_data[,data_vars,drop=FALSE])

# Create a fine grid
pixels_per_side = 200
bottom.left = apply(sp_points@coords,2,min)
top.right = apply(sp_points@coords,2,max)
margin = abs((top.right-bottom.left))/10
bottom.left = bottom.left-margin
top.right = top.right+margin
pixel.size = abs(top.right-bottom.left)/pixels_per_side
g = GridTopology(cellcentre.offset=bottom.left,
             cellsize=pixel.size,
             cells.dim=c(pixels_per_side,pixels_per_side))

# Clip the grid to the state regions
map_base_data = subset(map_data("state"), region %in% regions)
colnames(map_base_data)[match(c("long","lat"),colnames(map_base_data))] = c("longitude","latitude")
foo = function(x) {
  state = unique(x$region)
  print(state)
  Polygons(list(Polygon(x[,c("latitude","longitude")])),ID=state)
}
state_pg = SpatialPolygons(dlply(map_base_data, .(region), foo))
grid_points = SpatialPoints(g)
in_points = !is.na(over(grid_points,state_pg))
fit_points = SpatialPoints(as.data.frame(grid_points)[in_points,])

# Do kriging
krig = autoKrige(APPT~1, sp_df, new_data=fit_points)
interp_data = as.data.frame(krig$krige_output)
colnames(interp_data) = c("latitude","longitude","APPT_pred","APPT_var","APPT_stdev")

# Set up map plot
map_base_aesthetics = aes(x=longitude, y=latitude, group=group)
map_base = geom_polygon(data=map_base_data, map_base_aesthetics)
borders = geom_polygon(data=map_base_data, map_base_aesthetics, color="black", fill=NA)

nbin=20
ggplot(data=interp_data, aes(x=longitude, y=latitude)) + 
  geom_tile(aes(fill=APPT_pred),color=NA) +
  stat_contour(aes(z=APPT_pred), bins=nbin, color="#999999") +
  scale_fill_gradient2(low="blue",mid="white",high="red", midpoint=mean(interp_data$APPT_pred)) +
  borders +
  coord_equal() +
  geom_point(data=sub_data,color="black",size=0.3)

样条插值代码：

library(data.table)
library(ggplot2)
library(automap)
library(plyr)
library(akima)

# Data munging
sub_data = as.data.frame(fread("./Downloads/PRISM_1895_db_all.csv"))
coord_vars = c("latitude","longitude")
data_vars = setdiff(colnames(sub_data), coord_vars)
sp_points = SpatialPoints(sub_data[,coord_vars])
sp_df = SpatialPointsDataFrame(sp_points, sub_data[,data_vars,drop=FALSE])

# Clip the grid to the state regions
regions<- c("north dakota","south dakota","nebraska","kansas","oklahoma","texas",
            "minnesota","iowa","missouri","arkansas", "illinois", "indiana", "wisconsin")
map_base_data = subset(map_data("state"), region %in% regions)
colnames(map_base_data)[match(c("long","lat"),colnames(map_base_data))] = c("longitude","latitude")
foo = function(x) {
  state = unique(x$region)
  print(state)
  Polygons(list(Polygon(x[,c("latitude","longitude")])),ID=state)
}
state_pg = SpatialPolygons(dlply(map_base_data, .(region), foo))

# Set up map plot
map_base_aesthetics = aes(x=longitude, y=latitude, group=group)
map_base = geom_polygon(data=map_base_data, map_base_aesthetics)
borders = geom_polygon(data=map_base_data, map_base_aesthetics, color="black", fill=NA)

# Do spline interpolation with the akima package
fld = with(sub_data, interp(x = longitude, y = latitude, z = APPT, duplicate="median",
                            xo=seq(min(map_base_data$longitude), max(map_base_data$longitude), length = 100),
                            yo=seq(min(map_base_data$latitude), max(map_base_data$latitude), length = 100),
                            extrap=TRUE, linear=FALSE))
melt_x = rep(fld$x, times=length(fld$y))
melt_y = rep(fld$y, each=length(fld$x))
melt_z = as.vector(fld$z)
level_data = data.frame(longitude=melt_x, latitude=melt_y, APPT=melt_z)
interp_data = na.omit(level_data)
grid_points = SpatialPoints(interp_data[,2:1])
in_points = !is.na(over(grid_points,state_pg))
inside_points = interp_data[in_points, ]

ggplot(data=inside_points, aes(x=longitude, y=latitude)) + 
  geom_tile(aes(fill=APPT)) + 
  stat_contour(aes(z=APPT)) +
  coord_equal() + 
  scale_fill_gradient2(low="blue",mid="white",high="red", midpoint=mean(inside_points$APPT)) +
  borders

【讨论】：

这是一个很好的答案，正是我正在寻找的。但是，由于我需要一个包含更多州的更大区域，这需要一段时间才能运行，而且我需要绘制 75 年。你能帮忙应用一个样条曲线，让它运行得更快一点吗？
非常感谢您。很好的答案！你能解释一下xo=seq(1.02*min(longitude), max(longitude), length = 400), yo=seq(0.96*min(latitude), max(latitude), length = 400),我不完全确定为什么1.02和0.96 ....谢谢
那是一个丑陋的黑客。混乱我知道。我希望用于插值的网格包含所有状态。问题是min(latitude) 就在德克萨斯州南端的北部，所以我将它乘以 0.96 以将网格向南延伸一点。我现在意识到最好做xo=seq(min(map_base_data$longitude), max(map_base_data$longitude), length = 100)。我已经更改了答案中的代码以反映这一点。
好的，很好用。您将长度缩小到 100 是否有原因？我唯一不确定的是rep(x, times) 和rep(x, each)。我知道它会复制 fld$x 或 fld$y 中的数据，但为什么要使用 times 和 each？我终于完成了对代码的审查，这很棒。只是不确定最后一部分。再次感谢。
很难解释rep(x,times) 和rep(x,each) 在没有看到的情况下做了什么。阅读帮助文件；制作一个小矩阵；并玩弄它。这就是我所做的。我将长度缩小到因为这样会更快，尽管结果图像的分辨率会更低（100x100 而不是 400x400 平铺）。您可以根据自己的喜好调整网格大小。

【解决方案2】：

之前的答案可能不是您需要的最佳（或准确）。这有点小技巧：

gg <- ggplot() 
gg <- gg + geom_polygon(data=subset(map_data("state"), region %in% regions), 
                        aes(x=long, y=lat, group=group))
gg <- gg + geom_point(data=PRISM_1895_db, aes(x=longitude, y=latitude, color=APPT), 
                      size=5, alpha=1/15, shape=19)
gg <- gg + scale_color_gradient(low="#023858", high="#ece7f2")
gg <- gg + geom_polygon(data=subset(map_data("state"), region %in% regions), 
                        aes(x=long, y=lat, group=group), color="white", fill=NA)
gg <- gg + coord_equal()
gg

这需要将geom_point 中的size 更改为更大的绘图，但您会获得比stat_summary2d 行为更好的渐变效果，并且它传达的信息相同。

另一种选择是在您拥有的经度和纬度之间插入更多 APPT 值，然后将其转换为更密集的栅格对象并使用 geom_raster 绘制它，就像您提供的示例一样。

【讨论】：

hmmmm....这有点难以看出差异....是否有一个参数可以帮助清理这个？感谢您的回答
感谢您的回答。这似乎比我迄今为止所做的要好一些，所以我会玩这个。我认为您对插入一些附加值是正确的，但是我将其构建为一个函数，因此我可以运行多个大约 80 年的图，因此从这方面来看可能有点困难。这是初步研究，我只需要为我的顾问和演示文稿提供一些漂亮的情节来传达这个想法。最终，我将对这些数据进行建模以适应三角函数以获得每日平均值，然后您的栅格想法将是完美的。谢谢！
您认为您可以在这方面进一步帮助我吗？我没有得到任何关于赏金的进一步关注，并且非常感谢您的帮助。特别是，我认为插值以增加图表上的点数是正确的。但我不确定如何将它们应用到栅格对象，甚至不知道如何创建附加点。