Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add visualization #8

Open
cells2numbers opened this issue Sep 8, 2017 · 1 comment
Open

add visualization #8

cells2numbers opened this issue Sep 8, 2017 · 1 comment
Assignees

Comments

@cells2numbers
Copy link
Owner

Add visualization including

  • roseplot / windrose plot
  • trajectory plot 2d/3d
@cells2numbers
Copy link
Owner Author

Example code to load cell profiler output, track cells and plot windrose plot to visualize migration pattern.

# define column name for track label
strata <- 'TrackObjects_Label'

# read cellprofiler output
profiles <- read_csv("../profiles/201903_testsequence/IdentifyPrimaryObjects.csv") %>% 
  mutate(Metadata_timepoint = as.numeric(Metadata_timepoint)) %>%
  group_by(TrackObjects_Label) %>% 

# summarize tracks 
tracks <- track(profiles, strata = strata,  t_var = "Metadata_timepoint"  ) 

# windrose plot 
plot_windrose(tracks)
plot_windrose <- function(population, title_name = "windrose plot", legend = "speed in pixel/frame", dir_res = 30, spd_res = 0.5){
  # rotate track angles 
  population <- tracks %>%
    filter(Track_Length > 19) %>%
    mutate(Track_Angle = Track_Angle + pi/2) %>% # rotate all angles by 90 degree or pi/2
    mutate(Track_Angle = ifelse(Track_Angle > pi, Track_Angle - 2*pi, Track_Angle) ) %>%
    mutate(Track_Angle = ifelse(Track_Angle < 0, Track_Angle + 2*pi, Track_Angle) ) 

  
  h1 <- plot.windrose(
    spd = population %>% extract2("Track_Speed"), 
    dir = (180 * (population$Track_Angle) / pi), 
    spdmin = 0, 
    spdmax = max(population$Track_Speed), 
    spdres = spd_res, 
    dirres = dir_res, 
    title_name =  title_name,
    scale_name = legend
  )
} 

plot_windrose(tracks)

image

Code for the complete windrose plot comes from following https://stackoverflow.com/questions/17266780/wind-rose-with-ggplot-r/17266781#17266781

# WindRose.R
#
# following https://stackoverflow.com/questions/17266780/wind-rose-with-ggplot-r/17266781#17266781
require(RColorBrewer)

plot.windrose <- function(data,
  spd,
  dir,
  spdres = 2,
  dirres = 30,
  spdmin = 2,
  spdmax = 20,
  spdseq = NULL,
  palette = "YlGnBu",
  countmax = NA,
  debug = 0,
  scale_name = "speed in µm/s",
  title_name = "Distribution of angle and speed"){
  
  
  # Look to see what data was passed in to the function
  if (is.numeric(spd) & is.numeric(dir)){
    # assume that we've been given vectors of the speed and direction vectors
    data <- data.frame(spd = spd,
      dir = dir)
    spd = "spd"
    dir = "dir"
  } else if (exists("data")){
    # Assume that we've been given a data frame, and the name of the speed 
    # and direction columns. This is the format we want for later use.    
  }  
  
  # Tidy up input data ----
  n.in <- NROW(data)
  dnu <- (is.na(data[[spd]]) | is.na(data[[dir]]))
  data[[spd]][dnu] <- NA
  data[[dir]][dnu] <- NA
  
  # figure out the wind speed bins ----
  if (missing(spdseq)){
    spdseq <- seq(spdmin,spdmax,spdres)
  } else {
    if (debug >0){
      cat("Using custom speed bins \n")
    }
  }
  # get some information about the number of bins, etc.
  n.spd.seq <- length(spdseq)
  n.colors.in.range <- n.spd.seq - 1
  
  # create the color map
  spd.colors <- colorRampPalette(brewer.pal(min(max(3,
    n.colors.in.range),
    min(9,
      n.colors.in.range)),                                               
    palette))(n.colors.in.range)
  
  if (max(data[[spd]],na.rm = TRUE) > spdmax){    
    spd.breaks <- c(spdseq,
      max(data[[spd]],na.rm = TRUE))
    spd.labels <- c(paste(c(spdseq[1:n.spd.seq-1]),
      '-',
      c(spdseq[2:n.spd.seq])),
      paste(spdmax,
        "-",
        max(data[[spd]],na.rm = TRUE)))
    spd.colors <- c(spd.colors, "grey50")
  } else{
    spd.breaks <- spdseq
    spd.labels <- paste(c(spdseq[1:n.spd.seq-1]),
      '-',
      c(spdseq[2:n.spd.seq]))    
  }
  data$spd.binned <- cut(x = data[[spd]],
    breaks = spd.breaks,
    labels = spd.labels,
    ordered_result = TRUE)
  # clean up the data
  data. <- na.omit(data)
  
  # figure out the wind direction bins
  dir.breaks <- c(-dirres/2,
    seq(dirres/2, 360-dirres/2, by = dirres),
    360+dirres/2)  
  # 
  dir.labels <- c(paste(360-dirres/2,"-",dirres/2),
    paste(seq(dirres/2, 360-3*dirres/2, by = dirres),
      "-",
      seq(3*dirres/2, 360-dirres/2, by = dirres)),
    paste(360-dirres/2,"-",dirres/2))
  
  #dir.labels <- c(paste(0),
  #  paste(seq(dirres/2, 360-3*dirres/2, by = dirres)),
  #"-",
  #seq(3*dirres/2, 360-dirres/2, by = dirres)),
  #  paste(360-dirres/2,"-",dirres/2))
  
  # assign each wind direction to a bin
  dir.binned <- cut(data[[dir]],
    breaks = dir.breaks,
    ordered_result = TRUE)
  levels(dir.binned) <- dir.labels
  data$dir.binned <- dir.binned
  
  # Run debug if required ----
  if (debug>0){    
    cat(dir.breaks,"\n")
    cat(dir.labels,"\n")
    cat(levels(dir.binned),"\n")       
  }  
  
  # deal with change in ordering introduced somewhere around version 2.2
  if(packageVersion("ggplot2") > "2.2"){    
    #cat("Hadley broke my code\n")
    data$spd.binned = with(data, factor(spd.binned, levels = rev(levels(spd.binned))))
    spd.colors = rev(spd.colors)
  }
  
  # create the plot ----
  p.windrose <- ggplot(data = data,
    aes(x = dir.binned,
      fill = spd.binned)) +
    geom_bar() + 
    scale_x_discrete(drop = FALSE,
      labels = waiver()) +
    coord_polar(start = -((dirres/2)/360) * 2*pi) +
    scale_fill_manual(name = scale_name, 
      values = spd.colors,
      drop = FALSE) +
    theme(axis.title.x = element_blank()) +
    labs(title = title_name)
  
  # adjust axes if required
  if (!is.na(countmax)){
    p.windrose <- p.windrose +
      ylim(c(0,countmax))
  }
  
  # print the plot
  print(p.windrose)  
  
  # return the handle to the wind rose
  return(p.windrose)
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

1 participant