Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

ggplot2 - fail to apply color with scale_fill_manual inside a loop [closed]

Tags:

loops

r

ggplot2

I'm running a loop to get, at each sub setting of my data set, a map and apply a given palette (and respective legend) accordingly.

People tend to dislike the use of for() loops and maximize vectorization of their approaches. I don't know the best way to vectorize processes with this particular data set.

In this particular case, I'm handling a relatively large data set (distribution species Atlas) that is particularly complex since different methodologies were used and different options must be passed for each species, considering a particular season, different set of observations, etc. Species may be present at one season and missed into another (They may be a breeder, a resident or a migrant). Maps should be created for all cases (seasons), empty when absent. Additional data (besides those from field work) may be available and used. Map Legend must accommodate all variations, besides presenting variable in interest (abundances) in a custom discrete scale.

By running a loop I feel (to my limited expertise) I can more easily retain and control the several needed objects, while stepping into the flux I created to produce the pieces of interest and finally create sets of species distributions maps.

My Problem is that I'm storing each resulting ggplot in a list() object. Each species at each season will be stored in a list. The issue I'm facing is related to scale_fill_manual when used inside a loop.

The behavior is strange since I get the maps done but with colors applied only to the last ggplot output. Nonetheless all values still being correctly identified in the legend.

to exemplify:

Packages

if (!require(ggplot2)) install.packages("ggplot2",
    repos = "http://cran.r-project.org"); library(ggplot2)
if (!require(grid)) install.packages("grid",
    repos = "http://cran.r-project.org"); library(grid)
if (!require(RColorBrewer)) install.packages("RColorBrewer",
    repos = "http://cran.r-project.org"); library(RColorBrewer)
if (!require(reshape)) install.packages("reshape",
    repos = "http://cran.r-project.org"); library(reshape)

A simple example first

#Create a list of colors to be used with scale_manual
palette.l <- list()
palette.l[[1]] <- c('red', 'blue', 'green')
palette.l[[2]] <- c('pink', 'blue', 'yellow')
# Store each ggplot in a list object
plot.l <- list()
#Loop it
for(i in 1:2){
  plot.l[[i]] <- qplot(mpg, wt, data = mtcars, colour = factor(cyl)) +
    scale_colour_manual(values = palette.l[[i]])
}

In my session plot.l[1] will be painted with colors from palette.l[2].

My particular case

Functions

Arrange Plots

ArrangeGraph <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
  dots <- list(...)
  n <- length(dots)
  if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
  if(is.null(nrow)) { nrow = ceiling(n/ncol)}
  if(is.null(ncol)) { ncol = ceiling(n/nrow)}
  ## NOTE see n2mfrow in grDevices for possible alternative
  grid.newpage()
  pushViewport(viewport(layout=grid.layout(nrow,ncol)))
  ii.p <- 1
  for(ii.row in seq(1, nrow)) {
    ii.table.row <- ii.row
    if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
    for(ii.col in seq(1, ncol)) {
      ii.table <- ii.p
      if(ii.p > n) break
      print(dots[[ii.table]], vp=VPortLayout(ii.table.row, ii.col))
      ii.p <- ii.p + 1
    }
  }
}

ViewPort

VPortLayout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)

Species data sets

bd.aves.1 <- structure(list(quad = c("K113", "K114", "K114", "K114", "K114",...
due to limited body character number limit, please download entire code from
https://docs.google.com/open?id=0BxSZDr4eTnb9R09iSndzZjBMS28

Species list

list.esp.1 <- c("Sylv mela", "Saxi rube","Ocea leuc")#
# download from the above link

Some taxonomy and other data

txcon.1 <- structure(list(id = c(156L, 359L, 387L), grupo = c("Aves", "Aves",# 
# download from the above link

Seasons

kSeason.1 <- c("Inverno", "Primavera", "Outono")

A Sample Grid

grid500.df.1 <- structure(list(id = c("K113", "K113", "K113", "K113", "K113",#... 
# download from the above link

Additional Mapping elements

Shoreline

coastline.df.1 <- structure(list(long = c(182554.963670234, 180518, 178865.39,#...
# download from the above link

Labels placement adjustments

kFacx1 <- c(9000, -13000, -10000, -12000)

The R Code

for(i in listsp.1) { # LOOP 1 - Species
  # Set up objects 
  sist.i <- list() # Sistematic observations
  nsist.i <- list() # Non-Sistematic observations
  breaks.nind.1 <- list() # Breaks on abundances
  ## Grid and merged dataframe
  spij.1 <- list() # Stores a dataframe for sp i at season j
  ## Palette build
  classes.1 <- list()
  cllevels.1 <- list()
  palette.nind.1 <- list() # Color palette
  ## Maps
  grid500ij.1 <- list() # Grid for species i at season j
  map.dist.ij.1 <- NULL
  for(j in 1:length(kSeason.1)) { # LOOP 2 - Seasons
    # j assume each season: Inverno, Primavera, Outono
    # Sistematic occurences ===================================================
    sist.i.tmp <- nrow(subset(bd.aves.1, esp == i & cod_tipo %in% sistematica &
      periodo == kSeason.1[j]))
    if (sist.i.tmp!= 0) { # There is sistematic entries, Then:
      sist.i[[j]]<- ddply(subset(bd.aves.1,
                                 esp == i & cod_tipo %in% sistematica & 
                                   periodo == kSeason.1[j]),
                          .(periodo, quad), summarise, nind = sum(n_ind),
                          codnid = max(cod_nidi))
    } else { # No Sistematic entries, Then: 
      sist.i[[j]] <- data.frame('quad' = NA, 'periodo' = NA, 'nind' = NA, 
                                'codnid' = NA, stringsAsFactors = F)
    }
    # Additional Entries (RS1) e other non-sistematic entries  (biblio) =======
    nsist.tmp.i = nrow(subset(bd.aves.1, esp == i & !cod_tipo %in% sistematica &
      periodo == kSeason.1[j]))
    if (nsist.tmp.i != 0) { # RS1 and biblio entries, Then:
      nsist.i[[j]] <- subset(bd.aves.1,
                             esp == i & !cod_tipo %in% sistematica &
                               periodo == kSeason.1[j] & 
                               !quad %in% if (nrow(sist.i[[j]]) != 0) {
                                            subset(sist.i[[j]],
                                                   select = quad)$quad
                                          } else NA,
                             select = c(quad, periodo, cod_tipo, cod_nidi)
                             )
      names(nsist.i[[j]])[4] <- 'codnid'
    } else { # No RS1 and biblio entries, Then:      
        nsist.i[[j]] = data.frame('quad' = NA, 'periodo' = NA, 'cod_tipo' = NA,
                                'codnid' = NA, stringsAsFactors = F)
      }
    # Quantile breaks =========================================================
    if (!is.na(sist.i[[j]]$nind[1])) {
      breaks.nind.1[[j]] <- c(0,
                            unique(
                              ceiling(
                                quantile(unique(
                                  subset(sist.i[[j]], is.na(nind) == F)$nind), 
                                         q = seq(0, 1, by = 0.25)))))
    } else {
        breaks.nind.1[[j]] <- 0
      }
    # =========================================================================
    # Build Species dataframe and merge to grid
    # =========================================================================
    if (!is.na(sist.i[[j]]$nind[1])) { # There are Sistematic entries, Then:
      spij.1[[j]] <- merge(unique(subset(grid500df.1, select = id)),
                         sist.i[[j]],
                         by.x = 'id', by.y = 'quad', all.x = T)
      # Adjust abundances when equals to NA ===================================
      spij.1[[j]]$nind[is.na(spij.1[[j]]$nind) == T] <- 0
      # Break abundances to create discrete variable ==========================
      spij.1[[j]]$cln <- if (length(breaks.nind.1[[j]]) > 2) {
        cut(spij.1[[j]]$nind, breaks = breaks.nind.1[[j]], 
            include.lowest = T, right = F)
        } else {
            cut2(spij.1[[j]]$nind, g = 2)
          }
      # Variable Abundance ====================================================
      classes.1[[j]] = nlevels(spij.1[[j]]$cln)
      cllevels.1[[j]] = levels(spij.1[[j]]$cln)
      # Color Palette for abundances - isolated Zero class (color #FFFFFF) ====
      if (length(breaks.nind.1[[j]]) > 2) {
        palette.nind.1[[paste(kSeason.1[j])]] = c("#FFFFFF", brewer.pal(length(
          cllevels.1[[j]]) - 1, "YlOrRd"))
      } else {
          palette.nind.1[[paste(kSeason.1[j])]] = c(
            "#FFFFFF",  brewer.pal(3, "YlOrRd"))[1:classes.1[[j]]]
        }
        names(palette.nind.1[[paste(kSeason.1[j])]])[1 : length(
          palette.nind.1[[paste(kSeason.1[j])]])] <- cllevels.1[[j]]
      # Add RS1 and bilbio values to palette ==================================
      palette.nind.1[[paste(kSeason.1[j])]][length(
        palette.nind.1[[paste(kSeason.1[j])]]) + 1] <- '#CCC5AF'
      names(palette.nind.1[[paste(kSeason.1[j])]])[length(
        palette.nind.1[[paste(kSeason.1[j])]])] <- 'Suplementar'
      palette.nind.1[[paste(kSeason.1[j])]][length(
        palette.nind.1[[paste(kSeason.1[j])]]) + 1] <- '#ADCCD7'
      names(palette.nind.1[[paste(kSeason.1[j])]])[length(
        palette.nind.1[[paste(kSeason.1[j])]])] <- 'Bibliografia'
      # Merge species i dataframe to grid map =================================
      grid500ij.1[[j]] <- subset(grid500df.1, select = c(id, long, lat, order))
      grid500ij.1[[j]]$cln = merge(grid500ij.1[[j]],
                                 spij.1[[j]],
                                 by.x = 'id', by.y = 'id', all.x = T)$cln
      # Adjust factor levels of cln variable - Non-Sistematic data ============
      levels(grid500ij.1[[j]]$cln) <- c(levels(grid500ij.1[[j]]$cln), 'Suplementar',
                                      'Bibliografia')
      if (!is.na(nsist.i[[j]]$quad[1])) {
        grid500ij.1[[j]]$cln[grid500ij.1[[j]]$id %in% subset(
          nsist.i[[j]], cod_tipo == 'RS1', select = quad)$quad] <- 'Suplementar'
        grid500ij.1[[j]]$cln[grid500ij.1[[j]]$id %in% subset(
          nsist.i[[j]], cod_tipo == 'biblio', select = quad)$quad] <- 'Bibliografia'
      }
    } else { # No Sistematic entries, Then:
        if (!is.na(nsist.i[[j]]$quad[1])) { # RS1 or Biblio entries, Then:
          grid500ij.1[[j]] <- grid500df
          grid500ij.1[[j]]$cln <- '0'
          grid500ij.1[[j]]$cln <- factor(grid500ij.1[[j]]$cln)
          levels(grid500ij.1[[j]]$cln) <- c(levels(grid500ij.1[[j]]$cln),
                                          'Suplementar', 'Bibliografia')
          grid500ij.1[[j]]$cln[grid500ij.1[[j]]$id %in% subset(
            nsist.i[[j]], cod_tipo == 'RS1', 
            select = quad)$quad] <- 'Suplementar'
          grid500ij.1[[j]]$cln[grid500ij.1[[j]]$id %in% subset(
            nsist.i[[j]],cod_tipo == 'biblio', 
            select = quad)$quad] <- 'Bibliografia'
        } else { # No entries, Then:
            grid500ij.1[[j]] <- grid500df
            grid500ij.1[[j]]$cln <- '0' 
            grid500ij.1[[j]]$cln <- factor(grid500ij.1[[j]]$cln)
            levels(grid500ij.1[[j]]$cln) <- c(levels(grid500ij.1[[j]]$cln),
                                            'Suplementar', 'Bibliografia')      
          }
      } # End of Species dataframe build
    # Distribution Map for  species i at season j =============================    
    if (!is.na(sist.i[[j]]$nind[1])) { # There is sistematic entries, Then:
      map.dist.ij.1[[paste(kSeason.1[j])]] <- ggplot(grid500ij.1[[j]],
                                                  aes(x = long, y = lat)) +
        geom_polygon(aes(group = id, fill = cln), colour = 'grey80') +
        coord_equal() +
        scale_x_continuous(limits = c(100000, 180000)) +
        scale_y_continuous(limits = c(-4000, 50000)) +
        scale_fill_manual(
          name = paste("LEGEND",
                       '\nSeason: ', kSeason.1[j],
                       '\n% of Occupied Cells : ',
                         sprintf("%.1f%%", (length(unique(
                           grid500ij.1[[j]]$id[grid500ij.1[[j]]$cln != levels(
                           grid500ij.1[[j]]$cln)[1]]))/12)*100), # percent 
                        sep = ""
                       ),
          # Set Limits
          limits = names(palette.nind.1[[j]])[2:length(names(palette.nind.1[[j]]))],
          values = palette.nind.1[[j]][2:length(names(palette.nind.1[[j]]))],
          drop = F) +
          opts(
            panel.background = theme_rect(),
            panel.grid.major = theme_blank(),
            panel.grid.minor = theme_blank(),
            axis.ticks = theme_blank(),
            title = txcon.1$especie[txcon.1$esp == i],
            plot.title = theme_text(size = 10, face = 'italic'),
            axis.text.x = theme_blank(),
            axis.text.y = theme_blank(),
            axis.title.x = theme_blank(),
            axis.title.y = theme_blank(),
            legend.title = theme_text(hjust = 0,size = 10.5),
            legend.text = theme_text(hjust = -0.2, size = 10.5)
          ) +
          # Shoreline
          geom_path(inherit.aes = F, aes(x = long, y = lat),
                    data = coastline.df.1, colour = "#997744") +
          # Add localities
          geom_point(inherit.aes = F, aes(x = x, y = y),  colour = 'grey20',
                     data = localidades, size = 2) +
          # Add labels
          geom_text(inherit.aes = F, aes(x = x, y = y, label = c('Burgau',
                                                                 'Sagres')),
                    colour = "black",
                    data = data.frame(x = c(142817 + kFacx1[1], 127337 + kFacx1[4]),
                                      y = c(11886, 3962), size = 3))
    } else { # NO sistematic entries,then:
        map.dist.ij.1[[paste(kSeason.1[j])]] <- ggplot(grid500ij.1[[j]],
                                                    aes(x = long, y = lat)) +
          geom_polygon(aes.inherit = F, aes(group = id, fill = cln),
                       colour = 'grey80') +
          #scale_color_manual(values = kCorLimiteGrid) +
          coord_equal() +
          scale_x_continuous(limits = c(100000, 40000)) +
          scale_y_continuous(limits = c(-4000, 180000)) +
          scale_fill_manual(
            name = paste('LEGENDA',
                         '\nSeason: ', kSeason.1[j],
                         '\n% of Occupied Cells :',
                         sprintf("%.1f%%", (length(unique(
                           grid500ij.1[[j]]$id[grid500ij.1[[j]]$cln != levels(
                           grid500ij.1[[j]]$cln)[1]]))/12 * 100)), # percent 
                         sep = ''),
            limits = names(kPaletaNsis)[2:length(names(kPaletaNsis))],
            values = kPaletaNsis[2:length(names(kPaletaNsis))],
            drop = F) +
            opts(
              panel.background = theme_rect(),
              panel.grid.major = theme_blank(),
              panel.grid.minor = theme_blank(),
              title = txcon.1$especie[txcon.1$esp == i],
              plot.title = theme_text(size = 10, face = 'italic'),
              axis.ticks = theme_blank(),
              axis.text.x = theme_blank(),
              axis.text.y = theme_blank(),
              axis.title.x = theme_blank(),
              axis.title.y = theme_blank(),
              legend.title = theme_text(hjust = 0,size = 10.5),
              legend.text = theme_text(hjust = -0.2, size = 10.5)
            ) +
            # Add Shoreline
            geom_path(inherit.aes = F, data = coastline.df.1,
                      aes(x = long, y = lat),
                      colour = "#997744") +
            # Add Localities
            geom_point(inherit.aes = F, aes(x = x, y = y),
                       colour = 'grey20',
                       data = localidades, size = 2) +
            # Add labels
            geom_text(inherit.aes = F, aes(x = x, y = y,
                                           label = c('Burgau', 'Sagres')),
                      colour = "black",
                      data = data.frame(x = c(142817 + kFacx1[1],
                                              127337 + kFacx1[4],),
                                        y = c(11886, 3962)),
                      size = 3)
      } # End of Distribution map building for esp i and j seasons
  } # Fim do LOOP 2: j Estacoes
  # Print Maps
  png(file = paste('panel_species',i,'.png', sep = ''), res = 96, 
      width = 800, height = 800)
  ArrangeGraph(map.dist.ij.1[[paste(kSeason.1[3])]],
               map.dist.ij.1[[paste(kSeason.1[2])]],
               map.dist.ij.1[[paste(kSeason.1[1])]],
               ncol = 2, nrow = 2)
  dev.off()
  graphics.off()
} # End of LOOP 1

map.dist.ij.1[[paste(kSeason.1[3])]] is the only with color palette applied to polygons, but the legend items is well defined for each j map.

Output using R Code

enter image description here

As we see, Legends are OK but not colored.

Hope not missing anything. Sorry for some lost Portuguese terminology.

like image 264
Paulo E. Cardoso Avatar asked Nov 13 '22 23:11

Paulo E. Cardoso


1 Answers

Honestly, I have not looked much at your code for your specific problem--a bit too much to wade through!--but for your demo example, adding print(plot.l[[i]]) in your loop.

#Create a list of colors to be used with scale_manual
palette.l <- list()
palette.l[[1]] <- c('red', 'blue', 'green')
palette.l[[2]] <- c('pink', 'blue', 'yellow')

# Store each ggplot in a list object
plot.l <- list()

# Loop it
for(i in 1:2) {
  plot.l[[i]] <- qplot(mpg, wt, data = mtcars, colour = factor(cyl)) +
    scale_colour_manual(values = palette.l[[i]])
  print(plot.l[[i]]) ### Added to your loop
}

In the case of your minimal example, though, this also works (without first having to create an empty list to store your plots) and I think it at least looks a lot cleaner. I'm not sure if something similar can be adapted to suit your larger scenario.

#Create a list of colors to be used with scale_manual
palette.l <- list(c('red', 'blue', 'green'),
                  c('pink', 'blue', 'yellow'))

p <- qplot(mpg, wt, data = mtcars, colour = factor(cyl))

# Use lapply and "force" to get your plots in a list    
plot.l <- lapply(palette.l, 
                 function(x) { 
                   force(x)
                   p + scale_color_manual(values = x) 
                 })
like image 78
A5C1D2H2I1M1N2O1R2T1 Avatar answered Nov 15 '22 13:11

A5C1D2H2I1M1N2O1R2T1