地形图绘制

作者简介

勾蒙蒙，R语言资深爱好者。

个人公众号： R语言及生态系统服务。

前文传送门：

脏数据-数据量纲差异

在写论文或者是做学术PPT的时候，多数人都倾向于将科研图片做出不一样的效果，集“高大上”于一体，一方面能给受众以视觉冲击，另一方面漂亮的图片确实会给你的论文或者PPT增色，即使是最简单的样地地形数据图。

下图是我从一些论文上看到的关于样地地形情况或者采样点情况的介绍图，其实像实现这些图并不难，在ArcGis中当然能够实现，但有些许的复杂，本文着重对样地地形图在R语言中的实践加以详细介绍，以供大家互相交流学习！

640?wx_fmt=jpeg

ready go

##加载程序包

install.packages("sp")

library(sp)

##生成样方中心店海拔数据

XY<-expand.grid(x=seq(0,90,10),y=seq(0,90,10))

setwd("C:/Users/lx/Desktop")

Elevation<-read.csv("Elevation.csv")

Data<-cbind(XY,Elevation)

##排序

newdata<-Data[order(Data$x,-Data$y),]

##数据转化为矩阵

newdata1<-matrix(newdata$Elevation,ncol=10,nrow=10,byrow =FALSE)

一、简单的点图

par(mar=c(4,4,2,2))

coordinates(Data)<-c("x","y") ##设定坐标

plot(Data,xlab="横坐标",ylab="纵坐标")

axis(1,at=seq(0,100,10),labels=seq(0,100,10),line=0,tick=F)

axis(2,at=seq(0,100,10),labels=seq(0,100,10),line=0,tick=F)

box(lwd=2,col="red")

640?wx_fmt=jpeg

##生成点图

spplot(Data,"Elevation")

640?wx_fmt=jpeg

二、简单的面图

##矩阵图

install.packages("psych")

library(psych)

label<-seq(0,90,10)

cor.plot(newdata1,number=TRUE,show.legend=TRUE,zlim=c(70,130),labels=label,n.legend=5)#n.legend设置图例显示的级别

640?wx_fmt=jpeg

三、使用plot3D绘图

install.packages("plot3D")

library(plot3D)

##2D图

par(mfrow = c(3, 3), mar = c(3, 3, 3, 2))

contour2D(newdata1, lwd = 2, colkey = FALSE)#colkey设置图例TRUE/FALSE

contour2D(newdata1, lwd = 2)#colkey默认显示图例

image2D(newdata1, clab = "m")#clab设置图例的名称

image2D(newdata1, shade = 0.4)#shade设置阴影程度

image2D(newdata1, facets = FALSE)#facets设置每个小块是否填充颜色

image2D(newdata1, contour = TRUE)#contour在图中添加等高线

image2D(newdata1, rasterImage = TRUE, contour = list(lwd = 2, col = jet.col(11)))#设置等高线颜色

image2D(newdata1, theta = 30, NAcol = "black")#theta设置图形方位；NAcol设置空值的颜色，在image2D中其默认颜色为白色，在contour2D中默认为无色

image2D(newdata1, lighting = TRUE, rasterImage = TRUE,contour = list(col = "white", labcex = 0.8, lwd = 3, alpha = 0.5))#lighting设置图形被照亮的效果；rasterImage为TRUE，则会显示为栅格模式

640?wx_fmt=jpeg

##3D图

par(mfrow = c(3, 3), mar = c(2, 2, 2, 2))

persp(newdata1,expand=0.25)

persp(newdata1, theta = 40, phi = 40, col = "gold", border = NA, shade = 0.5,expand=0.25)#expand调整膨胀率;border设置表面线条颜色，默认为NULL，若为NA，则表面线条的绘制命令失效，

persp3D(z = newdata1, clab = "m",expand=0.25,border = NULL)

persp3D(z = newdata1, clab = "m", shade = 0.2,expand=0.25)

persp3D(z = newdata1, facets = FALSE,expand=0.25)

persp3D(z = newdata1, facets = FALSE, curtain = TRUE,expand=0.25)

persp3D(z = newdata1, col = "white", shade = 0.5,expand=0.25)

persp3D(z = newdata1, col = ramp.col(c("white", "black")), border = "red",expand=0.25)#border设置表面线条颜色，col设置表面的填充颜色

persp3D(z = newdata1, facets = FALSE, col = "darkblue",expand=0.25)

640?wx_fmt=jpeg

##背景和坐标轴设置

par(mfrow = c(3, 3), mar = c(1, 1, 1, 1))

persp3D(z = newdata1, col = "lightblue", shade = 0.5,expand=0.25)

persp3D(z = newdata1, col = "lightblue", shade = 0.5, box = FALSE,expand=0.25)#box = FALSE边框去除

persp3D(z = newdata1, col = "lightblue", shade = 0.5, ticktype = "detailed",expand=0.25)#ticktype设置坐标轴刻度,如为sample则只用箭头表示增加的方向，如第一个小图，如为detailed则跟2D图一样绘制对应的刻度

##bty介绍：f：显示全部边框；"b"：一般默认，只显示背面；"b2"：背面和网格线都显示；"g"：灰色背景白色的网格线；"bl"：黑色的背景；"bl2"黑色的背景灰色的网格线;"n"：不显示背景也不显示网格线；"u"：自己可以设置背景和网格线

persp3D(z = newdata1, col = "lightblue", shade = 0.5, bty = "f",expand=0.25)

persp3D(z = newdata1, col = "lightblue", shade = 0.5, bty = "b2",expand=0.25)

persp3D(z = newdata1, col = "lightblue", shade = 0.5, bty = "g",expand=0.25)

persp3D(z = newdata1, col = "lightblue", shade = 0.5, bty = "bl",expand=0.25)

persp3D(z = newdata1, col = "lightblue", shade = 0.5, bty = "u",col.panel = "yellow", col.grid = "grey",,expand=0.25)#col.panel设置背景颜色；col.grid设置网格线颜色

persp3D(z = newdata1, col = "lightblue", shade = 0.5,ticktype = "detailed", bty = "b2",,expand=0.25)

640?wx_fmt=jpeg

##调整视图和角度

par(mfrow = c(3, 2), mar = c(1, 1, 1, 1))

x <- 1:nrow(newdata1)

y <- 1:ncol(newdata1)

persp3D(x,y, z = newdata1, col = "lightblue", scale = FALSE,shade = 1, expand =0.05)#scale决定是否归一化

persp3D(x, y , z = newdata1, col = "lightblue", scale = FALSE,shade = 0.5, expand = 0.05, d =0.8)#d值调整透视的强度，d大于1将减少透视效果，小于1增强透视效果

persp3D(x, y , z = newdata1, col = "lightblue", scale = FALSE,shade = 0.5, expand = 0.05, r =10)#r值调整观察点到图中心点距离，很显然距离越远就越接近轴测图,距离越近透视效果越强烈

persp3D(x, y , z = newdata1, col = "lightblue", scale = FALSE,shade = 0.5, expand = 0.05, theta = 10,phi=-20)#定义视角的角度theta为方位,phi为角度

persp3D(x, y , z = newdata1, col = "lightblue", scale = FALSE,shade = 0.5, expand = 0.05, ltheta = 10,lphi=90)#如果ltheta和lphi为特定值，则图的表面被亮光着色，即在theta方位，phi角度上打光的效果

640?wx_fmt=jpeg

##绘制彩带图及直方图

par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

ix <- seq(1, nrow(newdata1), length.out = 20)

iy <- seq(1, ncol(newdata1), length.out = 20)

ribbon3D(z = newdata1[, iy],expand=0.25)#沿着列的彩带

ribbon3D(z = newdata1[ix, ], along = "y",curtain = TRUE, space = 0.9, shade =0.3,expand=0.25)#space设置间隔，其范围为0.1-0.9;

ribbon3D(z = newdata1[ix, iy], along = "xy",expand=0.25)#xy交织的彩带

hist3D(z = newdata1[ix,iy], shade = 0.1,expand=0.8)#直方图

640?wx_fmt=jpeg

#colkey设置

par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

persp3D(z = newdata1, log = "c", clab = c("m (Elevation)"),expand=0.2)#clab设置图例的名称

persp3D(z = newdata1, clab ="m",colkey = list(side = 4, length = 1, width = 0.5, cex.axis = 2),expand=0.2)#side设置放置位置（1234对应下左上右）；length设置长度,width设置宽度；cex.axis设置标签的大小

persp3D(z = newdata1, clab = c("height", "m"),colkey = list(length = 0.5, shift = -0.1),expand=0.2)

par(mar = c(4, 4, 2, 2))

image2D(z = newdata1, clab = "height, m",colkey = list(dist = -0.20, shift = 0.15,side = 3, length = 0.5, width = 0.5,cex.clab = 1.2, col.clab = "white", line.clab = 2,col.axis = "white", col.ticks = "white", cex.axis = 0.8))

640?wx_fmt=jpeg

##组合persp3D\image\contour

par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

ribbon3D(z = newdata1, zlim = c(-100, 200), image = TRUE)#彩带图+2D图

persp3D(z = newdata1, zlim = c(-100, 200), contour = TRUE)#3D图+地形图

persp3D(z = newdata1, zlim= c(-200, 200), phi = 30,contour = list(nlevels = 20, col = "red"),image = list(col = grey (seq(0, 1, length.out = 100))))#大满贯：3D图+2D图+地形图

persp3D(z = newdata1, contour = list(side = c("zmax", "z")), zlim= c(90, 300),

phi = 30, theta = 20, d = 10, box = FALSE)#3D图+地形图，地形图的位置2个，side设置地形图的位置，即side在Elevetion这一列滑动的位置

640?wx_fmt=jpeg

##更复杂的多图组合

x<-1:nrow(newdata1)

y<-1:ncol(newdata1)

#首先设置一个3D图，设置contour在3D图的下面

persp3D(x,y,z=newdata1,theta=10,phi=20,expand=0.01,box=FALSE,scale=FALSE,contour=TRUE,zlim=c(-100,200),clim=range(newdata1),plot=TRUE)

#在Elevation=100的地方添加一个2D图

image3D( x, y, z = 100, add = TRUE, clim = range(newdata1),colvar = newdata1, colkey = FALSE, facets = NA,plot=TRUE)

#3D地形图的添加

contour3D(x, y, z = 100,add = TRUE, clim = range(volcano), lwd = 3,colvar = newdata1, colkey = FALSE,plot=TRUE)

640?wx_fmt=jpeg

　往期精彩内容整理合集　

2017年R语言发展报告（国内）

R语言中文社区历史文章整理（作者篇）

R语言中文社区历史文章整理（类型篇）

640?wx_fmt=jpeg

公众号后台回复关键字即可学习

回复 R                  R语言快速入门及数据挖掘
回复 Kaggle案例  Kaggle十大案例精讲（连载中）
回复文本挖掘   手把手教你做文本挖掘
回复可视化   R语言可视化在商务场景中的应用
回复大数据         大数据系列免费视频教程
回复量化投资      张丹教你如何用R语言量化投资
回复用户画像      京东大数据，揭秘用户画像
回复数据挖掘     常用数据挖掘算法原理解释与应用
回复机器学习人工智能系列之机器学习与实践
回复爬虫            R语言爬虫实战案例分享