平滑 ggplot2 地图
Smoothing out ggplot2 map
以前的帖子
Problem/Question
我正在尝试整理一些数据以使用 ggplot2 进行映射。感谢@MrFlick 和@hrbrmstr，我取得了很大进步，但在我需要列出的状态上获得 "gradient" 效果时遇到了问题。
这里有一个例子，可以让您了解我在寻找什么：
**** 这正是我要实现的目标。
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
*** 编辑：对于那些感兴趣的人，数据集不可用，但我在我的网站上制作了一个 post，将此代码与位于 [=18= 的加利福尼亚数据的子集连接起来]
情节 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 <- 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, 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)
之前的回答可能不是您需要的最佳（或准确）答案。这有点骇人听闻：
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 绘制它，就像您提供的示例中那样.
CRAN spatial view 让我开始 "Kriging"。下面的代码在我的笔记本电脑上需要大约 7 分钟才能 运行 。您可以尝试更简单的插值（例如，某种样条曲线）。您还可以从 high-density 区域中删除一些位置。您不需要所有这些点来获得相同的热图。据我所知，没有简单的方法可以使用 ggplot2 创建真正的渐变（gridSVG 有一些选项，但与您在精美的 SVG 编辑器中找到的 "grid gradient" 完全不同） .
根据要求，这里是使用样条插值（快得多）。很多代码取自 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
以前的帖子
Problem/Question
我正在尝试整理一些数据以使用 ggplot2 进行映射。感谢@MrFlick 和@hrbrmstr，我取得了很大进步，但在我需要列出的状态上获得 "gradient" 效果时遇到了问题。
这里有一个例子，可以让您了解我在寻找什么：
**** 这正是我要实现的目标。
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
*** 编辑：对于那些感兴趣的人，数据集不可用，但我在我的网站上制作了一个 post，将此代码与位于 [=18= 的加利福尼亚数据的子集连接起来]
情节 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 <- 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, 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)
之前的回答可能不是您需要的最佳（或准确）答案。这有点骇人听闻：
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 绘制它，就像您提供的示例中那样.
CRAN spatial view 让我开始 "Kriging"。下面的代码在我的笔记本电脑上需要大约 7 分钟才能 运行 。您可以尝试更简单的插值（例如，某种样条曲线）。您还可以从 high-density 区域中删除一些位置。您不需要所有这些点来获得相同的热图。据我所知，没有简单的方法可以使用 ggplot2 创建真正的渐变（gridSVG 有一些选项，但与您在精美的 SVG 编辑器中找到的 "grid gradient" 完全不同） .
根据要求，这里是使用样条插值（快得多）。很多代码取自 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