Friday, May 3, 2013

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)





No comments:

Post a Comment

Note: Only a member of this blog may post a comment.