平滑ggplot2地图

问题描述

以前的帖子

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

按照要求，这里是使用样条插值（更快）的插值。许多代码来自在不规则网格上绘制轮廓。

克里格代码：

  library（data.table）
库（ggplot2）
库（automap）
 
＃数据转移
 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％ （'纬度'，'经度'，'APPT'）] 
 coord_vars = c（纬度，经度）
 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]）
 
＃创建一个细网格
 pixels_per_side = 200 
 bottom.left = apply（sp_分数@ 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））
 
＃将网格剪切到状态区域
 map_base_data =子集（map_data（state），region％％in％regions）
 colnames（map_base_data）[match（c（long，lat），colnames（map_base_data））] = c （x）{
 state = unique（x $ region）
 print（state）
 Polygons（list（Polygon（x，longitude），latitude）
 foo = [，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_p oint，]）
 
＃做克里格
 krig = autoKrige（APPT〜1，sp_df，new_data = fit_points）
 interp_data = as.data.frame（krig $ krige_output）
 colnames（interp_data）= c（纬度，经度，APPT_pred，APPT_var，APPT_stdev）
 
＃设置地图图
 map_base_aesthetics = aes （x =经度，y =纬度，组=组）
 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））+ 
 border + 
 coord_equal（）+ 
 geom_point（data = sub_data，color =black，size = 0.3）
  

样条插值的代码：

 库（data.ta （bb）b 
 
 
 
 
 
 
 
 
 $ b $ = 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]）
 
＃将网格剪辑到状态（北达科他州，南达科他州，内布拉斯加州，堪萨斯州，俄克拉荷马州，德克萨斯州，美元b，明尼苏达州，爱荷华州，密苏里州 ，阿肯色州，伊利诺伊州，印第安纳州，威斯康星州）
 map_base_data =子集（map_data（state），region％％in regions 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（列表（Polygon（x [，c（latitude，longitude）]）），ID = state）
} 
 state_pg = SpatialPolygons（dlply（map_base_data，。（region），foo））
 
＃设置地图绘制
 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）
 
＃使用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 $经度），长度= 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）
g rid_points = SpatialPoints（interp_data [，2：1]）$ ​​b $ b 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））+ 
 border 
  

Previous Posts

Cleaning up a map using geom_tile

Get boundaries to come through on states

Problem/Question

I'm trying to smooth out some data to map with ggplot2. Thanks to @MrFlick and @hrbrmstr, I've made a lot of progress, but am having problems getting a "gradient" effect over the states I need listed.

Here is an example to give you an idea about what I'm looking for :

**** This is exactly what I'm trying to achieve.

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

(1) How can I make the most of ggplot2 with my data?

(2) Is there a better method for achieving a gradient effect?

Goals

Goals I would like to achieve from this bounty is :

(1) Interpolate the data to build a raster object and then plot with ggplot2

(or, if more can be done with the current plot and the raster object is not a good strategy)

(2) Build a better map with ggplot2

Current Results

I've been playing around with a lot of these different plots, but am still not happy with the results for two reasons: (1) The gradient isn't telling as much as I'd like; and (2) The presentation could be improved, although I'm not sure how to do it.

As @hrbrmstr has pointed out, it might provide better results if I did some interpolating with the data to produce more data, and then fit those into a raster object and plot with ggplot2. I think this is what I should be after at this point, but I'm not sure how to do this given the data I have.

I've listed below the code I have done so far with the results. I really appreciate any help in this matter. Thank you.

Data sets

Here are two data sets :

(1) Full data set (175 mb) : PRISM_1895_db_all.csv (NOT AVAILABLE)

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

(2) Partial data set (14 mb) : PRISM_1895_db.csv (NOT AVAILABLE)

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

*** EDIT : To those interested, the data sets are not available, but I've made a post on my website that connects this code with a subset of California data at http://johnwoodill.com/pages/r-code.html

Plot 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()

Plot 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()

Plot 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)

解决方案

The CRAN spatial view got me started on "Kriging". The code below takes ~7 minutes to run on my laptop. You could try simpler interpolations (e.g., some sort of spline). You might also remove some of the locations from the high-density regions. You don't need all of those spots to get the same heatmap. As far as I know, there is no easy way to create a true gradient with ggplot2 (gridSVG has a few options but nothing like the "grid gradient" you would find in a fancy SVG editor).

As requested, here is interpolation using splines (much faster). Alot of the code is taken from Plotting contours on an irregular grid.

Code for kriging:

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)

Code for spline interpolation:

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

这篇关于平滑ggplot2地图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！