R语言绘制等值线和等高线

江南才子

　　
数据读取
　　data1 <- read.table(‘forams-Plio.txt’,sep = ‘\t’,header = TRUE)
head(data1)

str查看数据结构后，知道数据是242个，那么说明高是242个高，那咱们看242能被2和121进行相乘。涉及的矩阵为2成121。那么我想设计成30乘以8的矩阵，就可以删掉2个元素，暂时只是删掉一个。如下所示。


#设置高程值
　　str(data1)
zz=data1[-1,]
str(zz)
z1=zz$NumSp
z=matrix(z1, 10, 10)
str(z)
str(x)
#x坐标，每个网格为10米分辨率，方向由南向北
x <- 10*(1:nrow(z))
#y坐标，每个网格为10米分辨率，方向由东向西
y <- 10*(1:ncol(z))

二维可视化：栅格+等高线
　　par(mar=rep(0.5,4))
image(x, y, z, col=terrain.colors(100), axes=F)
contour(x, y, z, levels=seq(from=min(z), to=max(z), by=10),axes=T, add=T)
library(png)
logo <- readPNG(‘E:/logo.png’)
　　layout(z,heights=x,widths=y)

绘制等高线图二
　　数据描述：

install.packages("readxl")
install.packages("akima")
install.packages("fields")
install.packages("maps")
install.packages("spam")
install.packages("maps")
library(readxl)
library(akima)
library(fields)
library(maps)
library(spam)
data=read_excel("C:\\Users\\Administrator\\Desktop\\data.xlsx")
colnames(data) <- c("V1","V2","V3")
mydata=as.data.frame(data)
#数据转化 使用akima 包的 interp，内插数据
s <- interp(mydata$V1,mydata$V2,mydata$V3)
#绘制图像与图形主标题，副标题，计算平均值，Uniformity，极差
image.plot(s,main=
paste("Rs Mean is",round(mean(mydata$V3), 2),
", RSU= SD*100/MEAN%=",round(sd(mydata$V3)*100/mean(mydata$V3), 2),"%",
"\n","Max = ",round(max(mydata$V3), 2),"Min = ",round(min(mydata$V3), 2),"Range = Max-Min = ",
round(max(mydata$V3)-min(mydata$V3), 2)),
sub=paste(date(),", Map Created By Stone_Hou @ R"),
cex=2,xlab="X coordinate / mm",ylab="Y coordinate / mm",col.main = "#DD3333",
col.sub="gray",horizontal=TRUE)
#循环输出49个坐标点的图像add x,y point,
for(i in 1:49){
points(mydata$V1[i],mydata$V2[i],type="p", col="blue",bg = "blue", cex = 1.5)
}
#使用CircleSwitchSign来控制是否输出同心圆
CircleSwitchSign <- 0
if(CircleSwitchSign==1){
#循环输出五个半径分别为0,49,98,147,150的同心圆，包括原点
for(r in c(0,49,98,147,150)){
k = 1:360
x=r*sin(k/180*pi)
y=r*cos(k/180*pi)
lines(x,y,col="#444444", lwd=1.5)
}
}
#循环输出49个点的测试值，添加数据到每一个坐标点，偏移5个单位
for(i in 1:length(mydata$V3)){
info = paste(i,"_",round(mydata$V3[i],2))
xPoint=mydata$V1[i]
yPoint=mydata$V2[i]
k = 1:360
r = 150
#print(r)
x=r*sin(k/180*pi)
y=r*cos(k/180*pi)
#lines(x,y,col="red")
text(xPoint,yPoint,info,pos=3,cex=0.6, col = "black")
#text(xPoint,yPoint,paste(percent,"%"),cex=3)
#add text value
#text(x-10, y-10,mydata$V3, cex=0.8)
}
#使用SwitchSign来控制是否输出网格线
SwitchSign <- 0
if(SwitchSign==1){
#输出网格线，grid line
rx <- range(x <- mydata$V1)
ry <- range(y <- mydata$V2)
abline(h = pretty(rx, 10), v = pretty(ry, 10), col = "gray")
abline(h = pretty(rx, 10), v = pretty(ry, 10), col = "#999999")
## the x=0- and y=0-axis,
abline(h=0, v=0, col = "#555555",lwd=2)
## the x- and y-axis, and an integer grid
#abline(h=0, v=0, col = "gray60",lwd=5)
#text(1,0, "abline( h = 0 )", col = "gray60", adj = c(0, -.1))
#abline(h = -1:5, v = -2:3, col = "lightgray", lty=3)
#abline(a=1, b=2, col = 2)
#text(1,3, "abline( 1, 2 )", col=2, adj=c(-.1,-.1))
}
#使用contourSwitchSign来控制是否输出二维等高线图add contour map
contourSwitchSign <- 1
if(contourSwitchSign==1){
contour(s,add=TRUE,cex=1.5)
}
#关闭图像输出设备
#dev.off()
#lines(mydata$V3, col="red", type="o", pch=22, lty=2)

等温线图

library(MASS)
Sigma <- matrix(c(1,0.7,0.7,1),2,2)
Sigma
r=mvrnorm(n=1000, c(0,3), Sigma)
par(mfrow = c(2, 2))
#density plot
plot(density(r))
#kernel density estimate
bivn.kde <- kde2d(r[,1], r[,2], n = 50)
#perspective plot
persp(bivn.kde, phi = 45, theta = 30)
#contour plot
contour(bivn.kde)
#contour plot with image
image(bivn.kde,col = terrain.colors(100)); contour(bivn.kde, add = T)

library(MASS)
x<-rnorm(20);y<-rnorm(20)
contour(kde2d(x,y), col = 'red', drawlabel=F)
points(x,y)

x <- seq(-pi, pi, length = 50)
y <- x
f <- outer(x, y, function(x, y)cos(y)/(1+x^2))
contour(x, y, f)

x <- seq(-pi, pi, length = 50)
y <- x
f <- outer(x, y, function(x, y)cos(y)/(1+x^2))
contour(x, y, f, nlevels = 45)
fa <- (f-t(f))/2