Function Documentation

# python
def pythagorize(x, y):
  """
  Calculate and return hypothenuse given x and y
  Parameters
  =========
  x: float
     Distance in x direction
  y: float
     Distance in y direction
  """
  return(x*x + y*y) ** 0.5
print(pythagorize(3, 4))
## 5.0
# R
#' Calculate and return the hypothenuse given x and y
#' @param x distance in x direction
#' @param y distance in y direction
#' @return the diagonal
pythagorize <- function(x, y) {
  return(sqrt(x^2 + y^2))
}
pythagorize(3, 4)
## [1] 5

Import

Slope Analysis Data

Import the slope analysis data.

# R
data_slope_analysis <-
  # file path to slope analysis csv data
  file.path("data", "slope_analysis20180826.csv") %>%
  # read csv data
  read_csv(
    # define column types as safety check
    col_types = cols(
      `Run name` = col_character(),
      `Disturbance rate parameter` = col_double(),
      `Weathering rate parameter` = col_double(),
      `Slope angle` = col_double(),
      `Slope gradient` = col_double()
    )
  ) %>%
  # select and rename relevant columns
  select(
    angle = `Slope angle`,
    wprime = `Weathering rate parameter`,
    dprime = `Disturbance rate parameter`
  ) %>%
  # scale the dimensionless numbers
  mutate(
    wprime = 1e3 * wprime,
    dprime = 1e3 * dprime
  )

# print out data frame
data_slope_analysis

Dissolution Only Data

Import the dissolution only data.

# read dissolution analysis csv data
data_dissolution_only <-
  # file path
  file.path("data", "dissolution_only_results20180717.csv") %>%
  # read csv
  read_csv(
    # skip the first row which has landlab info
    skip = 1,
    # define columns for safety
    col_types = cols(
       `Dissolution rate parameter (1/yr)` = col_double(),
       `Gradient (m/m)` = col_double(),
       `Slope angle (deg)` = col_double()
    )
  ) %>%
  # select and rename relevant columns
  select(
    angle = `Slope angle (deg)`,
    dparam = `Dissolution rate parameter (1/yr)`
  ) %>%
  # scale identically to the model output
  mutate(
    wprime = 1e3 * dparam
  )

# print data
data_dissolution_only

Processing

Analytical solution

Calculate the analytical solution for the dissolution-only case.

analytical_dissolution_only <-
  data_frame(
    wprime = seq(
      from = min(data_dissolution_only$wprime),
      to = max(data_dissolution_only$wprime),
      length.out = 100),
    angle = 60 - 360 * wprime / pi
  )
analytical_dissolution_only

Combine data

Combine all the data into a uniform format for easy calculations and visualization. Introduce new information columns scenario and analytical to make it easy to use grouping and aesthetics later on.

all_data <-
  # combine data frames
  bind_rows(
    # each dprime is one scenario
    data_slope_analysis %>% mutate(scenario = as.character(dprime), analytical = FALSE),
    # dissolution only
    data_dissolution_only %>% mutate(scenario = "inf", analytical = FALSE),
    # analytical solution
    analytical_dissolution_only %>% mutate(scenario = "analytical", analytical = TRUE)
  ) %>%
  # focus on columns of interest
  select(analytical, scenario, angle, wprime) %>%
  # sort
  arrange(analytical, scenario)

# print out
all_data

Information

Spit out summary information about the data.

# summarize how many data points there are scenario and what the min and max angles are
all_data %>%
  group_by(scenario) %>%
  summarize(
    n_data_points = n(),
    angle_min = min(angle),
    angle_max = max(angle)
  )
# summarize at what value of the dimensionless weathering rate the facet angle first falls below 30 degrees (the effective angle of repose)
all_data %>%
  filter(angle < 30) %>%
  group_by(scenario) %>%
  summarize(lowest_w_below_30deg = min(wprime))

Visualization

Plot the facet angle for each value of \(d'\) as a function of weathering parameter (\(w'\)) and add the dissolution only scenario (\(d'\rightarrow \infty\)) and analytical solution for the latter. Also include the model’s effective 30\(^\circ\) angle of repose for reference.

# plot the combined data
all_data %>%
  # focus on data in the relevant range
  filter(between(wprime, 0.1, 100)) %>%
  # start plot object
  ggplot() +
  # define top level aesthetics
  aes(x = wprime, y = angle) +
  # add data points (i.e. everything that is not the analytical solution)
  geom_point(
    data = function(df) filter(df, !analytical), # data filter
    mapping = aes(shape = scenario), # aesthetics
    size = 3 # symbol size
  ) +
  # add analytical solution
  geom_line(
    data = function(df) filter(df, analytical), # data filter
    mapping = aes(color = scenario), # aesthetics
    size = 1 # line thickness
  ) +
  # add horizontal line for 30 degrees
  geom_hline(
    data = data_frame(y = 30, scenario = "eff_angle"),
    mapping = aes(yintercept = y, color = scenario), size = 1, linetype = 2
  ) +
  # switch x-axis to log scale
  scale_x_log10(expand = c(0, 0))

# make the basic plot above interactive
# conversion from ggplot to plotly has some issues but works okay for simple plots
# build plotly (also aesthetics based) from scratch for more complicated plots
ggplotly(p = ggplot2::last_plot(), dynamicTicks = TRUE)

Fancify Plot

# customize the aesthetics of the basic plot above
ggplot2::last_plot() +
  # shape scale & latex legend for symbols
  scale_shape_discrete(
    # replace plain text labels with rendered latex
    labels = function(x) {
      case_when(
        x == "inf" ~ "$d' \\rightarrow \\infty$",
        x != "inf" ~ sprintf("$d' = 10^{%.0f}$", log10(as.numeric(x)))
      ) %>% TeX()
    }
  ) +
  # color scale & latex legend for lines
  scale_color_manual(
    # selecting manual colors
    values = c("red", "dark green"),
    # replace plain text labels with rendered latex
    labels = function(x) {
      case_when(
        x == "analytical" ~ "$\\theta = 60° - 360° \\frac{w'}{\\pi}$",
        x == "eff_angle" ~ "$\\theta = 30°$"
      ) %>% TeX()
    }
  ) +
  # x axis scale
  scale_x_log10(
    labels = function(x) sprintf("$10^{%.0f}$", log10(x)) %>% TeX(),
    expand = c(0, 0)) +
  # y axis scale
  scale_y_continuous(labels = function(x) str_c(x, "°")) +
  # add labels
  labs(
    x = TeX("Dimensionless weathering rate parameter, $w'$"),
    y = "Facet angle",
    shape = "symbols",
    color = "lines"
  ) +
  # switch to better default theme
  theme_bw() +
  # add additional customization
  theme(
    text = element_text(size = 18), # font size
    legend.text.align = 0, # legend text alignment
    panel.grid.minor = element_blank() # remove minor grid lines
  )

Save Plot

ggsave(filename = "facet_angle_vs_wprime.pdf", plot = ggplot2::last_plot(), width = 8, height = 6)

Alternative visualization

This next figure isn’t meant for the paper, but is just here to check out what it looks like if you plot difference between fault dip and facet angle instead of the facet angle itself.

# simply switch out the y aesthetic on the previous plot
ggplot2::last_plot() %+% aes(y = 60 - angle) + labs(y = "60° - facet angle")

Interactive Plot

For more complex interactive plots, it’s useful to call plot_ly directly.

# start plotly object
plot_ly(
  x = ~ wprime,
  y = ~ angle
) %>%
  # add points
  add_markers(
    data = filter(all_data, between(wprime, 0.1, 100), !analytical),
    color = ~ str_c("d' = ", scenario),
    marker = list(size = 10)
  ) %>%
  # add lines
  add_lines(
    data = filter(all_data, between(wprime, 0.1, 100), analytical),
    color = ~ scenario,
    line = list(color = "black")
  ) %>%
  # add layout specifications
  layout(
    xaxis = list(title = "Dimensionless weathering rate parameter", type = "log"),
    yaxis = list(title = "facet angle")
  )

Animated Plot

# animation frames
frames <-
  data_frame(
    wprime_cutoff = 10^seq(from = -1, to = 2, by = 0.05),
    `w'` = round(wprime_cutoff, 3)
  )

# animated data
animated_data <- all_data %>%
  # focus on data between 0.1 and 100 wprime
  filter(between(wprime, 0.08, 100)) %>%
  # cross with the frames data frame
  crossing(frames) %>%
  # ensure each frame only goes to the cutoff
  filter(wprime <= wprime_cutoff)
# generate animated plot
# Note: also works by defining the frame aesthetic in ggplot and using ggplotly
  plot_ly(
    x = ~ wprime,
    y = ~ angle,
    frame = ~ `w'`
  )  %>%
  # add points
  add_markers(
    data = filter(animated_data, !analytical),
    color = ~ str_c("d' = ", scenario),
    marker = list(size = 10)
  ) %>%
  # add lines
  add_lines(
    data = filter(animated_data, analytical),
    color = ~ scenario,
    line = list(color = "black")
  ) %>%
  # add layout specifications
  layout(
    title = "<b>Wonderfuly Unnecessarily Animated Plot</b>",
    xaxis = list(title = "Dimensionless weathering rate parameter", type = "log"),
    yaxis = list(title = "facet angle")
  ) %>%
  # animation options
  animation_opts(
    frame = 100,
    transition = 0,
    redraw = FALSE
  )