##########################################

#### 新西兰的芒格法奥火山地形可视化   ####

##########################################

#首先加载volcano.csv数据，并且设置没有列头名称

volcanodem <- read.csv("E:/volcanodem", header=FALSE)

#将数据集转换为矩阵

volcano <- as.matrix(volcanodem)

#设置高程值

z <- volcano

#x坐标，每个网格为10米分辨率，方向由南向北

x <- 10*(1:nrow(z))   

#y坐标，每个网格为10米分辨率，方向由东向西

y <- 10*(1:ncol(z))    

## 二维可视化：栅格+等高线（见图片1）

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=F, add=T)

## 三维可视化（见图片2）

par(mar=rep(0,4))    

#三维网格模式显示（白膜）

persp(x,y,z,theta=120,phi=15,scale=F,axes=F)

## 带色带的三维模式渲染（见图片3）

#设置高程值为2倍拉伸

z <- 2 * volcano

x <- 10 * (1:nrow(z))

y <- 10 * (1:ncol(z))

## 创建一份新的用于绘制的网格数据集

#z0是用于设置栅格边缘的颜色

z0 <- min(z) - 20

#用z0的颜色，把整个栅格包裹起来

z <- rbind(z0, cbind(z0, z, z0), z0)

#用x的最小值和最大值，把x包裹起来，对应上面那个z0

x <- c(min(x) - 1e-10, x, max(x) + 1e-10)

#用y的最小值和最大值，把x包裹起来，对应上面那个z0

y <- c(min(y) - 1e-10, y, max(y) + 1e-10)

## 创建用于显示颜色的矩阵

#默认全部使用绿色

fcol <- matrix("green3", nr = nrow(z)-1, nc = ncol(z)-1)

#用灰色把所有的绿色都包裹起来,即设置四周的边界值

fcol[ , i2 <- c(1,ncol(fcol))] <- "gray"

fcol[i1 <- c(1,nrow(fcol)) , ] <- "gray"

## Take average of four neighboring values for palette

##将上面设置的默认色，用都取相邻的四个格网颜色的平均值进行替换

zi <- (volcano[ -1,-1] + volcano[ -1,-61] + volcano[-87,-1] + volcano[-87,-61])/4

pal <- terrain.colors(20)[cut(zi, quantile(zi, seq(0,1, len = 21)), include.lowest = TRUE)]

fcol[-i1,-i2] <- pal

## 绘图

par(mar=rep(0,4))

persp(x, y, z, theta=120, phi=15, col = fcol, scale = FALSE, shade = 0.4, border = NA)