R/plot_c4_aci_fit.R

In PhotoGEA: Photosynthetic Gas Exchange Analysis

plot_c4_aci_fit <- function(
    fit_results,
    identifier_column_name,
    x_name,
    plot_operating_point = TRUE,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    pcm_column_name = 'PCm',
    ...
)
{
    if (length(fit_results) != 3 || !is.exdf(fit_results$fits) ||
        !is.exdf(fit_results$parameters) || !(is.exdf(fit_results$fits_interpolated))) {
        stop('fit_results must be the output from fit_c4_aci')
    }

    # Make sure the x variable is acceptable
    if (!x_name %in% c(ci_column_name, pcm_column_name)) {
        stop('The x_name must be the Ci or PCm column name')
    }

    # Get the appropriate operating point name
    operating_x <- if (x_name == ci_column_name) {
        'operating_Ci'
    } else {
        'operating_PCm'
    }

    # Make sure the required variables are defined and have the correct units
    required_variables <- list()
    required_variables[[ci_column_name]]         <- 'micromol mol^(-1)'
    required_variables[[identifier_column_name]] <- NA
    required_variables[[pcm_column_name]]        <- 'microbar'

    # Don't throw an error if some columns are all NA
    check_required_variables(fit_results$fits_interpolated, required_variables, check_NA = FALSE)

    required_variables[[a_column_name]] <- unit_dictionary('A')

    # Don't throw an error if some columns are all NA
    check_required_variables(fit_results$fits, required_variables, check_NA = FALSE)

    required_variables <- list()
    required_variables[['operating_Ci']]       <- 'micromol mol^(-1)'
    required_variables[['operating_PCm']]      <- 'microbar'
    required_variables[['operating_An_model']] <- unit_dictionary('A')

    # Don't throw an error if some columns are all NA
    check_required_variables(fit_results$parameters, required_variables, check_NA = FALSE)

    # Choose line settings
    assim_cols <- multi_curve_colors()[1:5]
    assim_cols[1] <- '#676767'

    line_settings <- list(
        col = assim_cols,
        lwd = c(4, 2, 2, 2, 2),
        lty = c(1, 5, 5, 5, 5)
    )

    # Plot the fits, operating point, and raw data
    lattice::xyplot(
        An + Apc + Apr + Ar + Aj ~ fit_results$fits_interpolated[, x_name] | fit_results$fits_interpolated[, identifier_column_name],
        data = fit_results$fits_interpolated$main_data,
        type = 'l',
        par.settings = list(superpose.line = line_settings),
        auto.key = list(space = 'right', lines = TRUE, points = FALSE),
        xlab = paste(x_name, '[', fit_results$fits_interpolated$units[[x_name]], ']'),
        ylab = paste(
            'Net CO2 assimilation rate [', fit_results$fits_interpolated$units[['An']],
            ']\n(filled black circles: measured data used for fits',
            '\nopen black circles: measured data excluded from fits',
            '\nopen red circle: estimated operating point)'
        ),
        curve_ids = fit_results$fits_interpolated[, identifier_column_name],
        panel = function(...) {
            # Get info about this curve
            args <- list(...)
            curve_id <- args$curve_ids[args$subscripts][1]

            curve_parameters <-
                fit_results$parameters[fit_results$parameters[, identifier_column_name] == curve_id, ]

            curve_data <-
                fit_results$fits[fit_results$fits[, identifier_column_name] == curve_id, ]

            used_for_fit <- points_for_fitting(curve_data)

            # Plot the fit lines
            lattice::panel.xyplot(...)

            # Plot the operating point, if desired
            if (plot_operating_point) {
                lattice::panel.points(
                    curve_parameters[1, 'operating_An_model'] ~ curve_parameters[1, operating_x],
                    col = 'red',
                    pch = 1
                )
            }

            # Plot the measured data points
            lattice::panel.points(
                curve_data[used_for_fit, a_column_name] ~ curve_data[used_for_fit, x_name],
                col = 'black',
                pch = 16
            )

            lattice::panel.points(
                curve_data[!used_for_fit, a_column_name] ~ curve_data[!used_for_fit, x_name],
                col = 'black',
                pch = 1
            )
        },
        ...
    )
}