# 13.其他图形
# 13.1 绘制相关矩阵图
mtcars
mcor <- cor(mtcars)
# 输出mcor,保留两位小数
round(mcor,digits = 2)
# 如果数据含有不能用来计算系统的任何列，应该先将这些列剔除。
# 如果在原始数据中存在缺失值（NA)，得到的相关矩阵中也会有缺失值。
# 可以使用函数选项use="complete.obs"或者use="pairwise.complete.obs"
# 使用corrplot包来绘制相关矩阵图
library(corrplot)
corrplot(mcor)
# 使用颜色方块和黑色文本标签，并且上边的文本标签呈45°右倾
corrplot(mcor,method = "shade",shade.col = NA,tl.col = "black",tl.srt = 45)
# 对矩阵重新排序，这样相近的变量在图中会更近，其中使用的参数时order="AOE"（前两个特征向量的角排序）
# 生成一个淡一点的调色板
col <- colorRampPalette(c("#BB4444","#EE9988","#FFFFFF","#77AADD","#4477AA"))
corrplot(mcor,method = "shade",shade.col = NA,tl.col = "black",tl.srt = 45,col = col(200),addCoef.col = "black",cl.pos="no",order = "AOE")
# 13.2绘制函数曲线
# 使用stat_function()函数
library(ggplot2)
# 这个数据库仅仅用来设定范围
p <- ggplot(data.frame(x=c(-3,3)),aes(x=x))
p+stat_function(fun = dnorm)
# 某些函数需要格外的参数
p+stat_function(fun = dt,args = list(df=2))
# 绘制自定义的函数
myfun <- function(xvar){
1/(1+exp(-xvar+10))
}
ggplot(data.frame(x=c(0,20)),aes(x=x))+stat_function(fun = myfun)
# 计算函数值时默认使用了给出x范围内的101个点。
# 为了让曲线更加光滑，可以给stat_function()传递一个更大的n
ggplot(data.frame(x=c(0,20)),aes(x=x))+stat_function(fun = myfun,n=200)
# 13.3 在函数曲线下添加阴影
# 根据曲线函数定义一个新的函数，把x范围外对应的值替换为NA
# 在0<x<2时返回dnorm(x)，其他的时候返回NA
dnorm_limit <- function(x){
y <- dnorm(x)
y[x<0|x>2] <- NA
return(y)
}
# ggplot()使用“哑”数据
p <- ggplot(data.frame(x=c(-3,3)),aes(x=x))
p+stat_function(fun = dnorm_limit,geom = "area",fill="blue",alpha=0.2)+stat_function(fun = dnorm)
# 给这个函数传入的是一个向量，并不是一个单独的值。
# R中有第一类函数，可以写一个函数来返回一个闭包。也就是说，可以编写一个能够编写函数的函数
limitRange <- function(fun,min,max){
function(x){
y <- fun(x)
y[x<min|y>max] <-NA
return(y)
}
}
# 调用这个函数来生成另一个函数
dlimit <- limitRange(dnorm,0,2)
# 仅对0-2之间的输入返回输出值
dlimit(-2:4)
# 使用limitRange()来生成函数，并传递给stat_function():
p+stat_function(fun = dnorm)+stat_function(fun = limitRange(dnorm,0,2),geom = "area",fill="blue",alpha=0.2)
# limitRange()函数可以用来生成任何函数的“区间限制式”函数。
# 13.4绘制网络图
# 使用igraph包。
# 首先给graph()函数传递一个包含所有边的向量，然后绘制结果对象
install.packages("igraph")
library(igraph)
# 指定一个有向图的边
gd <- graph(c(1,2,2,3,2,4,1,4,5,5,3,6))
plot(gd)
# 一个无向图
gu <- graph(c(1,2,2,3,2,4,1,4,5,5,2,6),directed = FALSE)
# 不画标签
plot(gu,vertex.label=NA)
# 图对象的结构
str(gd)
str(gu)
# 在网络图中，节点的位置并不是由所给数据确定的，它们是随机放置的。
set.seed(200)
plot(gu)
# 从数据库直接生成图
# Fruchterman-Reingold布局算法。该算法的主要思想是所有节点之间都有电磁斥力，但是连接节点的边像弹簧一样，会把相应的节点拉在一起
library(gcookbook)
madmen2
# 从数据集中生成图对象
g <- graph.data.frame(madmen2,directed = TRUE)
# 移除多余的空白边
par(mar=c(0,0,0,0))
plot(g,layout=layout.fruchterman.reingold,vertex.size=8,edge.arrow.size=0.5,vertex.label=NA)
# 使用圆圈布局
g <- graph.data.frame(madmen,directed = FALSE)
par(mar=c(0,0,0,0))
plot(g,layout=layout.circle,vertex.size=8,vertex.label=NA)
# igraph的一个替代选择是Rgraphviz，该包是Graphviz(一个绘制网络图的开源库)的前端。
# 它的优点是能更方便地处理标签，而且也更容易控制布局。但是它比较难安装，其由Bioconductor系统库维护
# 13.5 在网络图中使用文本标签
# 可以给vertex.label参数传递一个命名向量
library(igraph)
library(gcookbook)
# 赋值madmen并删除偶数行
m <- madmen[1:nrow(madmen)%%2==1,]
g <- graph.data.frame(m,directed = FALSE)
# 输出节点名称
V(g)$name
plot(g,layout=layout.fruchterman.reingold,
vertex.size = 4, #让节点更小
vertex.label = V(g)$name, #设置标签
vertex.label.cex = 0.8, #小号字体
vertex.label.dist = 0.4, #标签和节点的位置错开
vertex.label.color = "black"
)
# 第二种方法:使用V()$xxx <-来代替vertex.xxx参数传递值。
V(g)$size <- 4
V(g)$label <- V(g)$name
V(g)$label.cex <- 0.8
V(g)$label.dise <- 0.4
V(g)$label.color <- "black"
# 设置这个图的属性
g$layout <- layout.fruchterman.reingold
plot(g)
# 设置边的属性，使用E()函数或者给edge.xxx参数传递相应的值
# 查看边
E(g)
# 将几个边的名字赋值为“M”
E(g)[c(2,11,19)]$label <- "M"
# 将所有边颜色设置为灰色，然后把其中的几个变为红色
E(g)$color <- "grey70"
E(g)[c(2,11,19)]$color <- "red"
plot(g)
# 13.6 如何绘制热图
# 使用geom_tile()或者geom_raster()，并将一个连续变量映射到fill上
presidents
str(presidents)
# 将其转化为ggplot()可用的数据框格式，其中的列都是数值形式
pres_rating <- data.frame(
rating = as.numeric(presidents),
year = as.numeric(floor(time(presidents))),
quarter = as.numeric(cycle(presidents))
)
pres_rating
# 基础图形
p <- ggplot(pres_rating,aes(x=year,y=quarter,fill=rating))
# 使用geom_tile()
p+geom_tile()
# 使用geom_raster()
p+geom_raster()
# geom_tile()和geom_raster()的结果看起来一样，实际上，它们是有区别的。
# 自定义热图的外观
p+geom_tile()+scale_x_continuous(breaks = seq(1940,1976,by=4))+scale_y_reverse()+scale_fill_gradient2(midpoint = 50,mid = "grey70",limits=c(0,100))
# 13.7 绘制三维散点图
# 使用rgl包，该包提供了OpenGL图形库的3D绘图接口。要画三维散点图，可使用plot3d()函数。
# 其输入参数可以是两种形式：（1）一个数据框，前三列分别表示x,y,z的坐标；（2）直接传递三个向量，分别表示x,y,z的坐标。
install.packages("rgl")
library(rgl)
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,type = "s",size = 0.75,lit=FALSE)
# 默认情况下，plot3d()函数使用立方体式的点。
# type="s"选择球形点，size=0.75让点变小，lit=FASLE关闭3D灯光（否则点会闪闪发亮）
# 添加数值线段来增强空间点的位置表达力度
# 交错出现两个向量的值。
interleave <- function(v1,v2) as.vector(rbind(v1,v2))
# 绘制点
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,xlab = "Weight",ylab="Displacement",zlab="MPG",size = .75,type = "s",lit=FALSE)
# 添加线段
segments3d(interleave(mtcars$wt,mtcars$wt),interleave(mtcars$disp,mtcars$disp),interleave(mtcars$mpg,min(mtcars$mpg)),alpha=0.4,col="blue")
# 不画坐标刻度和标签
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,xlab = "",ylab = "",zlab = "",axes=FALSE,size=.75,type = "s",lit=FALSE)
segments3d(interleave(mtcars$wt,mtcars$wt),interleave(mtcars$disp,mtcars$disp),interleave(mtcars$mpg,min(mtcars$mpg)),alpha=0.4,col="blue")
# 绘制盒子
rgl.bbox(color="grey50", # 表面颜色为grey60，黑色文本
emission="grey50", # 光照颜色为grey50
xlen = 0,ylen = 0,zlen = 0 # 不添加刻度
)
# 设置默认颜色为黑
rgl.material(color = "black")
# 在指定边添加坐标轴标签，可能的值类似于“x--”，“x-+”,"x+-"，“x++”
axes3d(edges = c("x--","y+-","z--"),nticks = 6, # 每个轴上6个刻度线
cex=.75#较小的字体
)
# 添加坐标标签，“line"指定标签和坐标轴的距离
mtext3d("Weight",edge = "x--",line = 2)
mtext3d("Displacement",edge = "y+-",line = 3)
mtext3d("MPG",edge = "z--",line = 3)
# 13.8在三维图上添加预测曲面
# 首先，需要定义一些功能函数来得到模型的预测值
# 给定一个模型，根据xvar和yvar预测zvar
# 默认为变量x和y的范围，生成16*16的网格
predictgrid <- function(model,xvar,yvar,zvar,res=16,type=NULL){
# 计算预测变量的范围，下面的代码对lm、glm以及其他模型方法都适用，
# 但针对其他模型方法时可能需要适当调整
xrange <- range(model$model[[xvar]])
yrange <- range(model$model[[yvar]])
newdata <- expand.grid(x=seq(xrange[1],xrange[2],length.out = res),
y=seq(yrange[1],yrange[2],length.out = res))
names(newdata) <- c(xvar,yvar)
newdata[[zvar]] <- predict(model,newdata = newdata,type = type)
newdata
}
# 将长数据框中的x,y,z转化为列表
# 其中x,y为行列值，z为矩阵
df2mat <- function(p,xvar=NULL,yvar=NULL,zvar=NULL){
if(is.null(xvar)) xvar <- names(p)[1]
if(is.null(yvar)) yvar <- names(p)[2]
if(is.null(zvar)) zvar <- names(p)[3]
x <- unique(p[[xvar]])
y <- unique(p[[yvar]])
z <- unique(p[[zvar]],nrow=length(y),ncol=length(x))
m <- list(x,y,z)
names(m) <- c(xvar,yvar,zvar)
m
}
# 交错出现两个向量的元素
interleave <- function(v1,v2) as.vector(rbind(v1,v2))
# 利用surface3d()函数在原散点图上添加网格式的预测图
library(rgl)
m <- mtcars
# 生成线性模型
mod <- lm(mpg~wt+disp+wt:disp,data=m)
# 根据wt和disp，得到mpg的预测值
m$pred_mpg <- predict(mod)
# 根据wt和disp的网格，得到mpg的预测值
mpgrid_df <- predictgrid(mod,"wt","disp","mpg")
mpgrid_list <- df2mat(mpgrid_df)
# 根据数据点绘图
plot3d(m$wt,m$disp,m$mpg,type = "s",size = 0.5,lit=FALSE)
# 添加预测点（较小）
spheres3d(m$wt,m$disp,m$pred_mpg,alpha=0.4,type = "s",size = 0.5,lit=FALSE)
# 添加表示误差的线段
segments3d(interleave(m$wt,m$wt),interleave(m$disp,m$disp),interleave(m$mpg,m$pred_mpg),alpha=0.4,col="red")
# 添加预测点网络
surface3d(mpgrid_list$wt,mpgrid_list$disp,mpgrid_list$mpg,alpha=0.4,front = "lines",back="lines")
# 调节图形的外观，逐一添加图形的各个组件
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,xlab = "",ylab = "",zlab = "",axes = FALSE,size = .5,type = "s",lit=FALSE)
# 添加预测点（较小）
spheres3d(m$wt,m$disp,m$pred_mpg,alpha=0.4,type="s",size=0.5,lit=FALSE)
# 添加误差线段
segments3d(interleave(m$wt,m$wt),interleave(m$disp,m$disp),interleave(m$mpg,m$pred_mpg),alpha=0.4,col="red")
# 添加预测值网络
surface3d(mpgrid_list$wt,mpgrid_list$disp,mpgrid_list$mpg,alpha=0.4,front="lines",back="lines")
# 绘制盒子
rgl.bbox(color="grey50", # 表面颜色为grey60,黑色文本
emission = "grey50", # 光照颜色为grey50
xlen = 0,ylen = 0,zlen = 0
)
# 对象默认色设置为黑色
rgl.material(color = "black")
# 在指定边添加坐标轴标签。可能的值类似于“x--”,"x-+","x+-"和“x++"
axes3d(edges = c("x--","y+-","z--"),nticks = 6, # 每个轴上6个刻度线
cex=.75 # 较小字体
)
# 添加坐标标签，“line"指定标签和坐标轴的距离。
mtext3d("Weight",edge = "x--",line=2)
mtext3d("Displacement",edge = "y+-",line=3)
mtext3d("MPG",edge = "z--",line = 3)
# 13.9 保存三维图
# 使用rgl.snapshot()来保存rgl包绘制的位图。它会精确捕捉屏幕上的图形
library(rgl)
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,type = "s",size = 0.75,lit=FALSE)
rgl.snapshot("3dplot.png",fmt = "png")
# 也可以使用rgl.postscript()保存为PostScript或PDF格式文件：
rgl.postscript("3dplot.pdf",fmt = 'pdf')
rgl.postscript("3dplot.ps",fmt = 'ps')
# PostScript和PDF输出文件并不支持rgl以来的OpenGL库的很多特性。
# 保存当前视角
view <- par3d("userMatrix")
# 恢复保存的视角
par3d(userMatrix = view)
# 将视角保存为代码，使用dput()函数，然后将输出复制粘贴到自己的代码中
dput(view)
view <- structure(c(1, 0, 0, 0, 0, 0.342020143325668, -0.939692620785909,
0, 0, 0.939692620785909, 0.342020143325668, 0, 0, 0, 0, 1), .Dim = c(4L,
4L))
par3d(userMatrix = view)
# 13.10 三维图动画
# 旋转三维图能够更完整、多方位地观察数据，可以在play3d()中使用spin3d()来生成三维动画
library(rgl)
plot3d(mtcars$wt,mtcars$disp,mtcars$mpg,type = "s",size = 0.75,lit=FALSE)
# play3d(spin3d())
# 默认情况下，图像会绕着z轴（竖直的轴）旋转。
# 绕x轴转动，每分钟4转，持续20秒钟
play3d(spin3d(axis = c(1,0,0),rpm = 4),duration = 20)
# 使用movie3d()来保存动画，方法和play3d()一样。它将会生成一系列.png格式的图片文件，每个文件代表一帧，然后利用ImageMagick软件提供的convert命令将这些文件合并转化为.gif动画文件
# 绕z轴转动，每分钟4转，持续15秒
movie3d(spin3d(axis = c(0,0,1),rpm = 4),duration = 15,fps = 50)
# 输出的文件将会被存放在一个临时文件夹中，地址名会在R窗口中输出。
# 如果不想利用ImageMagick将输出的图片装换为.gif。可以设置convert = FALSE，然后用其他软件将这一系列的.png文件转换为动画
# 13.11绘制谱系图
# 使用hclust()并画出它的结果
library(gcookbook)
# 得到2009年的数据
c2 <- subset(countries,Year==2009)
# 去掉含有NA的行
c2 <- c2[complete.cases(c2),]
# 随机选择25个国家
# (设定随机种子保证可重复性)
set.seed(201)
c2 <- c2[sample(1:nrow(c2),25),]
c2
# 去掉聚类中不需要的列
rownames(c2) <- c2$Name
c2 <- c2[,4:7]
c2
# 数据标准化
c3 <- scale(c2)
c3
# scale()函数默认是将每一列相对于其标准差进行标准化
hc <- hclust(dist(c3))
# 画树状图
plot(hc)
# 对齐文本
plot(hc,hang = -1)
# 聚类分析是把n维空间中把点分配到类的一种简单方法。层次聚类分析则将每组分成两个更小的组，以谱系展示
# 对于距离的计算，使用的是默认的方法——“euclidean"，这个方式计算的是点与点之间的欧氏距离；还有其他方法，比如”maximum","manhattan","canberra","binary"和"minkowski".
# hclust()函数提供了几种做聚类分析的方法，默认的方法是“complete”;其他可用的方法包括“ward”,"single","average","mcquitty","median"和”centroid“。
# 13.12 绘制向量场
# 使用geom_segment()函数
library(gcookbook)
isabel
# 每段都有一个起点和一个终点。用x和y值作为每段起点，然后在此基础上分别加上vx和vy的一个比例，以此作为终点，如果不按比例缩小这些取值，线条就会变得很长
islice <- subset(isabel, z==min(z))
library(ggplot2)
ggplot(islice,aes(x=x,y=y))+geom_segment(aes(xend=x+vx/50,yend=y+vy/50),size=0.25) # 线段0.25mm粗
# 向量场有两个问题：数据分辨率太高不太容易阅读，而且每段没有箭头表示方向。为了降低分辨率，定义一个函数every_n(),在数据的每n个值中保留一个，其他的去掉。
# 选择z取值等于z的最小值的部分数据
islice <- subset(isabel,z==min(z))
# 向量x中每“by”个只里面保留一个
every_n <- function(x,by=2){
x <- sort(x)
x[seq(1,length(x),by=by)]
}
# x和y每四个值保留一个
keepx <- every_n(unique(isabel$x),by=4)
keepy <- every_n(unique(isabel$y),by=4)
# 保留那些x值在keepx中并且y值在keepy中的数据
islicesub <- subset(islice,x%in%keepx&y%in%keepy)
# 使用arrow()，需要加载grid包
library(grid)
# 用子集画图，箭头的长度为0.1cm
ggplot(islicesub,aes(x=x,y=y))+geom_segment(aes(xend=x+vx/50,yend=y+vy/50),arrow = arrow(length = unit(0.1,"cm")),size=0.25)
# 箭头的一个影响是，短的向量会表现得比它实际长度的比例更大，这会导致曲解数据，为了减轻这种影响，把速度映射到其他属性上可能是有用的
# “speed”列包含z的部分，计算水平速度
islicesub$speedxy <- sqrt(islicesub$vx^2+islicesub$vy^2)
# 映射速度到透明度alpha
ggplot(islicesub,aes(x=x,y=y))+geom_segment(aes(xend=x+vx/50,yend=y+vy/50,alpha=speed),arrow = arrow(length = unit(0.1,"cm")),size=0.6)
install.packages("maps")
library(maps)
# 得到美国地图数据
usa <- map_data("usa")
# 把数据映射到颜色上，颜色从“grey80"到”darkred"
ggplot(islicesub,aes(x=x,y=y))+geom_segment(aes(xend=x+vx/50,yend=y+vy/50,colour=speed),arrow = arrow(length = unit(0.1,"cm")),size=0.6)+
scale_color_continuous(low="grey80",high = "darkred")+geom_path(aes(x=long,y=lat,group=group),data=usa)+coord_cartesian(xlim=range(islicesub$x),ylim = range(islicesub$y))
# isabel数据集是一个三维数据，画一个分面的图形
# x和y中每5个值保留1个，z中每两个值保留一个
keepx <- every_n(unique(isabel$x),by=5)
keepy <- every_n(unique(isabel$y),by=5)
keepz <- every_n(unique(isabel$z),by=2)
isub <- subset(isabel,x%in%keepx&y%in%keepy&z%in%keepz)
ggplot(isub,aes(x=x,y=y))+geom_segment(aes(xend=x+vx/50,yend=y+vy/50,colour=speed),arrow = arrow(length = unit(0.1,"cm")),size=0.5)+scale_color_continuous(low = "grey80",high = "darkred")+facet_wrap(~z)
# 13.13绘制QQ图
# 使用qqnorm()和正态分布比较。给qqnorm()一个数值向量，在此基础上用qqline()加上理论分布
library(gcookbook)
# height的QQ图
qqnorm(heightweight$heightIn)
qqline(heightweight$heightIn)
# age的QQ图
qqnorm(heightweight$ageYear)
qqline(heightweight$ageYear)
# ggplot2有一个stat_qq()函数，但是没有提供画QQ线的简单方法
# 13.14绘制经验累积分布函数图
# 使用stat_ecdf()
library(gcookbook)
# heightIn的ecdf
ggplot(heightweight,aes(x=heightIn))+stat_ecdf()
# ageYear的ecdf
ggplot(heightweight,aes(x=ageYear))+stat_ecdf()
# ECDF表明了在观察数据中，小于或等于给定x值得观测所占的比例。因为是经验分布，所以累积分布线在每个有一个或者更多观测值得x值处产生一个阶梯。
# 13.15 创建马赛克图
# 使用vcd包里的mosaic()函数。
UCBAdmissions
# 显示“平铺”后的列联表
ftable(UCBAdmissions)
dimnames(UCBAdmissions)
# 使用mosaic()函数，输入的公式要包含在分割数据中使用的变量
install.packages("vcd")
library(vcd)
# 按照先Admit然后Gender再Dept的顺序分割数据
mosaic(~Admit+Gender+Dept,data = UCBAdmissions)
mosaic(~Dept+Gender+Admit,data=UCBAdmissions,highlighting = "Admit",highlighting_fill=c("lightblue","pink"),direction=c("v","h","v"))
# 另一种可能的分割方向
mosaic(~Dept+Gender+Admit,data = UCBAdmissions,highlighting = "Admit",highlighting_fill=c("lightblue","pink"),direction=c("v","v","h"))
# 这个顺序比较男和女不是很容易
mosaic(~Dept+Gender+Admit,data = UCBAdmissions,highlighting = "Admit",highlighting_fill=c("lightblue","pink"),direction=c("v","h","h"))
# mosaicplot是另外一个画马赛克图的函数
# 13.16 绘制饼图
# 使用pie()函数
library(MASS)
# 得到fold变量每个水平的频数
fold <- table(survey$Fold)
fold
# 画饼图
pie(fold)
# 提供一个命名的向量，或者一个只有数值的向量和一个标签向量
pie(c(99,18,120),labels = c("L on R","Neither","R on L"))
# 13.17 创建地图
# 从maps包里面获取地图数，用geom_polygon()（可以用颜色填充）或者geom_path()（不能填充）绘制
# 经度和纬度默认是画在直角坐标系中的，可以用coord_map()指定一个投影。
# 默认的投影是“mercator"（墨卡托投影）,和直角坐标系不一样，墨卡托投影中纬度线之间的距离会逐渐发生变化
library(maps)
# 美国地图数据
states_map <- map_data("state")
ggplot(states_map,aes(x=long,y=lat,group=group))+geom_polygon(fill="white",colour="black")
# geom_path(没有填充)和墨卡托投影
install.packages("mapproj")
ggplot(states_map,aes(x=long,y=lat,group=group))+geom_path()+coord_map("mercator")
map_data("state")
# 世界地图数据
world_map <- map_data("world")
world_map
sort(unique(world_map$region))
euro <- map_data("world",region = c("UK","France","Netherlands","Belgium"))
# Map region to fill color
ggplot(euro,aes(x=long,y=lat,group=group,fill=region))+geom_polygon(colour="black")+scale_fill_brewer(palette = "Set2")+scale_y_continuous(limits = c(40,60))+scale_x_continuous(limits = c(-25,25))
# 从世界地图中得到新西兰地图数据
nzl <- map_data("world",region = "New Zealand")
nzl <- subset(nzl,long>0&lat>-48) # 剔除岛屿
ggplot(nzl,aes(x=long,y=lat,group=group))+geom_path()
# 从新西兰（nz)地图中得到新西兰地图数据
nz2 <- map_data("nz")
ggplot(nz2,aes(x=long,y=lat,group=group))+geom_path()
# 查看mapdata包获得更多地图数据，查看map()函数，快速产生地图，查看mapproject获得地图的投影方法
# 13.18绘制等值区域图
# 把变量值和地图数据合并在一起，然后把一个变量映射到fill上
# 把USArrests数据集转换成正确的格式
crimes <- data.frame(state=tolower(rownames(USArrests)),USArrests)
crimes
library(maps)
states_map <- map_data("state")
# 合并数据集
crime_map <- merge(states_map,crimes,by.x = "region",by.y = "state")
# 合并之后，顺序发生了变化，可能会导致多边形位置不对，所以要对数据排序
head(crime_map)
library(plyr) # 为了使用arrange()函数
# 按group,order排序
crime_map <- arrange(crime_map,group,order)
head(crime_map)
# 画图，把一列数值映射到fill上
ggplot(crime_map,aes(x=long,y=lat,group=group,fill=Assault))+geom_polygon(colour="black")+coord_map("polyconic")
# 展示变量值如何从某个中点值向外发散的，可以用scale_fill_gradient2()
ggplot(crimes,aes(map_id = state,fill=Assault))+geom_map(map=states_map,colour="black")+scale_fill_gradient2(low = "#559999",mid = "grey90",high = "#BB650B",midpoint = median(crimes$Assault))+expand_limits(x=states_map$long,y=states_map$lat)+coord_map("polyconic")
# 找到分位数的边界
qa <- quantile(crimes$Assault,c(0,0.2,0.4,0.6,0.8,1.0))
qa
# 加入一个分位数类别的列
crimes$Assault_q <- cut(crimes$Assault,qa,labels = c("0-20%","20-40%","40-60%","60-80%","80-100%"),include.lowest = TRUE)
crimes
# 产生一个5个离散取值的调色板
pal <- colorRampPalette(c("#559999","grey80","#BB650B"))(5)
pal
ggplot(crimes,aes(map_id=state,fill=Assault_q))+geom_map(map = states_map,colour="black")+scale_fill_manual(values = pal)+expand_limits(x=states_map$long,y=states_map$lat)+coord_map("polyconic")+labs(fill="Assault Rate\nPercentile")
# 另一个画等值区域图（Choropleth Map)的方法是使用geom_map()，它不需要把变量取值和地图数据结合起来，因此可以比之前的方法更快地生成地图
# crimes中‘state'列要和state_map中的’region‘列匹配
ggplot(crimes,aes(map_id = state,fill=Assault))+geom_map(map=states_map)+expand_limits(x=states_map$long,y=states_map$lat)+coord_map("polyconic")
# 13.19创建空白背景的地图
# 首先保存下面的主题
# 创建一个去掉了很多背景元素的主题
theme_clean <- function(base_size=12){
require(grid) # unit()函数需要
theme_grey(base_size)%+replace%
theme(
axis.title = element_blank(),
axis.text = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0,"cm"),
axis.ticks.margin = unit(0,"cm"),
panel.margin = unit(0,"lines"),
plot.margin = unit(c(0,0,0,0),"lines"),
complete=TRUE
)
}
# 然后把它加在地图上。
ggplot(crimes,aes(map_id=state,fill=Assault_q))+geom_map(map = states_map,colour="black")+scale_fill_manual(values = pal)+expand_limits(x=states_map$long,y=states_map$lat)+coord_map("polyconic")+labs(fill="Assault Rate\nPercentile")+theme_clean()
# 13.20 基于空间数据格式（shapefile)创建地图
# maptools包中的readShapePoly()载入空间数据文件，用fortify()把数据转化为数据框的格式，然后画图
install.packages("maptools")
library(maptools)
# 载入空间数据并转化为数据框
uk_shp <- readShapePoly("GBR_adm/GBR_adm2.shp")
uk_map <- fortify(uk_shp)
ggplot(uk_map,aes(x=long,y=lat,group=group))+geom_path()
# Esri shapefile是一种很常用的地图数据，readShapePoly()函数读入空间数据文件，返回一个SpatialPolygonsDataFrame对象
uk_shp <- readShapePoly("GBR_adm/GBR_adm2.shp")
# 查看该对象的格式
str(uk_shp)
# 把它转为常规的数据框
uk_map <- fortify(uk_shp)
uk_map
# 实际上，也可以直接在ggplot()中输入SpatialPolygonsDataFrame对象，这里自动使用了fortify()：
ggplot(uk_shp,aes(x=long,y=lat,group=group))+geom_path()
# 数据集和其他空间数据集在http://www.gadm.org/网站上下载。
扫一扫
04-12
3131
3131
03-20
1859
1859
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
