RG #110: 3D scatter plot with multiple series in Y axis

X = seq(1, 100, 5)
Y = seq (1, 100, 5)
Z = rnorm (length (X), 10, 2)
data1 <- data.frame (X, Y, )
data2 <- data.frame (X, Y, Z1 = Z - 5)
data3 <- data.frame (X, Y, Z1 = Z - 3)

require(scatterplot3d)
s3d <- scatterplot3d(data1, color = "blue", pch = 19, xlim=NULL, ylim=NULL, zlim= c(0, 20))
s3d$points3d(data2, col = "red", pch = 18)
s3d$points3d(data3, col = "green4", pch = 17)

 

RG#104: Rootogram plot

require(latticeExtra)
require(lattice)

library(lattice)

set.seed(1234)
x1 <- rpois(1000, lambda = 25)

plt <- rootogram(~x1, dfun = function(x) dpois(x, lambda = 25))
plt 

 

 

RG#109:small plot(s) with in a big plot

require(ggplot2)
library(gridBase)
plot(cos, -pi, 2*pi, ylim = c(-1.3, 1.5), col = "red")
myd <- data.frame (X = 1:10, Y = c(3, 4, 8, 7, 2, 1, 9, 4, 2, 3))
qp <- qplot(X, Y, data=myd) + theme_bw()
print(qp, vp=viewport(.65, .65, .25, .25))

 library(lattice)
library(gridBase)
plot.new()
pushViewport(viewport())
set.seed(1234)
xvars <- rnorm(25, 5, 1)
yvars <- rnorm(25, 5, 1)
xyplot(yvars~xvars,  xlim = c(0, 10), ylim = c(0, 10) )
pushViewport(viewport(x=.6,y=.85,width=.20,height=.15,just=c("left","top")))
grid.rect()
par(plt = gridPLT(), new=TRUE)
plot(xvars,yvars)
popViewport(2)

 

RG#107: Plot 3d horizontal lines (bars) over map (world and US example)

library("maps")
require(ggplot2)
library(ggsubplot)
world.map <- map("world", plot = FALSE, fill = TRUE)
world_map <- map_data("world")
require(lattice)
require(latticeExtra)

 # Calculate the mean longitude and latitude per region (places where subplots are plotted)

library(plyr)
cntr <- ddply(world_map,.(region),summarize,long=mean(long),lat=mean(lat))


# example data
 myd <- data.frame (region = c("USA","China","USSR","Brazil", "Australia","India", "Nepal", "Canada",
                                "South Africa", "South Korea", "Philippines", "Mexico", "Finland",
                                 "Egypt", "Chile", "Greenland"),
               frequency = c(501, 350, 233, 40, 350, 150, 180, 430, 233, 120, 96, 87, 340, 83, 99, 89))



subsetcntr  <- subset(cntr, region %in% c("USA","China","USSR","Brazil", "Australia","India", "Nepal", "Canada",
                                "South Africa", "South Korea", "Philippines", "Mexico", "Finland",
                                 "Egypt", "Chile", "Greenland"))

simdat <- merge(subsetcntr, myd)
colnames(simdat) <- c( "region","long","lat", "myvar" )


panel.3dmap <- function(..., rot.mat, distance, xlim,
     ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled) {
       scaled.val <- function(x, original, scaled) {
      scaled[1] + (x - original[1]) * diff(scaled)/diff(original)
     }
       m <- ltransform3dto3d(rbind(scaled.val(world.map$x,
           xlim, xlim.scaled), scaled.val(world.map$y, ylim,
          ylim.scaled), zlim.scaled[1]), rot.mat, distance)
        panel.lines(m[1, ], m[2, ], col = "green4")
      }


p2 <- cloud(myvar ~ long + lat, simdat, panel.3d.cloud = function(...) {
         panel.3dmap(...)
          panel.3dscatter(...)
 }, type = "h", col = "red", scales = list(draw = FALSE), zoom = 1.1,
            xlim = world.map$range[1:2], ylim = world.map$range[3:4],
          xlab = NULL, ylab = NULL, zlab = NULL, aspect = c(diff(world.map$range[3:4])/diff(world.map$range[1:2]),
          0.3), panel.aspect = 0.75, lwd = 2, screen = list(z = 30,
          x = -60), par.settings = list(axis.line = list(col = "transparent"),
            box.3d = list(col = "transparent", alpha = 0)))
 
print(p2)

# Over US map
library("maps")
state.map <- map("state", plot = FALSE, fill = FALSE)
require(lattice)
require(latticeExtra)

 # data
 state.info <- data.frame(name = state.name, long = state.center$x,
      lat = state.center$y)

set.seed(123)
state.info$yvar<- rnorm (nrow (state.info), 20, 5)


panel.3dmap <- function(..., rot.mat, distance, xlim,
     ylim, zlim, xlim.scaled, ylim.scaled, zlim.scaled) {
       scaled.val <- function(x, original, scaled) {
      scaled[1] + (x - original[1]) * diff(scaled)/diff(original)
     }
       m <- ltransform3dto3d(rbind(scaled.val(state.map$x,
           xlim, xlim.scaled), scaled.val(state.map$y, ylim,
          ylim.scaled), zlim.scaled[1]), rot.mat, distance)
        panel.lines(m[1, ], m[2, ], col = "grey40")
      }

pl <- cloud(yvar ~ long + lat, state.info, subset = !(name %in%
       c("Alaska", "Hawaii")), panel.3d.cloud = function(...) {
         panel.3dmap(...)
          panel.3dscatter(...)
 }, col = "blue2",  type = "h", scales = list(draw = FALSE), zoom = 1.1,
            xlim = state.map$range[1:2], ylim = state.map$range[3:4],
          xlab = NULL, ylab = NULL, zlab = NULL, aspect = c(diff(state.map$range[3:4])/diff(state.map$range[1:2]),
          0.3), panel.aspect = 0.75, lwd = 2, screen = list(z = 30,
          x = -60), par.settings = list(axis.line = list(col = "transparent"),
            box.3d = list(col = "transparent", alpha = 0)))
 print(pl)

RG#106: add satellite imagery, or open street maps to your plots using openmap package (bing, mapquest)

library(OpenStreetMap)
library(rgdal)




# get world map
map <- openmap(c(70,-179), c(-70,179))

plot(map)




bingmap <- openmap(c(70,-179), c(-70,179), type = "bing")

plot(bingmap)







# types available 

# type = c("osm", "osm-bw", "maptoolkit-topo", "waze", "mapquest", "mapquest-aerial", "bing", "stamen-toner", "stamen-terrain", "stamen-watercolor", "osm-german", "osm-wanderreitkarte", "mapbox", "esri", "esri-topo", "nps", "apple-iphoto", "skobbler", "cloudmade-<id>", "hillshade", "opencyclemap", "osm-transport", "osm-public-transport", "osm-bbike", "osm-bbike-german")

#zoom maps, plot a portion

# upperLeft, lowerRight

lat <- c(43.834526782236814, 30.334953881988564)

lon <- c(-85.8857421875, -70.0888671875)

southest <- openmap(c(lat[1],lon[1]),c(lat[2],lon[2]),zoom=7,'osm')

plot(southest) 










library(UScensus2000tract)

data(south_carolina.tract)

lat <- c(35.834526782236814, 30.334953881988564)

lon <- c(-85.8857421875, -70.0888671875)

southest <- openmap(c(lat[1],lon[1]),c(lat[2],lon[2]),zoom=7,'osm')

south_carolina.tract <- spTransform(south_carolina.tract,osm())

plot(southest)
plot(south_carolina.tract,add=TRUE,col=(south_carolina.tract@data$med.age>35)+4)







plot(southest)
plot(south_carolina.tract,add=TRUE,col=(south_carolina.tract@data$med.age>55)+4)












RG#105: Line range or Crossbar plot

set.seed (1234)
require(ggplot2)
Xv <- c(rnorm (500, 10,3), rnorm (500, 50, 20), rnorm (500, 70, 20))
Yv <- c(rnorm (500, 10,3), rnorm (500, 70, 5), rnorm (500, 30, 5))
myd <- data.frame (Xv, Yv )

m <- ggplot(myd, aes(x = Xv, y = Yv)) +
  geom_point() + geom_density2d() + theme_bw()
  
m


myd <- data.frame (Xv = c("A", "B", "C", "E", "F"), Yv = c( 35, 35, 69, 15, 60),
                    errs = c(5, 10, 3, 5, 12))

se <- ggplot(myd, aes(Xv, Yv,
  ymin = Yv - errs, ymax= Yv + errs, colour = Xv))
  
# line range plot 

se + geom_linerange() + theme_bw()

se + geom_linerange()  + coord_flip() + theme_bw()

se + geom_crossbar(width = 0.5) + theme_bw()

RG#104: 2d density plots

require(ggplot2)

set.seed (1234)
Xv <- c(rnorm (500, 10,3), rnorm (500, 50, 20), rnorm (500, 70, 20))
Yv <- c(rnorm (500, 10,3), rnorm (500, 70, 5), rnorm (500, 30, 5))
myd <- data.frame (Xv, Yv )

m <- ggplot(myd, aes(x = Xv, y = Yv)) +
  geom_point() + geom_density2d() + theme_bw()
  
m

RG#103: Combing different types of plot in trellis type

require(lattice)

require(latticeExtra)

#xy plot, conditioned from quakes data, with histogram 
qk <- xyplot(lat ~ long | equal.count(depth, 2), quakes,
    aspect = "iso", pch = "+", cex = 1.5, xlab = NULL, ylab = NULL)
qk <- c(qk, depth = histogram(quakes$depth), layout = c(3, 1))

update(qk, scales = list(at = list(NA, NA, NA, NA),
                          y = list(draw = FALSE)))





# suppress scales on the first 2 panels
update(qk, scales = list(at = list(NULL, NULL, NA),
                          y = list(draw = FALSE)))

RG#102: Double Y axis trellis plot (weather data example)

require(latticeExtra)

require (lattice)

data(SeatacWeather)

tempatures <- xyplot(min.temp + max.temp ~ day | month,
               data = SeatacWeather, type = "l", layout = c(3, 1))
rainfall <- xyplot(precip ~ day | month, data = SeatacWeather, type = "h", lwd = 4)

doubleYScale(tempatures, rainfall, style1 = 0, style2 = 3, add.ylab2 = TRUE,
   text = c("min. T", "max. T", "rain"), columns = 3)

RG#101: Plot country map with annotation of regions

## Load required packages
library(maptools)
library(raster)

 # Download data from gadm.org using getData function 
admNPL <- getData('GADM', country='NPL', level=3)
plot(admNPL, bg="lightgreen", axes=T)

##Plot 
plot(admNPL, lwd=10, border="skyblue", add=T)
plot(admNPL,col="yellow", add=T)
grid()
box()

invisible(text(getSpPPolygonsLabptSlots(admNPL), labels=as.character(admNPL$NAME_3), cex=0.4, col="black", font=2))
mtext(side=3, line=1, "District Map of Nepal", cex=2)
mtext(side=1, "Longitude", line=2.5, cex=1.1)
mtext(side=2, "Latitude", line=2.5, cex=1.1)

RG#100: Trellis map plot with heatmap colors

require(maps)
library(mapproj)

worldmap <- map('world', plot = FALSE, fill = FALSE,  projection = "azequalarea")
country = worldmap$names
set.seed(1234)
var.2010 = rnorm (length (country), 20, 10)
var.2011 = var.2010*1.1 + rnorm (length (country), 0, 1)
var.2012 = var.2011*0.98 + rnorm (length (country), 0, 4)
var.2013 = var.2011*0.98 + rnorm (length (country), 0, 30)
worldt <- data.frame (country, var.2010, var.2011, var.2012, var.2013)

mapplot(country ~ var.2011, worldt, map = map("world",     plot = FALSE, fill = TRUE))

mapplot(country ~ var.2010 + var.2011 + var.2012 + var.2013, data = worldt, map = map("world",     plot = FALSE, fill = TRUE))

# trellis plot for country maps not available in maps package:

 require(maptools)
# get the map; may need sometime to be loaded 
con <- url("http://gadm.org/data/rda/NPL_adm3.RData")
print(load(con))
close(con)

# from your data file working directory 
## load ("NPL_adm3.RData")

# data 
districts = gadm$NAME_3
set.seed(1234)
var1 <- rnorm (length (districts), 100, 30)
var2 <- rnorm (length (districts), 100, 30)
 myd <- data.frame (districts, var1, var2)    

# US county level map 

uscountymap <- map('county', plot = FALSE, fill = FALSE,  projection = "azequalarea")
county = uscountymap$names
set.seed(1234)
var.2010 = rnorm (length (county), 50, 10)
var.2011 = var.2010*1.1 + rnorm (length (county), 0, 5)
var.2012 = var.2011*0.98 + rnorm (length (county), 0, 10)
var.2013 = var.2011*1.2 + rnorm (length (county), 0, 15)
uscounty <- data.frame (county, var.2010, var.2011, var.2012, var.2013)

mapplot(county ~ var.2010 + var.2011 + var.2012 + var.2013, data = uscounty, map = map("county",    plot = FALSE, fill = TRUE))

# US state level map 
usstmap <- map('state', plot = FALSE, fill = FALSE,  projection = "azequalarea")
state = usstmap$names
set.seed(1234)
var.2010 = rnorm (length (state), 50, 10)
var.2011 = var.2010*1.1 + rnorm (length (state), 0, 5)
var.2012 = var.2011*0.98 + rnorm (length (state), 0, 10)
var.2013 = var.2011*1.2 + rnorm (length (state), 0, 15)
usst <- data.frame (county, var.2010, var.2011, var.2012, var.2013)

mapplot(state ~ var.2010 + var.2011 + var.2012 + var.2013, data = usst, map = map("state",    plot = FALSE, fill = TRUE), colramp = colorRampPalette(c("green", "purple")))

RG#99: cloud 3D bars with heatmap

require(lattice)

require(latticeExtra)

data(VADeaths)

cloud(VADeaths, panel.3d.cloud = panel.3dbars,
      xbase = 0.4, ybase = 0.4, zlim = c(0, max(VADeaths)),
      scales = list(arrows = FALSE, just = "right"), xlab = NULL, ylab = NULL,
      col.facet = level.colors(VADeaths, at = do.breaks(range(VADeaths), 20),
                               col.regions = cm.colors,
                               colors = TRUE),
      colorkey = list(col = cm.colors, at = do.breaks(range(VADeaths), 20)),
      screen = list(z = 40, x = -30))

RG#98: Horizon plot (time series data)

require(latticeExtra)

#example data 
set.seed(123)
mydat <- ts(matrix(cumsum(rnorm(150 * 10)), ncol = 10))
colnames(mydat) <- paste("TS", letters[1:10], sep = "-")

#simple line plot
xyplot(mydat, scales = list(y = "same"))


# panel with different origin and scale:
horizonplot(mydat, layout = c(1,12), colorkey = TRUE) +  

 layer(panel.scaleArrow(x = 0.99, digits = 1, col = "grey",

                         srt = 90, cex = 0.7)) +
  layer(lim <- current.panel.limits(),
    panel.text(lim$x[1], lim$y[1], round(lim$y[1],1), font = 2,

        cex = 0.7, adj = c(-0.5,-0.5), col = "blue")) 

RG#97: Error bar plot with significance (line connecting) - publication purpose

myd <- data.frame (X = c(1:12,1:12),
                   Y = c(8, 12, 13, 18,  22, 16, 24, 29,  34, 15, 8, 6,
                         9, 10, 12, 18, 26, 28, 28, 30, 20, 10, 9, 9),
                   group = rep (c("A-group", "B-group"), each = 12),
                   error = rep (c(2.5, 3.0), each = 12))

plot1 = ggplot(data = myd, aes(x=X, y=Y, fill=group, width=0.8) ) +
     geom_errorbar(aes(ymin=Y, ymax=Y+error, width = 0.2), position=position_dodge(width=0.8)) +
     geom_bar(stat="identity", position=position_dodge(width=0.8)) +
     geom_bar(stat="identity", position=position_dodge(width=0.8), colour="black", legend=FALSE) +
     scale_fill_manual(values=c("grey70", "white")) +
     scale_x_discrete("X", limits=c(1:12)) +
     scale_y_continuous("Y (units)", expand=c(0,0), limits = c(0, 40), breaks=seq(0, 40, by=5)) + ggtitle ("My nice plot") +
     theme_bw() +
    theme( plot.title = element_text(face="bold", size=14),
          axis.title.x = element_text(face="bold", size=12),
          axis.title.y = element_text(face="bold", size=12, angle=90),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.text.y=element_text(angle=90, hjust=0.5),
          legend.title = element_blank(),
          legend.position = c(0.85,0.85),
          legend.key.size = unit(1.5, "lines"),
          legend.key = element_rect()
     )

print(plot1)

   # connecting segments
   plot1 + geom_segment(aes(x=5.7, y=16, xend=5.7, yend=28.5), col = "gray80")+
   geom_segment(aes(x=5.7, y=28.5, xend=6.1, yend=28.5), col = "gray80")      +
   geom_segment(aes(x=6.1, y=28.5, xend=6.1, yend=28), col = "gray80") +
      annotate("text", x=5.7, y=31.2, label="***")

RG#96: Basic point and line graph with error bars (publication purpose)

myd <- data.frame (X = c(1:12,1:12),
                   Y = c(8, 12, 13, 18,  22, 16, 24, 29,  34, 15, 8, 6,
                         9, 10, 12, 18, 26, 28, 28, 30, 20, 10, 9, 9),
                   group = rep (c("A-group", "B-group"), each = 12),
                   error = rep (c(2.5, 3.0), each = 12))
                   
require(ggplot2)
require(grid)
# line and point plot
f1 = ggplot(data = myd, aes(x = X, y = Y, group = group) )  # lesion becomes a classifying factor
f2 <- f1 + geom_errorbar(aes(ymin = Y - error, ymax = Y + error), width=0.3) +
geom_line() + geom_point(aes(shape=group, fill=group), size=5)

 f3 <- f2 +  scale_x_continuous("X (units)", breaks=1:12) +
     scale_y_continuous("Y (units)", limits = c(0, 40), breaks=seq(0, 40, by = 5)) +
     scale_shape_manual(values=c(24,21)) +
     scale_fill_manual(values=c("white","black")) +
     stat_abline(intercept=0, slope=0, linetype="dotted") +
     annotate("text", x=11, y=10, label="X") +
     theme_bw()

   optns <- theme (
          plot.title = element_text(face="bold", size=14),
          axis.title.x = element_text(face="bold", size=12),
          axis.title.y = element_text(face="bold", size=12, angle=90),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = c(0.2,0.8),
          legend.title = element_blank(),
          legend.text = element_text(size=12),
          legend.key.size = unit(1.5, "lines"),
          legend.key = element_blank()
     )
f3 +  ggtitle ( "MY awsome plot for publication") + optns

RG#95: Interactive Biplot

# data

set.seed(1234)

P <- vector()
DF <- as.data.frame(matrix(rep(NA, 100), nrow=10))
names(DF) <- c(paste("M",1:10, sep=""))
for(i in 1:10) {
DF[,i] <- rnorm(10, 10, 3)
}

rownames (DF) <- paste("O", 1:10, sep = "")

require(BiplotGUI)

Biplots(Data = DF, groups = rep(1, nrow(DF)),PointLabels = rownames(DF),
AxisLabels = colnames(DF))

# you can work in the interactive menu window

RG#94: Hive plot (diagram)

require(HiveR)

# test random data, nodes = 200, edge 200   

td1 <- ranHiveData(nx = 3, nn = 200, ne = 200)

plotHive(td1, bkgnd = "cornsilk") 

RG#94: Arch Diagram

# 'devtools' need to be installed if not installed 
install.packages("devtools") 

# load devtools
library(devtools)

# install 'arcdiagram'
install_github('arcdiagram',  username='gastonstat')

 library(arcdiagram)

# star graph with 20 nodes
star_graph = graph.star(20, mode="out")

# extract edgelist
staredges = get.edgelist(star_graph)

# default arc diagram
arcplot(staredges)

 # different ordering, arc widths, arc colors, and node sizes
set.seed(1234)
orderv <- sample (1:20)
labelv <- paste("nd",1:20,sep="-")
archlinewd <- 4*runif(20,.5,2)
colarcs <- heat.colors (19)
 nodesize <- runif(20,1,3)
arcplot(staredges, ordering=orderv, labels= labelv,
   lwd.arcs= archlinewd, col.arcs= colarcs,
   show.nodes=TRUE, pch.nodes=21, cex.nodes= nodesize,
   col.nodes="lightseagreen", bg.nodes="midnightblue", lwd.nodes= 2)

RG#93: Add countour or heat map plot to XY scatter plot

# data
set.seed(1234)
n <- 10000
X = rnorm (n, 10, 4)
Y = X*1.5 + rnorm (n, 0, 8)

## colour brewing
library(RColorBrewer)
g = 11
my.cols <- rev(brewer.pal(g, "RdYlBu"))
#compute 2D kernel density

# kernel density using MASS 
library(MASS)
z <- kde2d(X, Y, n=50)
plot(X, Y, xlab="X", ylab="Y", pch=19, cex=.3, col = "gray60")
contour(z, drawlabels=FALSE, nlevels=g, col=my.cols, add=TRUE, lwd = 2)
abline(h=mean(Y), v=mean(X), lwd=2, col = "black")
legend("topleft", paste("r=", round(cor(X, Y),2)), bty="n")

## estimate the z counts
prob <- c(.99, .95, .90, .8, .5, .1, 0.05)
dx <- diff(z$x[1:2])
dy <- diff(z$y[1:2])
sz <- sort(z$z)

c1 <- cumsum(sz) * dx * dy

levels <- sapply(prob, function(x) {
              approx(c1, sz, xout = 1 - x)$y })

plot(X,Y, col = "gray80", pch = 19, cex = 0.3)
contour(z, levels= round (levels,7), add=T, col = "red")

# smooth scatter
require(KernSmooth)
smoothScatter(X, Y, nrpoints=.3*n, colramp=colorRampPalette(my.cols), pch=19, cex=.3, col = "green1")




 

RG#91: Plot bar or pie chart over world map using rworldmap package

require (rworldmap)
require(rworldxtra)

# get world map 

plot(getMap(resolution = "high" ))

##getting example data

dataf <- getMap()@data 

mapBars( dataf, nameX="LON", nameY="LAT" , nameZs=c('GDP_MD_EST',

'GDP_MD_EST','GDP_MD_EST') , mapRegion='asia' , symbolSize=2  ,

 barOrient = 'horiz' )

mapBars( dataf, nameX="LON", nameY="LAT" , nameZs=c('GDP_MD_EST','GDP_MD_EST',

'GDP_MD_EST') , mapRegion='asia' , symbolSize=3  , barOrient = 'vert' ,  

oceanCol = "blue1", landCol = "lightgreen")

mapPies( dataf,nameX="LON", nameY="LAT", nameZs=c('GDP_MD_EST','GDP_MD_EST',

'GDP_MD_EST','GDP_MD_EST'),mapRegion='asia', oceanCol = "lightseagreen",

 landCol = "gray50")

RG#90: fluctutation diagram: graphical representation of a contingency table

set.seed (1234)
myd <- data.frame (x1 = rnorm (1000, 15, 5), x2 = sample (c("A", "B", "C"), 1000, replace = TRUE), x3 = sample (c(1,2,2), 1000, replace = TRUE))


# fluctuation plot
require (ggplot2)

ggfluctuation(table(myd$x2, myd$x3))+theme_bw()

ggfluctuation(table(myd$x2, myd$x3), type="colour")+theme_bw()

ggfluctuation(table(cut (myd$x1,4), myd$x3), type="colour")+theme_bw()

ggfluctuation(table(cut (myd$x1,5), myd$x3))+theme_bw()

RG#89: Raster in R

 require(grid)
grid.raster(matrix(colors()[1:600], ncol=20))

 

set.seed(1234)

mt <- matrix (sample(c("red","red1", "yellow", "purple", "green1", "green4", "blue"), 10000, replace = TRUE), ncol = 100)

grid.raster(mt)

rgb.palette <- colorRampPalette(c("red", "orange", "blue"),
                                space = "rgb") 

bg.palette <- colorRampPalette(c("blue", "green"), space = "rgb")

gr.palette <-  colorRampPalette(c("green", "red"),
                                space = "rgb")

colrs <- matrix(c(rgb.palette(20),bg.palette(20),gr.palette(20)),   nrow=6, ncol=10)

grid.newpage()

 grid.raster(colrs)
  grid.raster(colrs, interpolate=FALSE)

# raster in ggplot2

require(ggplot2)

# Generate data
funp <- function (n,r=2) {
 xv <- seq(-r*pi, r*pi, len=n)
 df1 <- expand.grid(x=xv, y=xv)
 df1$r <- sqrt(df1$x^2 + df1$y^2)
 df1$z <- cos(df1$r^2)*exp(-df1$r/6)
 df1
}
qplot(x, y, data = funp(1000), fill = z, geom = "raster") + theme_bw()