Usa R para recrear la gráfica de contorno hecha en Igor

plot ggplot2

Question

Esta gráfica de contorno, realizada con el programa Igor, es popular en estudios de química atmosférica y contaminación:

Intento recrearlo con R para un amigo que quiere dejar de usar Igor, y no podemos entenderlo. Aquí está el conjunto de datos (los mismos datos utilizados para hacer la trama con Igor), y esto es lo que tengo hasta ahora para hacer la trama con R:

# read in the data
dat <- read.csv("contour_plot_data.csv")

# focus on the untransformed values
dat <- dat[, 1:108]

# get Diameter value from col names
Diameter <- as.numeric(gsub("X", "", names(dat)[-1]))

# interpolate between the Diameter values for a smoother contour,
# a seperate interpolation for each row (date value)
# this takes a moment or two...

interp <-  seq(min(Diameter), max(Diameter), 0.2)
dat_interp <- data.frame(matrix(0, ncol =  length(interp), nrow = nrow(dat)))
for(i in 1:nrow(dat)){
  # get the values from row i
  vec <- unlist(dat[i, 2:108], use.names = FALSE)
  # compute loess interpolations
  lo <- loess(vec ~ Diameter)
  # predict interpolated values
  pr <- predict(lo, newdata = data.frame(Diameter = interp))
  # store in a data frame
  df <- data.frame(ct = unname(pr), Diameter = interp)
  # add as new row to new data frame
  dat_interp[i, ] <- df$ct
  print(i) # so we can see that it''s working
}

# add date col and col names to the interpolated data
names(dat_interp) <- interp
dat_interp$date <- as.character(dat$Time)

# melt data into long format
# see http://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
library(tidyr)
gather_cols <- interp
dat_long <- gather_(dat_interp, "Diameter", "dN_dlogDp", gather_cols)

# we want diameter as a numeric
dat_long$Diameter <- as.numeric(as.character(dat_long$Diameter))
# we want date as a date format
x <-  as.character(dat_long$date)
date_ <- as.Date(x, format = "%d/%m/%Y") 
time_ <- gsub(" ", "", substr(x, nchar(x) - 4, nchar(x)))
dat_long$date_time <- as.POSIXct(paste0(date_, " ", time_))

# The Igor plot seems to use log dN_dlogDp values, so let''s get those
dat_long$dN_dlogDp_log <- log10(dat_long$dN_dlogDp)
dat_long$dN_dlogDp_log <- ifelse(dat_long$dN_dlogDp_log == "NaN", 0, dat_long$dN_dlogDp_log)

# get on with plottong...
library(ggplot2)   
library(scales)

labels_breaks <- seq(0, max(Diameter), 100)
mytheme <- theme_bw(base_size = 14) +  theme(aspect.ratio = 1/4) 
ggplot(dat_long, aes(y = Diameter, x = date_time,  fill=dN_dlogDp_log)) +
  geom_raster(interpolate = TRUE)  +
  scale_fill_gradientn(name=expression(log(dN/dlogD[p])), colours = rainbow(7)) +
  scale_y_continuous(expand = c(0,0), breaks = labels_breaks ) +
  scale_x_datetime(expand = c(0,0), breaks = date_breaks("12 hours")) +
    ylab("Diameter (nm)") +
  xlab("Date and time") +
  mytheme

Mi argumento podría ser un poco más fino con etiquetas y marcas, etc. Sin embargo, mi pregunta principal es por qué mi relleno de contorno se ve tan diferente de la trama de Igor. La escala parece invertida, y la interpolación se ve muy diferente.

¿Cómo puedo hacer que mi trama se parezca más a la trama de Igor?

Tenga en cuenta que estas otras preguntas mías están estrechamente relacionadas con la tarea de recrear esta trama:

• interpolación geom_raster con escala de registro
• Parcela de densidad 2d para categorías

Y después de hacer esta pregunta, he mantenido una esencia actualizada de código R que combina detalles de las respuestas a estas preguntas, y replica con éxito estas gráficas (salida de ejemplo incluida en la esencia). Esa esencia está aquí: https://gist.github.com/benmarwick/9a54cbd325149a8ff405 .

ACTUALIZAR Ahora he hecho un paquete que producirá estas parcelas: https://github.com/benmarwick/smps

Answer 1

Puedo acercarme mucho más a la trama de Igor utilizando akima::interp lugar de loess para la interpolación:

# read in the data
dat <- read.csv("contour_plot_data.csv")

# focus on the untransformed values
dat <- dat[, 1:108]

# get Diameter value from col names
Diameter <- as.numeric(gsub("X", "", names(dat)[-1]))

# melt data into long format
# see http://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
library(tidyr)
dat_long <- gather(dat, "Diameter", "dN_dlogDp", 2:108)

# we want diameter as a numeric
dat_long$Diameter <- as.numeric(gsub("X", "", dat_long$Diameter ))
# we want time as a date-formatted variable
x <-  as.character(dat_long$Time)
date_ <- as.Date(x, format = "%d/%m/%Y") 
time_ <- gsub(" ", "", substr(x, nchar(x) - 4, nchar(x)))
dat_long$Time <- as.POSIXct(paste0(date_, " ", time_))

# The Igor plot seems to use log dN_dlogDp values, so let''s get those
dat_long$dN_dlogDp_log <- log10(dat_long$dN_dlogDp)
dat_long$dN_dlogDp_log <- ifelse(dat_long$dN_dlogDp_log == "NaN" |
                                   dat_long$dN_dlogDp_log == "-Inf"  , 0, dat_long$dN_dlogDp_log)


# interpolate between the  values for a smoother contour
# this takes a moment or two...

library(akima)
xo <- with(dat_long, seq(min(Time), max(Time), 120))
yo <- with(dat_long, seq(min(Diameter), max(Diameter), 0.5))
dat_interp <- with(dat_long, interp(Time, Diameter, dN_dlogDp_log, xo = xo, yo = yo) )

# get on with plotting...

# make into a data frame for ggplot
dat_interp_df <-  data.frame(matrix(data = dat_interp$z, ncol = length(dat_interp$y), nrow = length(dat_interp$x)))
names(dat_interp_df) <- dat_interp$y
dat_interp_df$Time <- as.POSIXct(dat_interp$x, origin = "1970-01-01")

# wide to long
dat_interp_df_long <- gather(dat_interp_df, "Diameter", "dN_dlogDp_log", 1:(ncol(dat_interp_df)-1))
dat_interp_df_long$Diameter <- as.numeric(as.character(dat_interp_df_long$Diameter))

# plot
library(ggplot2) 
library(scales)
y_labels_breaks <- seq(0, max(Diameter), 100)
ggplot(dat_interp_df_long, aes(y = Diameter, x = Time,  fill = dN_dlogDp_log)) +
  geom_raster(interpolate = TRUE)  +
  scale_fill_gradientn(name=expression(log(dN/dlogD[p])), colours = rev(rainbow(50))) +
  scale_y_continuous(expand = c(0,0), breaks = y_labels_breaks ) +
  scale_x_datetime(expand = c(0,0), breaks = date_breaks("1 day"))

Pero todavía hay una gran diferencia en el mapeo de colores, con la trama Igor que tiene bandas anchas con límites definidos, y mi trama tiene menos bandas de color y límites difusos entre ellas. Así que supongo que no tengo el método de interpolación que usa la trama Igor.

ACTUALIZADO después de experimentar con un grupo de rampas de color, encontré una colorRamps::blue2green2red bastante buena en colorRamps::blue2green2red . También puse un poco de esfuerzo aquí en elegantes tickmarks:

# plot
library(ggplot2) 
library(scales) # for date_breaks
library(colorRamps) # for blue2green2red

# function for minor tick marks
every_nth <- function(x, nth, empty = TRUE, inverse = FALSE) 
{
  if (!inverse) {
    if(empty) {
      x[1:nth == 1] <- ""
      x
    } else {
      x[1:nth != 1]
    }
  } else {
    if(empty) {
      x[1:nth != 1] <- ""
      x
    } else {
      x[1:nth == 1]
    }
  }
}

# add tick marks every two hours
start_date <- min(dat_interp_df_long$Time)
end_date <-  max(dat_interp_df_long$Time)
date_breaks_2h <-  seq(from = start_date, to = end_date, by = "2 hours")
date_breaks_1_day <- seq(from = start_date, to = end_date, by = "1 day")
multiple <- length(date_breaks_2h) / length(date_breaks_1_day)

insert_minor <- function(major_labs, n_minor) {labs <- 
  c( sapply( major_labs, function(x) c(x, rep("", multiple) ) ) )
labs[1:(length(labs)-n_minor)]}


y_labels_breaks <- seq(0, max(Diameter), 100)
mytheme <- theme_bw(base_size = 14) +  theme(aspect.ratio = 1/5)
ggplot(dat_interp_df_long, aes(y = Diameter, x = Time,  fill = dN_dlogDp_log)) +
  geom_raster(interpolate = TRUE)  +
  scale_fill_gradientn(name=expression(log(dN/dlogD[p])), colours = blue2green2red(100)) +
  scale_y_continuous(expand = c(0,0), 
                     labels = every_nth(y_labels_breaks, 2, inverse = TRUE),
                     breaks = y_labels_breaks) +
  scale_x_datetime(expand = c(0,0), 
                   breaks=date_breaks_2h, 
                   labels=insert_minor(format(date_breaks_1_day, "%d %b"),  
                                       length(date_breaks_1_day))) +
  xlab("Day and time") +
  ylab("Diameter (nm)") +
  mytheme

El gradiente verde-azul es un poco diferente de la trama Igor. ¡Tengo muy poco verde en absoluto! Tal vez una mayor experimentación con rampas de color podría mejorar el juego allí.

Para obtener el eje y en una escala de registro, se requiere un esfuerzo adicional. Tenemos que usar geom_rect y ajustar los tamaños de cada rectángulo para que quepan en la escala de log:

##################  y-axis with log scale ###########################
# get visually diminishing axis ticks
base_breaks <- function(n = 10){
  function(x) {
    axisTicks(log10(range(x, na.rm = TRUE)), log = TRUE, n = n)
  }
}

# Now with log axis, we need to replace the ymin and ymax
distance <- diff((unique(dat_interp_df_long$Diameter)))/2
upper <- (unique(dat_interp_df_long$Diameter)) + c(distance, distance[length(distance)])
lower <- (unique(dat_interp_df_long$Diameter)) - c(distance[1], distance) 

# Create xmin, xmax, ymin, ymax
dat_interp_df_long$xmin <- dat_interp_df_long$Time - 1000 # default of geom_raster is 0.5
dat_interp_df_long$xmax <- dat_interp_df_long$Time + 1000
idx <- rle(dat_interp_df_long$Diameter)$lengths[1]
dat_interp_df_long$ymin <- unlist(lapply(lower, function(i) rep(i, idx)))
dat_interp_df_long$ymax <- unlist(lapply(upper, function(i) rep(i, idx)))


ggplot(dat_interp_df_long, aes(y = Diameter, x = Time, 
                               xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, 
                               fill = dN_dlogDp_log)) +
  geom_rect()  +
  scale_fill_gradientn(name=expression(log(dN/dlogD[p])), colours = blue2green2red(1000)) +
  scale_y_continuous(expand = c(0,0), 
                     trans = log_trans(), breaks = base_breaks()) +
  scale_x_datetime(expand = c(0,0), 
                   breaks=date_breaks_2h, 
                   labels=insert_minor(format(date_breaks_1_day, "%d %b"),  
                                       length(date_breaks_1_day))) +
  xlab("Day and time") +
  ylab("Diameter (nm)") +
  mytheme

ACTUALIZACIÓN Después de algunos experimentos con rampas de color, he encontrado una coincidencia bastante cercana:

# adjust the colour ramp to match the Igor plot (their colour ramp is pretty uneven! lots of red and blue, it seems.)
colfunc <- colorRampPalette(c( rep("red", 3), 
                               rep("yellow", 1), 
                               rep("green", 2), 
                               "cyan", 
                               rep("blue", 3), 
                               "purple"))
y_labels_breaks <- seq(0, max(Diameter), 100)
mytheme <- theme_bw(base_size = 14) +  theme(aspect.ratio = 1/5)
ggplot(dat_interp_df_long, aes(y = Diameter, x = Time,  fill = dN_dlogDp_log)) +
  geom_raster(interpolate = TRUE)  +
  scale_fill_gradientn(name=expression(log(dN/dlogD[p])), colours = rev(colfunc(100))) +
  scale_y_continuous(expand = c(0,0), 
                     labels = every_nth(y_labels_breaks, 2, inverse = TRUE),
                     breaks = y_labels_breaks) +
  scale_x_datetime(expand = c(0,0), 
                   breaks=date_breaks_2h, 
                   labels=insert_minor(format(date_breaks_1_day, "%d %b"),  
                                       length(date_breaks_1_day))) +
  xlab("Day and time") +
  ylab("Diameter (nm)") +
  mytheme

El código de esta publicación también se encuentra en https://gist.github.com/benmarwick/9a54cbd325149a8ff405.

ACTUALIZAR Ahora he hecho un paquete que producirá estas parcelas: https://github.com/benmarwick/smps