Paul Maier Research

Diversity Interpolation

Paul Maier — Apr 22, 2015, 10:05 PM

# Some examples of 3D plots
# In this case our 3D "image" is an interpolation
# of gene diversity (π) in Yosemite Toads, YNP
# This is just for fun but feel free to leave a comment

library(raster)

Loading required package: sp

library(rgdal)

rgdal: version: 0.8-16, (SVN revision 498)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj

library(rgl)
library(plot3D)

Warning: package 'plot3D' was built under R version 3.1.2

YOSEdiv <- raster("~/Desktop/Scratch/YOSEKINGS 4/YOSEpi.tif")
YOSEsurf <- as.matrix(YOSEdiv)
SEKIdiv <- raster("~/Desktop/Scratch/YOSEKINGS 4/Kingspi2.tif")
SEKIsurf <- as.matrix(SEKIdiv)

persp(YOSEsurf, col = "green", phi=30,theta=60, ticktype="simple", box=T, 
      xlab='Long', ylab='Lat', zlab='Average π')

persp3D(z=YOSEsurf, phi=30,theta=60, xlab='Long', ylab='Lat', zlab='Average π')

persp3D(z = volcano, zlim = c(-60, 200), phi = 20,
        colkey = list(length = 0.2, width = 0.4, shift = 0.15,
        cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90,
        xlab = "latitude",ylab = "longitude", zlab = "Average π",
        clab = c("","height","m"), bty = "f", plot = FALSE)

persp3D(z = YOSEsurf, col = "white", shade = 0.1, plot = FALSE,
        xlab = "latitude",ylab = "longitude", zlab = "Average π")
contour3D(z = YOSEsurf, colvar = YOSEsurf, lwd = 2,
          add = TRUE, dDepth = 0.1, col = "black")

persp3D(z = YOSEsurf, zlim = c(0.06, 0.16), col = "white", 
        xlab = "latitude",ylab = "longitude", zlab = "Average π",
        shade = 0.1, d = 2, plot = FALSE)
contour3D(z = YOSEsurf, colvar = YOSEsurf, lwd = 2, add = TRUE, 
          nlevels = 20, clab = c("Average", "π"), plot = FALSE,
          colkey = list(at = seq(90, 190, length.out = 5)))
contour3D(z = 0.16, colvar = YOSEsurf, lwd = 2, col = "grey",
          add = TRUE, nlevels = 5)

#rect(-100, 20, -40, 80, lwd = 3)
zlim <- c(0.07,0.14)
# Actual bathymetry, 4 times increased resolution, with contours
persp3D(z = YOSEsurf, xlab = "latitude", bty = "bl2",
        ylab = "longitude", zlab = "Average π", clab = "Yosemite π IDW Interpolated",
        expand = 1, phi = 20, theta = 35, resfac = 2, cex.lab=1,
        contour = list(col = "grey", side = c("zmin", "z")),
        zlim = zlim, colkey = list(side = 1, length = 0.25))

zlim <- c(0.07,0.24)
#ylim <- c(-1,2)
# Actual bathymetry, 4 times increased resolution, with contours
persp3D(z = SEKIsurf, xlab = "latitude", bty = "bl2",
        ylab = "longitude", zlab = "Average π", clab = "Kings π IDW Interpolated",
        expand = 1, d = 20, phi = 20, theta = 35, resfac = 2, cex.lab=2.5,
        contour = list(col = "grey", side = c("zmin", "z")),
        zlim = zlim, #ylim = ylim,
        colkey = list(side = 1, length = 0.25))