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

刺骨的言语ヽ痛彻心扉 2021-12-19 23:03 1030阅读 0赞

## 数据读取 ##

data1 <- read.table(‘forams-Plio.txt’,sep = ‘\\t’,header = TRUE)  
head(data1)  
![在这里插入图片描述][20190630235232889.png]  
str查看数据结构后，知道数据是242个，那么说明高是242个高，那咱们看242能被2和121进行相乘。涉及的矩阵为2成121。那么我想设计成30乘以8的矩阵，就可以删掉2个元素，暂时只是删掉一个。如下所示。  
![在这里插入图片描述][20190630235316268.png]

## \#设置高程值 ##

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)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70]

## 绘制等高线图二 ##

数据描述：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 1]

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)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 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)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 3]

##  ##

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 4]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 5]

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 6]

[20190630235232889.png]: /images/20211219/46be2f6d69ba439a8e9762b65e74629e.png
[20190630235316268.png]: /images/20211219/f6e61c39443a43a482291b2a9620e016.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70]: /images/20211219/037950c6fcdb4c6dac7e5c98e86e4f7d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 1]: /images/20211219/4ca180a5831b405291b2cf9d15d2b348.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 2]: /images/20211219/f512489a7997490ea13f9950b17fa33c.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 3]: /images/20211219/f5e5f8d885c84a09a9b57be5082bd0a7.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 4]: /images/20211219/ec678262bd264bcc9b63bcc8a557dbbf.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 5]: /images/20211219/f5b733989a8246caba51ba9d4fa2d2d9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3RhbmRlbGlu_size_16_color_FFFFFF_t_70 6]: /images/20211219/3f5401bd6ca14f5892ae72831dad904f.png