R graph gallery: May 2013

Wednesday, May 22, 2013

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)

Friday, May 3, 2013

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)

Thursday, May 2, 2013

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)

Wednesday, May 1, 2013

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="***")

Labels

Wednesday, May 22, 2013

Friday, May 3, 2013

Thursday, May 2, 2013

Wednesday, May 1, 2013