用R重建Igor的等高线图

Question

这个等高线图是用Igor程序制作的，在大气化学和污染研究中很流行：

​

​

我试着用R为一个想停止使用Igor的朋友重新创建它，但是我们不能完全得到它。下面是数据集 (使用Igor制作这个图所用的相同数据)，下面是我到目前为止用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

​

​

我的情节可以用更巧妙的标签和勾标等来完成。然而，我的主要问题是为什么我的轮廓填充看起来与Igor图如此不同。比例似乎颠倒了，插值看起来很不一样。

我怎样才能让我的情节看起来更像伊戈尔的情节？

请注意，我的这些其他问题与重新创建这个情节的任务密切相关：

• 具有对数尺度的栅格插值
• 类别的2d密度图

在我提出这个问题之后，我一直在更新R代码的要点，它结合了对这些问题的回答的细节，并成功地复制了这些图( gist中包含的示例输出)。这个要点就在这里：https://gist.github.com/benmarwick/9a54cbd325149a8ff405。

更新--我现在做了一个包，它将生成这些情节：https://github.com/benmarwick/smps

Answer 1

我可以使用akima::interp而不是loess来进行插值，从而更接近Igor图：

# 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"))

​

​

但是在颜色映射上仍然有很大的差别，Igor图有宽的带和尖锐的边界，而我的图有较少的颜色带和模糊的边界。所以我想我还没有伊戈尔图所用的插值方法。

更新了，在试验了一系列彩色坡道之后，我在colorRamps::blue2green2red中找到了一个很好的匹配。我在这里也花了一点精力在花哨的刻痕上：

# 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

​

​

绿色-蓝色梯度与Igor图仍有一点不同。我一点绿色都没有！也许对彩色坡道的进一步试验可能会改善那里的匹配。

要使y轴在一个日志尺度上，需要额外的努力。我们必须使用geom_rect并调整每个矩形的大小以适应日志比例尺：

##################  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

​

​

在对彩色坡道进行了一些实验之后，更新了，我发现了一个非常接近的匹配：

# 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

​

​

这篇文章中的代码也在https://gist.github.com/benmarwick/9a54cbd325149a8ff405

更新--我现在做了一个包，它将生成这些情节：https://github.com/benmarwick/smps