Package 'EnvRtype'

Description

Returns the URL template get_AEZ uses when no local file is given and the bundled raster is unavailable.

Usage

aez_default_url(layer = "aez57")

Arguments

Details

This URL is unverified. It was constructed from the documented GAEZ v4 portal structure but could not be tested against a live FAO server, and FAO has reorganised these paths between releases. It is exposed as a function, rather than hidden inside get_AEZ(), so it can be inspected and overridden without editing the package:

  get_AEZ(env, lat, lon, url = "https://.../correct/path.tif")

For reproducible analyses, rely on the bundled raster or download the layer once and pass file=.

Value

A single character string.

See Also

get_AEZ

Examples

aez_default_url()

Description

Returns the class table used to decode GAEZ agro-ecological zone raster values into human-readable zone names.

Usage

aez_legend(classes = 57)

Arguments

Details

Verify before relying on these labels. FAO has revised AEZ class definitions between GAEZ releases. The table here follows the GAEZ v4 documentation, but if your raster came from a different release the codes may not line up. get_AEZ will warn when a raster contains codes that are absent from the legend – that warning is the signal to check.

Value

A data.frame with code, label, thermal, moisture and arable (logical; FALSE for water, urban, barren and the constraint classes).

See Also

get_AEZ

Examples

leg <- aez_legend()
head(leg)
table(leg$thermal)

Description

Returns the long phenotype table, so a sim_met object can be handed straight to modelling functions or inspected as a plain frame.

Usage

## S3 method for class 'sim_met'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A data.frame with one row per plot and columns env, gid, rep and value.

See Also

sim_met

Description

Coerce a Simulated Envirome Diversity Matrix to a Plain Matrix

Usage

## S3 method for class 'sim_envirome'
as.matrix(x, ...)

Arguments

Value

A numeric matrix with no sim_envirome class or attributes.

See Also

sim_envirome

Description

Strips the sim_W class and its diagnostic attributes, returning the bare $q \times n_{var}$ covariable matrix for use with env_kernel, W_matrix consumers, or base matrix ops.

Usage

## S3 method for class 'sim_W'
as.matrix(x, ...)

Arguments

Value

A numeric matrix (environments in rows, covariables in columns) with no sim_W class or attributes.

See Also

sim_W

Description

Best planting date per genotype

Usage

best_planting_date(
  scan,
  dates = NULL,
  date_col = "date",
  genotypes = NULL,
  group = NULL,
  exclude_flagged = TRUE,
  envelope_level = c("gross", "mild")
)

Arguments

Value

data.frame: 'gid', 'best_date', 'yHat', 'lower', 'upper', 'position', 'n_dates_considered'.

Examples

## A scan is normally produced by scan_untested_envs(); here we build a
## minimal one by hand -- 3 genotypes across 6 candidate sowing windows --
## so the example is self-contained.
gids <- paste0("G", 1:3)
envs <- paste0("sow", 1:6)
mu <- matrix(c(6, 5.5, 5), 3, 6) +
      matrix(rep(c(0, .4, .7, .6, .2, -.3), each = 3), 3, 6)  # peaks mid-window
dimnames(mu) <- list(gids, envs)
scan <- list(
  yHat_matrix = mu, yHat_lower = mu - 0.5, yHat_upper = mu + 0.5,
  interpolation = data.frame(
    env               = envs,
    position          = c("interpolation", "interpolation", "edge",
                          "interpolation", "mild_extrapolation",
                          "gross_extrapolation"),
    mahalanobis_ratio = c(0.6, 0.8, 1.4, 0.9, 2.1, 3.6),
    prop_outside      = c(0, 0, 0.1, 0, 0.3, 0.7),
    weight_negativity = 0, stringsAsFactors = FALSE))
class(scan) <- "scan_untested_envs"

## Map each scanned environment to a calendar sowing date
dates <- data.frame(env = envs,
                    date = as.Date("2024-04-01") + c(0, 14, 28, 42, 56, 70))

## Best *supported* sowing date per genotype (gross extrapolations dropped)
best_planting_date(scan, dates)

## Rank every date, including the out-of-envelope ones
best_planting_date(scan, dates, exclude_flagged = FALSE)

## Stricter: also drop mild extrapolations before ranking
best_planting_date(scan, dates, envelope_level = "mild")

Description

Soil profile for clay content (g/kg) around Nairobi, downloaded from SoilGrids (https://soilgrids.org/). Shipped as a GeoTIFF under inst/extdata and read with terra.

Format

A single-layer GeoTIFF raster read with terra::rast().

Examples

f <- system.file("extdata", "clay_5_15.tif", package = "EnvRtype")
clay <- terra::rast(f)
terra::plot(clay)

Description

Delineates mega-environments by clustering the genetic correlation among environments, i.e. how similarly genotypes rank from one environment to another. Environments in which genotypes respond alike are placed in the same cluster, which is the operational definition of a mega-environment in a target population of environments (TPE).

The input is either a fitted kernel_model (carrying the environment genetic-correlation matrix in $genetic_cor_env, produced with keep_effects = TRUE) or any square environment correlation matrix you already have (for example from env_cor or an external GxE analysis).

Usage

cluster_environments(
  object,
  k = NULL,
  method = c("hclust", "pam", "kmeans"),
  distance = c("one_minus", "sqrt"),
  hclust_linkage = "ward.D2",
  k_range = 2:6,
  seed = NULL
)

Arguments

Details

The genetic correlation $r_g$ is turned into a distance and the chosen algorithm is run on that distance:

• distance = "one_minus" uses $d = 1 - r_g$ (default).

• distance = "sqrt" uses $d = \sqrt{1 - r_g}$, which spreads out highly correlated environments.

Correlations are clamped to $[-1, 1]$ and the distance matrix is symmetrised before clustering. Three algorithms are available through method: agglomerative hierarchical clustering ("hclust", controlled by hclust_linkage), partitioning around medoids ("pam"), and "kmeans" run on a classical MDS embedding of the distance matrix.

When k is NULL the number of clusters is selected automatically by maximising the average silhouette width over the candidate values in k_range; supplying k skips this search. At least three environments are required.

Value

A list with:

clusters

named integer vector giving the cluster of each environment.

k

the number of clusters used (selected or supplied).

method

the clustering algorithm used.

distance

the symmetric distance matrix that was clustered.

silhouette

per-environment silhouette widths plus their average (higher is a tighter, better-separated solution).

medoids

representative environment of each cluster ("pam"/"kmeans"; NULL for "hclust").

hclust

the hclust tree for method = "hclust", otherwise NULL.

table

a data.frame (env, cluster, sil_width) ordered by cluster and decreasing silhouette width.

See Also

kernel_model (fit with keep_effects = TRUE to obtain $genetic_cor_env), kernel_model_clustered to fit a model informed by the clusters, env_cluster for clustering from an environmental covariable matrix, and env_cor.

Examples

## ------------------------------------------------------------------
## 1. Cluster a genetic-correlation matrix directly (no model needed)
## ------------------------------------------------------------------
## Six environments forming two obvious groups
R <- matrix(0.2, 6, 6)
R[1:3, 1:3] <- 0.9
R[4:6, 4:6] <- 0.85
diag(R) <- 1
dimnames(R) <- list(paste0("E", 1:6), paste0("E", 1:6))

## Automatic number of clusters via silhouette width
cl <- cluster_environments(R)
cl$k
cl$table
cl$silhouette["average"]

## ------------------------------------------------------------------
## 2. Force a specific number of clusters
## ------------------------------------------------------------------
cluster_environments(R, k = 2)$clusters

## ------------------------------------------------------------------
## 3. Try different algorithms and distances
## ------------------------------------------------------------------
cluster_environments(R, k = 2, method = "pam", seed = 1)$medoids
cluster_environments(R, method = "kmeans", k_range = 2:4, seed = 1)$clusters
cluster_environments(R, distance = "sqrt",
                     hclust_linkage = "complete")$clusters


## ------------------------------------------------------------------
## 4. Full workflow from phenotypes, kernels and a fitted model
## ------------------------------------------------------------------
data("maizeYield"); data("maizeG"); data("maizeWTH")

ECs <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
                var.id = c("FRUE", "PETP", "SRAD", "T2M_MAX"))
K <- get_kernel(K_G = list(G = maizeG),
                K_E = list(W = env_kernel(env.data = ECs)[[2]]),
                data = maizeYield, model = "RNMM")

## keep_effects = TRUE is required for $genetic_cor_env
fit <- kernel_model(y = "value", env = "env", gid = "gid",
                    data = maizeYield, random = K,
                    keep_effects = TRUE)

megaenv <- cluster_environments(fit)
megaenv$table

## 5. Feed the clusters into a cluster-aware model
fit_cl <- kernel_model_clustered(y = "value", data = maizeYield,
                                 random = K, env = "env", gid = "gid",
                                 clusters = megaenv, use = "block")

Description

Tabulates how much of a scanned area falls in each envelope class, and which covariates most often drive sites outside. The second table answers "what should I go and measure".

Usage

coverage_summary(scan)

Arguments

Value

list of 'by_position' and 'by_covariate' data.frames.

Description

Turns each dense kernel (list(Kernel = <n x n>, Type = "D"|"BD")) into a cached spectral decomposition stored on $decomp. Call once and reuse across any number of kernel_model fits: kernels depend only on genotypes and environments – never on the phenotype vector – so masking y for a cross-validation fold does NOT invalidate the decomposition.

Usage

decompose_kernels(
  K,
  ne = NULL,
  tol = 1e-10,
  scale_kernels = TRUE,
  keep_var = 1,
  verbose = TRUE
)

Arguments

Details

Factored kernels from get_kernel arrive pre-decomposed and do not need this function.

Value

K with $decomp on each element and a "kd_meta" attribute.

See Also

get_kernel, kernel_model, undecompose_kernels

Description

Clusters environments into homogeneous groups ("mega-environments") from an environmental covariable matrix (W) using K-means. The number of groups can be set by the user or selected automatically by maximising the average silhouette width. An environmental relationship (kinship) matrix is also returned for convenience.

Usage

env_cluster(
  W,
  k = NULL,
  k.max = 10,
  scale = TRUE,
  nstart = 25,
  seed = NULL,
  verbose = TRUE
)

Arguments

Details

Mega-environments are repeatable groups of locations that rank genotypes similarly. Here they are approximated by grouping environments that share a similar enviromic profile, which is useful to stratify a Target Population of Environments (TPE) before multi-environment prediction. The silhouette width $s(i)$ contrasts within- and between-cluster distances; the k maximising its average is retained.

Value

A list with clusters (named integer vector), k, silhouette (data.frame of average silhouette per tested k, or NULL when k was supplied), kmeans (the fitted object) and env.kinship (an environment x environment correlation matrix).

Author(s)

Germano Costa Neto

See Also

W_matrix, env_kernel, env_target_importance, T_matrix

Examples

data("maizeWTH")

W <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
              var.id = c("T2M", "T2M_MAX", "PRECTOT"), statistic = "mean")

## Automatic number of mega-environments (chosen by silhouette)
mega <- env_cluster(W, seed = 1)
mega$k
mega$clusters
mega$silhouette

## Force a fixed number of groups
mega2 <- env_cluster(W, k = 2, seed = 1)
mega2$clusters

## Environmental relationship (kinship) matrix
round(mega$env.kinship, 2)

## Use the grouping downstream (e.g. as a fixed effect or to stratify CV)
table(mega$clusters)

Description

Performs a copula analysis of an environment x covariable matrix (W) and returns a matrix of copula indices that can replace the raw covariables. Sklar's theorem lets the joint distribution of the environmental variables be split into (i) the marginal behaviour of each variable and (ii) the dependence structure linking them (the copula). By moving to the copula scale, each environment is described by where it sits in the multivariate environmental distribution rather than by raw physical units.

Usage

env_copula(
  W,
  index = c("pobs", "joint", "survival", "kendall", "all"),
  groups = NULL,
  ties.method = "average",
  normalize = TRUE,
  verbose = TRUE
)

Arguments

Details

Why rank-based. All indices are computed from ranks, so they are invariant to any monotone transformation of the marginals: storing precipitation in mm or log(mm) gives the same result, and no distributional assumption is imposed. This is the main practical gain over a z-scored W, where one extreme site can dominate the Euclidean geometry.

Pseudo-observations. For $q$ environments and variable $j$,

$u_{ij} = r_{ij}/(q+1)$

with $r_{ij}$ the rank of environment $i$. Dividing by $q+1$ keeps values strictly inside $(0,1)$, avoiding boundary problems in downstream density evaluation.

Empirical copula. The joint index is the empirical copula evaluated at each environment's own pseudo-observation, a nonparametric estimator of $C$ that requires no copula family to be chosen.

Kendall function. $K_C$ is the CDF of $Z = C(U)$, estimated by the Genest-Rivest empirical construction. It collapses a $d$-dimensional dependence structure into one interpretable number per environment and underpins multivariate return periods in hydrology.

Sample-size caveat. These are rank statistics across environments, so their resolution is bounded by $q$. With five environments each pseudo-observation can only take values in 1/6,...,5/6 and the joint copula is very coarse; the indices become informative with dozens of environments. A warning is issued when $q < 10$.

Scope. The dependence modelled is that among environmental variables across the environment panel. Geographic coordinates are not used and no spatial autocorrelation (variogram / distance decay) is fitted. Use groups to control which variables share a copula.

Value

A numeric matrix with $q$ environments in rows (rownames preserved from W) and the requested copula indices in columns. Attributes "copula.index" and "copula.groups" record the index produced and the block structure used.

Author(s)

Germano Costa Neto

References

Sklar, A. (1959). Fonctions de repartition a n dimensions et leurs marges. Publ. Inst. Statist. Univ. Paris 8, 229-231.

Genest, C., Rivest, L.-P. (1993). Statistical inference procedures for bivariate Archimedean copulas. JASA 88, 1034-1043.

Salvadori, G., De Michele, C., Kottegoda, N.T., Rosso, R. (2007). Extremes in Nature: An Approach Using Copulas. Springer.

See Also

W_matrix, env_kernel, env_cluster

Examples

data("maizeWTH")
env.data <- maizeWTH[maizeWTH$daysFromStart < 100, ]

W <- W_matrix(env.data = env.data, var.id = c("T2M", "T2M_MAX", "PRECTOT"),
              statistic = "mean")

## 1. Marginal copula scale (pseudo-observations)
round(env_copula(W, index = "pobs"), 3)

## 2. Joint non-exceedance probability per environment
env_copula(W, index = "joint")

## 3. Compound-stress (joint exceedance) index
env_copula(W, index = "survival")

## 4. Kendall multivariate return level
env_copula(W, index = "kendall")

## 5. Blocks: heat and water handled as separate copulas
env_copula(W, index = "all",
           groups = list(heat  = c("T2M_mean", "T2M_MAX_mean"),
                         water = "PRECTOT_mean"))

## 6. Straight from W_matrix() via the 'copula' argument
Wc <- W_matrix(env.data = env.data, var.id = c("T2M", "T2M_MAX", "PRECTOT"),
               statistic = "mean", copula = "kendall",
               center = FALSE, scale = FALSE)

## 7. Copula-scale kernel for reaction-norm models
Kc <- env_kernel(env.data = env_copula(W, index = "pobs"), gaussian = TRUE)
round(Kc$envCov, 2)

Description

Accessory generic that pulls the environment correlation out of a sim_met object or a bare matrix, so a hand-off to sim_W never depends on internal structure.

Usage

env_cor(x, type = c("target", "realised"), ...)

## S3 method for class 'sim_met'
env_cor(x, type = c("target", "realised"), ...)

## S3 method for class 'matrix'
env_cor(x, type = c("target", "realised"), ...)

## S3 method for class 'sim_W'
env_cor(x, type = c("target", "realised"), ...)

## Default S3 method:
env_cor(x, type = c("target", "realised"), ...)

## S3 method for class 'sim_envirome'
env_cor(x, type = c("target", "realised"), ...)

Arguments

Value

A $q \times q$ correlation matrix.

See Also

sim_met, sim_W

Examples

met <- sim_met(maizeG, n_env = 6, seed = 1)
env_cor(met, "target")
env_cor(met, "realised")

Description

Pearson correlation among all environmental indices, as in Fig. 1b of Della Coletta et al. (2023), which motivated the dimensionality reduction: the 867 indices were highly correlated.

Usage

env_cor_heatmap(
  W,
  order = c("factor", "interval", "none"),
  max.index = 400L,
  title = "Correlation among environmental indices",
  verbose = TRUE
)

Arguments

Details

The mean absolute off-diagonal correlation is printed as a one-number summary of how redundant the index set is – the quantity that justifies running a PCA at all.

Value

Invisibly, the correlation matrix that was plotted.

See Also

env_indices, env_pca

Examples

W <- env_indices(maizeWTH, end.day = 130)
R <- env_cor_heatmap(W)

Description

Expands an environment-level W matrix to the $n$ observations of a phenotypic data.frame, so that each row of the output corresponds to one genotype-in-environment record.

Usage

env_expand(env.data, df.pheno, env.id = "env", skip = 3, verbose = TRUE)

Arguments

Value

A numeric matrix of $n$ rows (observations) by $k$ covariables.

See Also

W_matrix, env_kernel

Examples

data("maizeYield"); data("maizeWTH")
W <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
              var.id = c("T2M", "PRECTOT"), statistic = "mean")
Wn <- env_expand(env.data = W, df.pheno = maizeYield, env.id = "env")
dim(Wn)

Description

Converts daily environmental data into the environment x index matrix used by Della Coletta et al. (2023): each weather factor is summarised within fixed-width time intervals from planting, producing one column per factor x interval combination.

Usage

env_indices(
  env.data,
  env.id = "env",
  day.id = "daysFromStart",
  var.id = NULL,
  interval = 3L,
  end.day = 151L,
  statistic = c("mean", "sum", "min", "max", "median"),
  statistic.by = NULL,
  id.names = NULL,
  sep = "_",
  require.complete = TRUE,
  verbose = TRUE
)

Arguments

Details

Why a fixed end day. Environments differ in season length. Without a common cut-off, environments contribute different numbers of intervals and the matrix is ragged. The paper fixed the last interval at 151 days after planting; the same device is used here via end.day.

Interval labelling. With interval = 3 and end.day = 151, days 1-3 form interval 1, days 4-6 interval 2, and so on; the final partial window (day 151) forms interval 51. This reproduces the 51 intervals of the paper, and with 17 factors gives the 867 indices reported.

Mean vs sum. The paper does not state the within-interval statistic. A mean is used by default because it is defined for every factor, but for accumulating quantities such as precipitation a sum is more interpretable; statistic.by allows a per-factor choice.

Value

A numeric matrix with environments in rows and $n_{factor} \times n_{interval}$ indices in columns. Attributes:

"factors"

the factor names used

"intervals"

the interval labels

"interval.mid"

named numeric vector of interval mid-point days

"interval.range"

data.frame of interval start/end days

"n.missing"

number of missing cells

Author(s)

Implementation following Della Coletta et al. (2023)

References

Della Coletta R., Liese S.E., Fernandes S.B., Mikel M.A., Bohn M.O., Lipka A.E., Hirsch C.N. (2023). Linking genetic and environmental factors through marker effect networks to understand trait plasticity. GENETICS 224(4), iyad103.

See Also

env_pca, env_cor_heatmap

Examples

data("maizeWTH")

## 3-day intervals across the shared season. The bundled sample does not
## reach the paper's 151 days in every environment, so end.day is trimmed to
## a fully observed window (the default 151 would leave gaps and error out).
W <- env_indices(maizeWTH, env.id = "env", day.id = "daysFromStart",
                 end.day = 130)
dim(W)
attr(W, "intervals")[1:5]

## Sum precipitation within each window, average everything else
W2 <- env_indices(maizeWTH, statistic = "mean",
                  statistic.by = c(PRECTOT = "sum"), end.day = 130)

## Weekly windows instead
W3 <- env_indices(maizeWTH, interval = 7, end.day = 140)

Description

Returns environmental kinships for reaction-norm models. The output is a list containing varCov, the relatedness among environmental covariables, and envCov, the relatedness among environments. Linear (GBLUP-like), nonlinear Gaussian and arc-cosine deep kernels are supported.

Usage

env_kernel(
  env.data,
  Y = NULL,
  is.scaled = TRUE,
  sd.tol = 10,
  digits = 5,
  tol = 0.001,
  merge = FALSE,
  Z_E = NULL,
  stages = NULL,
  env.id = "env",
  gaussian = FALSE,
  h.gaussian = NULL,
  deep.kernel = FALSE,
  deep.layers = 1,
  QC = FALSE,
  verbose = TRUE
)

Arguments

Details

Three kernel methods are available.

Linear (GB). $K = WW' / \mathrm{trace}(WW')/n$, the environmental analogue of the GBLUP kernel of VanRaden (2008).

Gaussian (GK). $K = \exp(-h\,d/q)$, with $d$ the squared Euclidean distance and $q$ its median, so the bandwidth is scale-free.

Arc-cosine deep (DK). Emulates a deep neural network with deep.layers hidden layers (Cuevas et al. 2019); no bandwidth parameter is required.

Note that quality control is now controlled by the explicit QC argument. In earlier versions it ran only when is.scaled = FALSE, so with the default is.scaled = TRUE no filtering happened despite sd.tol being documented as active.

Value

A list with two elements: varCov, the covariable x covariable relatedness, and envCov, the environment x environment relatedness. When stages is supplied, each element is itself a named list of kernels, one per stage.

Author(s)

Germano Costa Neto

References

Cuevas J. et al. (2019). Deep kernel for genomic and near infrared predictions in multi-environment breeding trials. G3 9(9), 2913-2924.

Costa-Neto G., Fritsche-Neto R., Crossa J. (2021). Nonlinear kernels, dominance, and envirotyping data increase the accuracy of genome-based prediction in multi-environment trials. Heredity 126, 92-106.

See Also

W_matrix, get_kernel, env_cluster, T_matrix

Examples

data('maizeYield'); data("maizeWTH")

## Environmental covariable matrix (environments x covariables)
W.cov <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
                  var.id = c("T2M", "T2M_MAX", "PRECTOT"), statistic = "mean")

## 1. Linear (co)variance kernel
K.lin <- env_kernel(env.data = W.cov, gaussian = FALSE)
round(K.lin$envCov, 2)

## 2. Nonlinear Gaussian kernel
K.gau <- env_kernel(env.data = W.cov, gaussian = TRUE)
round(K.gau$envCov, 2)

## 3. Arc-cosine deep kernel with 2 hidden layers
K.deep <- env_kernel(env.data = W.cov, deep.kernel = TRUE, deep.layers = 2)
round(K.deep$envCov, 2)

## 4. One kernel per development stage
stages   <- c('VE', 'V1_V6', 'V6_VT', 'VT_R1')
interval <- c(0, 7, 30, 65, 90)
W.stage  <- W_matrix(env.data = maizeWTH, var.id = c('FRUE', 'PETP'),
                     by.interval = TRUE, time.window = interval,
                     names.window = stages)
K.stage <- env_kernel(env.data = W.stage, stages = stages, gaussian = TRUE)
names(K.stage$envCov)

## 5. Feeding a reaction-norm model
KE <- list(W = env_kernel(env.data = W.cov)$envCov)

Description

Reproduces Fig. 6c of Della Coletta et al. (2023): for a chosen PC, the loading of every index is plotted against its position in the growing season, with one loess curve per weather factor. This answers "which weather factors, at which point in the season, drive this PC".

Usage

env_loading_curve(
  x,
  pc = 1,
  factors = NULL,
  span = 0.4,
  points = TRUE,
  facet = FALSE,
  title = NULL,
  verbose = TRUE
)

Arguments

Details

Interpretation. A factor whose curve departs strongly from zero over a stretch of the season contributes heavily to that PC during that window. Sign is arbitrary up to the usual PCA reflection: only relative signs within a PC are meaningful.

Loess caveat. The smoother is cosmetic. With 51 intervals per factor the curve is well determined, but with few intervals it can imply structure that the underlying points do not support – keep points = TRUE to see the evidence.

Value

Invisibly, a data.frame of the plotted loadings (index, factor, day, loading), so the curves can be re-plotted or tabulated. When ggplot2 is available the ggplot object is attached as the "plot" attribute.

See Also

env_pca, env_pc_associate

Examples

pc <- env_pca(env_indices(maizeWTH, end.day = 130))

## All factors on one panel
env_loading_curve(pc, pc = 1)

## Temperature factors only, one panel each
env_loading_curve(pc, pc = 2,
                  factors = c("T2M", "T2M_MAX", "T2M_MIN"), facet = TRUE)

## Recover the plotted values
L <- env_loading_curve(pc, pc = 1)
head(L[order(-abs(L$loading)), ])

Description

Pearson correlation between each environmental PC and an external per-environment quantity – in the paper, module eigenmarkers; equally, a trait mean, heritability or predictive ability per environment. P-values are FDR-adjusted across PCs, as in the paper.

Usage

env_pc_associate(
  x,
  y,
  n.pc = NULL,
  method = "pearson",
  p.adjust = "BH",
  verbose = TRUE
)

Arguments

Details

Power warning. With $q$ environments each correlation has $q - 2$ degrees of freedom. At the paper's $q = 9$ even $r = 0.86$ gives an unadjusted $p \approx 0.003$, and a handful of PCs tested against a handful of modules exhausts the evidence quickly. Treat these associations as hypothesis- generating, not confirmatory, and prefer many environments.

FDR does not cover module selection. The adjustment applied here accounts for testing several PCs. It does not account for the eigenmarker being PC1 of a module that was itself chosen by clustering the same effect matrix. That selection inflates correlations and is invisible to any p-value computed from the observed data. See the worked example below for a permutation check, and for a demonstration that the permutation check is itself insufficient.

Value

A data.frame with column, PC, r, p, p.adj and n, sorted by adjusted p-value.

Examples

pc <- env_pca(env_indices(maizeWTH, end.day = 130))

## Trait mean per environment
ybar <- tapply(maizeYield$value, maizeYield$env, mean)
env_pc_associate(pc, ybar)



## ----------------------------------------------------------------------
## WORKED EXAMPLE: marker effect networks -> eigenmarkers -> environmental PCs
##
## The full analysis chain of Della Coletta et al. (2023), on simulated data
## where the answer is known in advance so the method can be checked.
##
## Planted truth:
##   module A markers respond to mid-season TEMPERATURE
##   module B markers respond to late-season PRECIPITATION
##   module C markers have real but environment-INDEPENDENT effects
## ----------------------------------------------------------------------

## --- 0. Simulate markers, weather and phenotypes ---------------------
set.seed(2026)
n.gid <- 160
n.marker <- 90
n.env <- 20
envs <- paste0("E", seq_len(n.env))
drivers <- scale(matrix(rnorm(n.env * 2), n.env, 2))
rownames(drivers) <- envs
M <- scale(matrix(rbinom(n.gid * n.marker, 2, 0.5), n.gid, n.marker))
dimnames(M) <- list(paste0("G", seq_len(n.gid)),
                    paste0("m", seq_len(n.marker)))
beta <- matrix(0.2, n.marker, n.env)
beta[1:30, ] <- rep(0.2 + 0.3 * drivers[, 1], each = 30)
beta[31:60, ] <- rep(0.2 + 0.3 * drivers[, 2], each = 30)
signal <- M %*% beta
Y <- signal + matrix(rnorm(length(signal), sd = 0.3), n.gid, n.env)
dimnames(Y) <- list(rownames(M), envs)

wth <- expand.grid(daysFromStart = 1:151, env = envs)
ie <- match(wth$env, envs)
day <- wth$daysFromStart
wth$T2M <- 22 + 4 * sin(pi * day / 151) +
  3 * drivers[ie, 1] * (day >= 50 & day <= 100) +
  rnorm(nrow(wth), sd = 0.5)
wth$PRECTOT <- exp(1 + 0.6 * drivers[ie, 2] * (day >= 100) +
                     rnorm(nrow(wth), sd = 0.1))

## --- 1. Marker effects per environment (RR-BLUP) ----------------------
## Solve via the n x n system: with p markers >> n genotypes this is far
## cheaper than the p x p form and gives the same ridge solution.
rrblup_effects <- function(y, Z, h2 = 0.5) {
  yc <- as.numeric(y) - mean(y)
  lambda <- ncol(Z) * (1 - h2) / h2
  a <- solve(tcrossprod(Z) + lambda * diag(nrow(Z)), yc)
  as.numeric(crossprod(Z, a))
}
beta.hat <- sapply(colnames(Y), function(e) rrblup_effects(Y[, e], M, h2 = 0.55))
dimnames(beta.hat) <- list(colnames(M), colnames(Y))

## --- 2. Marker effect network -> modules ------------------------------
## Markers whose effects do not VARY across environments cannot covary with
## anything; screening them out first avoids NA correlations and spurious
## modules. The paper used WGCNA; this is the same idea in base R.
eff.sd <- apply(beta.hat, 1, sd)
Bk     <- beta.hat[eff.sd > quantile(eff.sd, 0.35), , drop = FALSE]
adj    <- abs(cor(t(Bk)))^6                 # soft-thresholding power
mods   <- cutree(hclust(as.dist(1 - adj), method = "average"), h = 0.92)
mods   <- mods[mods %in% as.integer(names(table(mods))[table(mods) >= 15])]

## --- 3. Eigenmarkers --------------------------------------------------
## PC1 of a module's effect profiles, the analogue of a WGCNA module
## eigengene. Sign is arbitrary in PCA, so orient it deterministically --
## otherwise a module can flip between runs and reverse every downstream
## correlation.
eigenmarker <- function(B.mod) {
  X  <- scale(t(B.mod))
  em <- prcomp(X, center = FALSE, scale. = FALSE)$x[, 1]
  if (mean(cor(em, X)) < 0) em <- -em
  em
}
ids <- sort(unique(mods))
EM  <- sapply(ids, function(m) eigenmarker(Bk[names(mods)[mods == m], , drop = FALSE]))
dimnames(EM) <- list(colnames(beta.hat), paste0("ME", ids))

## --- 4. Environmental PCA and association -----------------------------
W  <- env_indices(wth, interval = 3, end.day = 151,
                  statistic.by = c(PRECTOT = "sum"))
pc <- env_pca(W)
assoc <- env_pc_associate(pc, EM, n.pc = 9, p.adjust = "BH")
head(assoc)

## --- 5. Name the drivers of an associated PC (paper Fig. 6c) ----------
k <- as.integer(sub("PC", "", assoc$PC[1]))
L <- data.frame(factor = pc$meta$factor, day = pc$meta$day,
                loading = pc$loadings[, k])
sort(tapply(abs(L$loading), L$factor, mean), decreasing = TRUE)
env_loading_curve(pc, pc = k)

## --- 6. Permutation null, and why it is NOT enough --------------------
## Permuting environment labels of the eigenmarker breaks the
## environment<->module link while preserving both structures.
perm.max <- replicate(500, {
  EMp <- EM[sample(nrow(EM)), , drop = FALSE]
  rownames(EMp) <- rownames(EM)
  max(abs(env_pc_associate(pc, EMp, n.pc = 9, verbose = FALSE)$r))
})
quantile(perm.max, 0.95)   # compare observed |r| with the permutation null

## A spurious module can pass both FDR adjustment and this permutation test.
## The permutation preserves the module definitions, so it cannot detect a
## module that was spuriously assembled in the first place. Guard against
## this with module stability across resampled genotypes, or by requiring
## an independent set of environments -- not with permutation alone.

Description

Runs prcomp on the environment x index matrix, as in Della Coletta et al. (2023), and returns scores, loadings and variance explained in a form the plotting functions in this module can consume.

Usage

env_pca(
  W,
  center = TRUE,
  scale. = TRUE,
  rank. = NULL,
  drop.constant = TRUE,
  verbose = TRUE
)

Arguments

Details

Scaling. The paper reports only that prcomp was used. Because the 17 factors are on incomparable scales – radiation in MJ/m2/day, precipitation in mm, FRUE unitless in [0,1] – an unscaled PCA is dominated by the highest-variance factor and the leading PCs largely describe that one variable. scale. = TRUE is therefore the default here. Set scale. = FALSE for a covariance-matrix PCA.

Rank. With $q$ environments and $k \gg q$ indices, at most $q - 1$ PCs carry variance. With the 9 environments of the paper that is 8, so the reported "ninth PC" explaining $<0.01\%$ is effectively numerical noise. Expect very few informative PCs unless you have many environments.

Value

A list of class "env_pca":

scores

environment x PC matrix (prcomp$x)

loadings

index x PC matrix (prcomp$rotation)

variance

data.frame with PC, sd, prop, cum

prcomp

the underlying prcomp object

dropped

indices removed as constant

meta

data.frame mapping each index to its factor, interval and day

See Also

env_indices, env_pca_biplot, env_loading_curve

Examples

data("maizeWTH")
W  <- env_indices(maizeWTH, end.day = 130)
pc <- env_pca(W)
head(pc$variance)
env_pca_scree(pc)
env_pca_biplot(pc)

Description

Scatter of environments on two PCs, reproducing Fig. 1c of Della Coletta et al. (2023), where each point is one growth environment labelled by its code.

Usage

env_pca_biplot(
  x,
  pc = c(1, 2),
  label = TRUE,
  group = NULL,
  title = "Environments in PC space",
  verbose = TRUE
)

Arguments

Value

A ggplot object when ggplot2 is available (invisibly), otherwise NULL after drawing a base plot.

See Also

env_pca, env_pca_scree

Examples

pc <- env_pca(env_indices(maizeWTH, end.day = 130))
env_pca_biplot(pc)
env_pca_biplot(pc, pc = c(1, 3))

Description

Plots the proportion of variance explained per PC, the visualisation the paper used to decide how many PCs to retain.

Usage

env_pca_scree(
  x,
  n.pc = NULL,
  cumulative = TRUE,
  title = "Variance explained by environmental PCs",
  verbose = TRUE
)

Arguments

Value

A ggplot object when ggplot2 is available (invisibly), otherwise NULL after drawing a base plot.

See Also

env_pca

Examples

pc <- env_pca(env_indices(maizeWTH, end.day = 130))
env_pca_scree(pc)

Description

Accumulates daily thermal time from planting and cuts each environment's season at crop-specific developmental thresholds, so that every downstream interval-based routine compares LIKE DEVELOPMENTAL STATES across sites rather than like calendar dates.

The returned stage column is a drop-in replacement for the interval produced by time.window in W_matrix(), summaryWTH() and env_typing().

Usage

env_phenology(
  env.data,
  env.id = "env",
  day.id = "daysFromStart",
  Tmax = "T2M_MAX",
  Tmin = "T2M_MIN",
  crop = NULL,
  stages = NULL,
  Tbase = NULL,
  Tupper = NULL,
  method = c("gdd", "baskerville", "trapezoid"),
  Tmaxdev = 45,
  photoperiod = FALSE,
  lat = "LAT",
  daylength.id = NULL,
  p.opt = NULL,
  p.sens = 0.02,
  day.type = NULL,
  planting.id = NULL,
  merge = TRUE,
  verbose = TRUE
)

Arguments

Details

Which method. "gdd" is the convention and what param_temperature() already computes; use it for comparability. "baskerville" is materially better where daily range straddles Tbase – high-elevation or early-spring sites – because the simple mean then understates accumulation. "trapezoid" is the only option that penalises supra-optimal heat, relevant for heat-stress work.

Photothermal time. Maize, sorghum and rice are quantitative short-day plants: long days delay flowering. With photoperiod = TRUE the daily increment is scaled by a linear daylength factor, which keeps tropical and temperate plantings on a common developmental scale. The default p.sens = 0.02 is deliberately mild; calibrate it before relying on the absolute stage dates.

Caveat. These thresholds are cultivar-generic. Maturity group alone moves maize R6 by 300-500 C d. Treat the stage boundaries as a consistent relative scale, not as a validated phenology prediction, unless you have calibrated stages against observed dates.

Value

A data.frame with these columns added:

dTT

daily thermal time (C d), or photothermal time if photoperiod = TRUE

cumTT

cumulative thermal time from planting

stage

ordered factor, the developmental stage of that day

dstage

relative position within the stage, in [0, 1)

stage_day

day counter within the stage (1-based)

The attribute "phenology" carries the stage table, the method and the cardinal temperatures. Environments whose season ends before the last threshold get NA for the stages never reached, and a warning names them – silently recycling the final stage would fabricate development that the weather does not support.

Author(s)

Germano Costa Neto

References

McMaster, G.S. & Wilhelm, W.W. (1997) Growing degree-days: one equation, two interpretations. Agric. For. Meteorol. 87(4), 291-300.

Baskerville, G.L. & Emin, P. (1969) Rapid estimation of heat accumulation from maximum and minimum temperatures. Ecology 50(3), 514-517.

Abendroth, L.J. et al. (2011) Corn Growth and Development. PMR 1009, Iowa State University Extension.

See Also

param_temperature, W_matrix, env_risk_profile, summaryWTH

Examples

data("maizeWTH")

## 1. Default maize staging on cumulative GDD
ph <- env_phenology(maizeWTH, crop = "maize")
table(ph$env, ph$stage)

## 2. Photothermal time, daylength taken from param_radiation()'s N column
ph2 <- env_phenology(processWTH(maizeWTH), crop = "maize",
                     photoperiod = TRUE, daylength.id = "N")

## 3. Custom stages and the Baskerville-Emin method
ph3 <- env_phenology(maizeWTH,
                     stages = c(veg = 500, flow = 1200, fill = 2200),
                     Tbase = 8, Tupper = 32, method = "baskerville")

## 4. Hand the stage factor to the existing covariable machinery.
##    This is the point of the function: W columns now align on
##    development instead of on calendar date.
W <- W_matrix(env.data = ph, var.id = c("T2M", "PRECTOT"),
              by.interval = TRUE, time.window = NULL)

Description

Quantifies how OFTEN each stress envirotype occurs at a site, by developmental stage, across many historical seasons. Where env_typing() describes the environments you actually observed, env_risk_profile() estimates the underlying climatology those environments were drawn from – the frequency, not the realisation.

This is the quantity a breeding programme needs to weight its testing network: a mega-environment that appears in 1 season out of 20 should not carry the same weight as one that appears in 12.

Usage

env_risk_profile(
  env.data,
  site.id = "env",
  date.id = "YYYYMMDD",
  var.id = NULL,
  years = NULL,
  planting = NULL,
  season.length = 150,
  cardinals = NULL,
  by.stage = TRUE,
  crop = "maize",
  stages = NULL,
  Tbase = NULL,
  Tupper = NULL,
  method = "gdd",
  min.days = 30L,
  return.seasons = FALSE,
  verbose = TRUE
)

Arguments

Details

Two frequencies are reported and they answer different questions. freq is the proportion of SEASONS in which the class occurred at least once – the probability of encountering that stress. mean_days is the expected NUMBER OF DAYS per season – the intensity. A site can have freq = 1.0 for mild heat while averaging only 3 days of it; that is a very different target from one averaging 40 days.

Return period is simply 1 / freq and inherits all the usual caveats: it is an average recurrence interval estimated from a finite record, not a schedule. With 20 years of data a return period above ~10 seasons is barely distinguishable from noise. The function will not print a return period longer than the record.

Stationarity. Frequencies computed over 1990-2024 describe that period. If the site is warming, they are already an underestimate of future heat risk. Use project_risk() with a scenario to see the shift.

Value

An object of class "env_risk": a list with

frequency

data.frame of site x stage x variable x class with n_years, n_occur, freq (probability a season shows that class at all), and mean_days (expected days per season)

severity

data.frame of site x stage x variable with the mean, sd and 10th/90th percentile of the variable across seasons

return_period

data.frame with 1 / freq, the expected number of seasons between occurrences of each stress class

seasons

per-site-year detail if return.seasons = TRUE

call, years, cardinals

provenance

Author(s)

Germano Costa Neto

References

Chapman, S.C. et al. (2000) Using crop simulation to generate genotype by environment interaction effects for sorghum in water-limited environments. Aust. J. Agric. Res. 51(2), 209-221.

Heinemann, A.B. et al. (2015) Drought impact on rainfed common bean production areas in Brazil. Agric. For. Meteorol. 225, 1-12.

See Also

env_phenology, env_typing, tpe_weights, project_risk

Examples

## 30 seasons of weather for three sites
sites <- data.frame(env = c("PIRA", "SETE", "NMAC"),
                    lat = c(-22.7, -19.4, -21.2),
                    lon = c(-47.6, -44.2, -45.2))

wth <- get_weather(env.id = sites$env, lat = sites$lat, lon = sites$lon,
                   start.day = "1994-01-01", end.day = "2024-12-31")
wth <- processWTH(wth)

## Risk of heat and drought at flowering, sown 15 October each year
rp <- env_risk_profile(wth, site.id = "env", planting = "10-15",
                       var.id = c("T2M_MAX", "PETP"),
                       crop = "maize", by.stage = TRUE)
rp

## How often is flowering heat-stressed at each site?
subset(rp$frequency, stage == "R1" & class == "heat_severe")

## Weights for a target population of environments
w <- tpe_weights(rp, stage = "R1", var = "PETP")

Description

Groups environments by their target metrics (Y: phenotypic means, heritabilities, selection accuracies, etc.) using K-means, then trains a Random Forest to predict those clusters from the environmental covariable matrix (W). The Random Forest variable importance ranks which environmental variables best explain the target-based grouping of environments. The number of clusters can be supplied by the user or optimised automatically with the average silhouette width.

Usage

env_target_importance(
  Y,
  W,
  k = NULL,
  k.max = 10,
  scale.Y = TRUE,
  ntree = 1000,
  importance.type = c("gini", "permutation"),
  nstart = 25,
  seed = NULL,
  save = FALSE,
  dir.path = NULL,
  verbose = TRUE
)

Arguments

Details

This is a supervised counterpart of env_cluster: instead of grouping environments by their enviromic profile, environments are grouped by what the breeder cares about (yield level, heritability, accuracy) and the enviromic data are then asked to explain that grouping. Variables at the top of importance are candidate drivers of the observed environmental stratification and good candidates for env_threshold_gxe.

Permutation importance ("permutation") is generally preferred when covariables differ in scale or in number of categories, since MeanDecreaseGini is biased towards high-cardinality predictors.

Value

A list with:

• clusters: named integer vector of the cluster of each environment.

• k: the number of clusters used.

• silhouette: data.frame of average silhouette width per tested k (NULL if k was given).

• kmeans: the kmeans object.

• randomForest: the fitted randomForest object.

• importance: data.frame of environmental variables ranked by importance (descending).

Author(s)

Germano Costa Neto

References

Costa-Neto, G., da Matta, D., Fernandes, I. K., & Heinemann, A. B. (2023). Environmental clusters defining breeding zones for tropical irrigated rice in Brazil. Agronomy Journal, 115(5). doi:10.1002/agj2.21481

See Also

env_cluster, W_matrix

Examples

if (requireNamespace("randomForest", quietly = TRUE)) {
  data("maizeYield"); data("maizeWTH")

  ## Environment-level target: mean grain yield per environment
  Y <- tapply(maizeYield$value, maizeYield$env, mean)

  ## Environmental covariables (environments x covariables)
  W <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
                var.id = c("T2M", "T2M_MAX", "PRECTOT"), statistic = "mean")

  ## Which environmental variables best explain the target-based grouping?
  out <- env_target_importance(Y = Y, W = W, seed = 1)
  out$k
  out$clusters
  head(out$importance)

  ## Permutation importance and a fixed number of clusters
  out2 <- env_target_importance(Y = Y, W = W, k = 2,
                                importance.type = "permutation", seed = 1)
  head(out2$importance)

  ## Multi-trait target (e.g. mean and heritability per environment)
  Ymat <- cbind(mean = Y, sd = tapply(maizeYield$value, maizeYield$env, stats::sd))
  out3 <- env_target_importance(Y = Ymat, W = W, seed = 1)
  out3$importance
}

Description

Returns environmental typologies (envirotypes) that can be used as envirotype markers. Typologies are given by cardinals (discrete intervals for each variable), by empirical quantiles, or – when envirotype_mining = TRUE – learned from the data: for each variable a K-means clustering is run over a grid of candidate k values and the optimal number of envirotypes is selected by the Calinski-Harabasz (variance ratio) criterion.

Usage

env_typing(
  env.data,
  var.id,
  env.id,
  cardinals = NULL,
  days.id = NULL,
  time.window = NULL,
  names.window = NULL,
  quantiles = NULL,
  id.names = NULL,
  by.interval = FALSE,
  scale = FALSE,
  format = NULL,
  ratio = FALSE,
  joint = FALSE,
  joint.by.interval = TRUE,
  combine.features = FALSE,
  combine.max.order = 2L,
  envirotype_mining = FALSE,
  k.range = 2:10,
  nstart = 25,
  iter.max = 100,
  seed = 1234,
  plot_envirotypes = FALSE,
  min.window = 2L,
  verbose = TRUE
)

Arguments

Details

The Calinski-Harabasz index for a partition into $k$ groups of $n$ observations is

$CH(k) = \frac{BGSS/(k-1)}{WGSS/(n-k)}$

where BGSS is the between-group and WGSS the within-group sum of squares. The $k$ maximising CH is retained. Clusters are relabelled in increasing order of their centres so that Envirotype_<var>_00001 is always the lowest range of the variable. For joint (multivariate) mining there is no natural ordering, so centroids are ordered along their first principal component and variables are standardised within the window before clustering so that no feature dominates the Euclidean distance.

Value

If envirotype_mining = FALSE (default) and plot_envirotypes = FALSE: a data.frame (or matrix, if format = 'wide') of envirotype frequencies.

Otherwise a list with:

typologies

data.frame (long) or matrix (wide) of envirotype frequencies.

envirotype_description

data.frame describing the numeric range captured by each envirotype: variable, label, cluster index, n, min, max, mean, sd, median and the lower/upper cut boundaries.

mining_summary

data.frame with the optimal k, the CH index at the optimum and the full CH profile over k.range.

feature_combinations

only when combine.features = TRUE: a data.frame mapping each generic combination id to its number of features (n_features), the features themselves and the time intervals.

plot

the plot object, or NULL.

Author(s)

Germano Costa Neto

See Also

W_matrix, env_kernel, T_matrix

Examples

data("maizeWTH")
env.data <- maizeWTH[maizeWTH$daysFromStart < 100, ]

## 1. Generic time intervals
env_typing(env.data = env.data, env.id = 'env', var.id = 'T2M', by.interval = TRUE)

## 2. Specific time intervals with names
env_typing(env.data = env.data, env.id = 'env', var.id = 'T2M',
           by.interval = TRUE,
           time.window  = c(0, 15, 35, 65, 90),
           names.window = c('1-initial growing', '2-leaf expansion I',
                            '3-leaf expansion II', '4-flowering', '5-grain filling'))

## 3. Cardinal (ecophysiological) thresholds
env_typing(env.data = env.data, env.id = 'env',
           var.id    = c('T2M', 'PRECTOT'),
           cardinals = list(T2M = c(0, 9, 22, 32, 45), PRECTOT = c(0, 5, 10, 25, 100)))

## 4. Data-driven univariate mining (K-means + Calinski-Harabasz)
out <- env_typing(env.data = env.data, env.id = 'env',
                  var.id = c('T2M', 'PRECTOT'),
                  envirotype_mining = TRUE, k.range = 2:8)
out$typologies
out$envirotype_description
out$mining_summary

## 5. Joint (multi-variable) envirotypes per development stage
outj <- env_typing(env.data = env.data, env.id = 'env',
                   var.id = c('T2M', 'PRECTOT'),
                   by.interval = TRUE, time.window = c(0, 30, 60, 90),
                   envirotype_mining = TRUE, joint = TRUE)
head(outj$envirotype_description)

## 6. Relative frequencies, wide format (ready for env_kernel)
ET <- env_typing(env.data = env.data, env.id = 'env', var.id = 'T2M',
                 format = 'wide', ratio = TRUE)

## 7. Envirotypes for feature combinations (singles + pairs by default)
outc <- env_typing(env.data = env.data, env.id = 'env',
                   var.id = c('T2M', 'PRECTOT', 'VPD'),
                   envirotype_mining = TRUE, combine.features = TRUE)
outc$feature_combinations        # generic id -> which variables
head(outc$envirotype_description)

## include triples too (combine.max.order = 3), or NULL for all orders
outc3 <- env_typing(env.data = env.data, env.id = 'env',
                    var.id = c('T2M', 'PRECTOT', 'VPD'),
                    envirotype_mining = TRUE, combine.features = TRUE,
                    combine.max.order = 3)

Description

Genomes to fields (G2F) is a multi-institutional, public collaborative to develop information and tools that support the translation of maize (Zea mays L.) genomic information into relevant phenotypes for the benefit of growers, consumers, and society. Building on existing maize genome sequence resources, the project focuses on developing approaches to improve phenomic predictability and facilitate the development and deployment of tools and resources that help address fundamental problems of sustainable agricultural productivity. Specific projects within G2F involve collaboration from research fields such as genetics, genomics, plant physiology, agronomy, climatology and crop modeling, computational sciences, statistics, and engineering. See https://bmcresnotes.biomedcentral.com/articles/10.1186/s13104-020-4922-8.

Here we used the geographic coordinates from 93 environments (diverse locations across 4 years)

Usage

data(G2F_2014_17)

Format

A data.frame containing environmental identifications (env), geographic coordinates (lat, lon) and the time window (start, end) given in YYYY-MM-DD.

data.frame file

Examples

require(EnvRtype)
data(G2F_2014_17)
head(G2F_2014_17)

# let's use the first 10 locations only
data  = G2F_2014_17
env   = data$env[1:10]
lon   = data$lon[1:10]
lat   = data$lat[1:10]
start = data$start[1:10]
end   = data$end[1:10]


# get_weather() queries the NASA POWER web service
weather.data =
  EnvRtype::get_weather(env.id = env,
                        lat = lat,
                        lon = lon,
                        start.day = start,
                        end.day = end)

head(weather.data)

Description

Extracts the FAO/IIASA Global Agro-Ecological Zone (GAEZ v4) class at a set of coordinates. By default it reads the GAEZ 57-class raster bundled with the package, so it works offline with no arguments beyond the coordinates. Companion to get_weather and get_soil.

Usage

get_AEZ(
  env.id,
  lat,
  lon,
  file = NULL,
  url = aez_default_url(),
  dir.path = NULL,
  legend = NULL,
  buffer = 0,
  overwrite = FALSE,
  timeout = 1800L,
  verbose = TRUE
)

Arguments

Details

Where the raster comes from. GAEZ v4 has no documented per-point API: it publishes whole-globe GeoTIFFs through a JavaScript portal. get_AEZ is therefore a raster-extract tool and resolves its raster in three steps – an explicit file, else the layer bundled under inst/extdata, else a download from url. The bundled layer is the GAEZ v4 57-class AEZ historical baseline (CRUTS32_Hist_8110, ~5 arc-min), which makes the default call fully offline and reproducible.

The default URL is unverified. It is only used when both file and the bundled raster are unavailable. FAO has reorganised GAEZ download paths between releases and the portal is JavaScript-rendered, so no stable path could be confirmed. If a download is triggered and fails:

1. open https://gaez.fao.org/pages/data-access-download,

2. download the AEZ classification layer manually,

3. pass the file with file = "path/to/layer.tif".

Coordinates must be WGS84 lon/lat. They are reprojected to the raster's CRS automatically, but only if the raster declares one. A raster with an undefined CRS triggers an error rather than a silent mis-extraction.

Interpreting buffer. A single-cell lookup at ~9 km resolution can be misleading near a zone boundary. With buffer > 0 the returned majority_share shows how homogeneous the neighbourhood is; a value near 1 means the site sits well inside one zone, while 0.4 means the classification is effectively arbitrary at that location.

Value

A data.frame with one row per environment:

env

environment id

lat, lon

coordinates as supplied

aez_code

raster cell value

aez_label, thermal, moisture, arable

from the legend, when one is used

majority_share

when buffer > 0, the proportion of cells in the neighbourhood holding the majority class – a purity measure

The raster path used is attached as the "source" attribute.

References

FAO & IIASA (2021). Global Agro-Ecological Zones (GAEZ v4). Rome, FAO. https://gaez.fao.org

See Also

aez_legend, get_weather, get_soil

Examples

sites <- data.frame(
  env = c("Ames_IA", "Lincoln_NE", "Piracicaba_BR"),
  lat = c(42.03, 40.81, -22.71),
  lon = c(-93.62, -96.68, -47.63))

## Default: uses the raster bundled with the package (offline)
aez <- get_AEZ(sites$env, sites$lat, sites$lon)
aez

## Check zone purity in a 25 km neighbourhood
aez25 <- get_AEZ(sites$env, sites$lat, sites$lon, buffer = 25000)
aez25[, c("env", "aez_label", "majority_share")]

## Raw codes, no labelling
get_AEZ(sites$env, sites$lat, sites$lon, legend = NA)


## Not run: 
## A different GAEZ layer you downloaded yourself
get_AEZ(sites$env, sites$lat, sites$lon, file = "GAEZ_AEZ_scenario.tif")

## End(Not run)

Description

Returns the 19 standard bioclimatic indices (BIO1-BIO19) for each site, from the 1970-2000 WorldClim 2.1 normals.

Usage

get_bioclim(
  env.id = NULL,
  lat = NULL,
  lon = NULL,
  vars = 1:19,
  resolution = c("2.5m", "5m", "10m", "30s"),
  dir.path = NULL,
  verbose = TRUE
)

Arguments

Value

data.frame with env, LAT, LON and one column per requested index. Units follow WorldClim: temperature indices in C, precipitation in mm, BIO3/BIO4/BIO15 dimensionless or percent.

What these are, and are not

Bioclim variables are LONG-TERM CLIMATE NORMALS, not the weather of your trial year. They describe the climate a site is drawn from, so they are appropriate for characterising a TPE, clustering locations, or as static site covariates. They are NOT a substitute for get_weather() + env_phenology() when you want to explain what happened in a particular season: a site's BIO1 is identical in a drought year and a wet year. Mixing normals and realised weather in one kernel without saying which is which is a common and serious error.

References

Fick, S.E. & Hijmans, R.J. (2017) WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37(12), 4302-4315. doi:10.1002/joc.5086

O'Donnell, M.S. & Ignizio, D.A. (2012) Bioclimatic predictors for supporting ecological applications in the conterminous United States. U.S. Geological Survey Data Series 691.

Examples

sites <- data.frame(env = c("PIRA", "SETE"),
                    lat = c(-22.7, -19.4), lon = c(-47.6, -44.2))
bc  <- get_bioclim(env.id = sites$env, lat = sites$lat, lon = sites$lon)
bc2 <- get_bioclim(sites$env, sites$lat, sites$lon, vars = c(1, 12))

Description

Retrieves downscaled CMIP6 monthly climate (WorldClim 2.1) for a set of sites and converts it into a DAILY weather series compatible with processWTH(), env_phenology() and water_balance(), by applying the monthly change signal to an observed baseline series (delta-change downscaling).

Usage

get_climate_scenario(
  env.id = NULL,
  lat = NULL,
  lon = NULL,
  baseline = NULL,
  model = "MPI-ESM1-2-HR",
  models = c("IPSL-CM6A-LR", "MPI-ESM1-2-HR", "MRI-ESM2-0", "UKESM1-0-LL"),
  scenario = "ssp245",
  period = "2041-2060",
  resolution = c("2.5m", "5m", "10m", "30s"),
  method = c("delta", "raw"),
  dir.path = NULL,
  baseline.id = "env",
  baseline.date = "YYYYMMDD",
  precip.scale = c("multiplicative", "additive"),
  verbose = TRUE
)

Arguments

Details

What delta-change does and does not do. The monthly mean change from the GCM is added to (temperature) or multiplied onto (precipitation) the observed daily series. The result therefore keeps the OBSERVED day-to-day sequence, wet-day frequency and spell structure, and shifts only the monthly means. This is the standard, defensible way to drive a crop model from monthly projections (Hawkins et al. 2013), but note the consequence: changes in variability, in dry-spell length, and in extreme-event frequency are NOT represented. A scenario that in reality brings the same total rain in fewer, heavier events will look identical to one that brings it evenly. Heat-stress day counts under delta-change are driven purely by the mean shift, so they are a lower bound.

If extremes matter to your question, you need daily GCM output and a quantile-mapping bias correction, not this function.

Model choice is not innocuous. Across CMIP6, equilibrium climate sensitivity ranges roughly 1.8-5.6 K, and several high-sensitivity models are now regarded as implausibly hot. A single-model result is one draw from a wide distribution. Prefer model = "ensemble" and report the spread; a projection whose sign flips across models is not a finding.

Network access. WorldClim is an external host. In restricted environments the download will fail; pre-populate dir.path with the GeoTIFFs and the function will use the cache.

Value

For method = "delta", a data.frame in get_weather() layout (env, YYYYMMDD, daysFromStart, T2M, T2M_MAX, T2M_MIN, PRECTOT, ...) with the scenario signal applied, plus scenario, model and period columns. For method = "raw", a monthly data.frame of projected tmin, tmax and prec per site. The attribute "deltas" holds the monthly change factors that were applied, so the perturbation is auditable.

Author(s)

Germano Costa Neto

References

Fick, S.E. & Hijmans, R.J. (2017) WorldClim 2. Int. J. Climatol. 37(12), 4302-4315.

Hawkins, E. et al. (2013) Calibration and bias correction of climate projections for crop modelling. Agric. For. Meteorol. 170, 19-31.

O'Neill, B.C. et al. (2016) The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geosci. Model Dev. 9, 3461-3482.

See Also

get_weather, project_risk, env_phenology

Examples

sites <- data.frame(env = c("PIRA", "SETE"),
                    lat = c(-22.7, -19.4), lon = c(-47.6, -44.2))

base <- get_weather(env.id = sites$env, lat = sites$lat, lon = sites$lon,
                    start.day = "2015-10-01", end.day = "2016-03-31")

## Mid-century, intermediate scenario, 4-model ensemble.
## The default 'models' set covers all four SSPs.
fut <- get_climate_scenario(env.id = sites$env, lat = sites$lat,
                            lon = sites$lon, baseline = base,
                            model = "ensemble", scenario = "ssp245",
                            period = "2041-2060")

## GFDL-ESM4 is fine for ssp126 or ssp370, just not ssp245/ssp585:
fut126 <- get_climate_scenario(env.id = sites$env, lat = sites$lat,
                               lon = sites$lon, baseline = base,
                               model = "GFDL-ESM4", scenario = "ssp126",
                               period = "2041-2060")

## Drop straight into the existing pipeline
fut <- processWTH(fut)
ph  <- env_phenology(fut, crop = "maize")

## Compare thermal time now vs mid-century
tapply(ph$dTT, ph$env, sum)

Description

Returns elevation (m) for each site. Elevation is a strong covariate for temperature and radiation, so it is often worth carrying alongside the weather-derived envirotypes – but note it is STATIC: it adds between-site structure, never within-season dynamics.

Usage

get_elevation(
  env.id = NULL,
  lat = NULL,
  lon = NULL,
  resolution = c("2.5m", "5m", "10m", "30s"),
  dir.path = NULL,
  verbose = TRUE
)

Arguments

Details

Sites in the sea or outside the land mask return NA with a warning naming them, rather than silently propagating NA into downstream kernels.

Value

data.frame with env, LAT, LON, ALT (m).

References

Fick, S.E. & Hijmans, R.J. (2017) WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37(12), 4302-4315. doi:10.1002/joc.5086

Examples

sites <- data.frame(env = c("PIRA", "SETE"),
                    lat = c(-22.7, -19.4), lon = c(-47.6, -44.2))
alt <- get_elevation(env.id = sites$env, lat = sites$lat, lon = sites$lon)

Description

Builds the genomic, envirotypic and genotype-by-environment kernels required by kernel_model. Given a genomic relationship matrix (p x p genotypes) and optionally an environmental relationship matrix (q x q environments), it returns one kernel per model term, already spectrally decomposed and ready to sample from.

Two computational engines are available and produce numerically identical results:

engine = "factored" (default)

Never allocates the n x n kernel. Each kernel is represented by its eigenvectors U (n x r) and eigenvalues s, obtained from the small p x p / q x q inputs. Memory is O(n*r) instead of O(n^2), which is what makes large multi-environment trials feasible.

engine = "dense"

Builds the explicit n x n matrix and attaches it as $Kernel. Slower and quadratic in memory, but needed if you want to inspect or post-process the kernel itself.

Set verbose = TRUE to trace each step with timings; verbose = FALSE (default) is entirely silent.

Usage

get_kernel(
  K_G = NULL,
  K_E = NULL,
  data = NULL,
  model = "RNMM",
  intercept.random = FALSE,
  engine = c("factored", "dense"),
  env = "env",
  gid = "gid",
  y = "value",
  tol = 1e-10,
  keep_var = 1,
  scale_kernels = TRUE,
  verbose = FALSE
)

Arguments

Value

A named list of kernels, class "kernel_set". Each element is a list with:

Type

"D" (dense) or "BD" (block-diagonal).

decomp

list of blocks, each list(s, U, nr, deltav, type, pos, scaled, tol).

diag_mean

mean(diag(K)), computed in O(n*r).

Kernel

the n x n matrix – only when engine = "dense".

Attributes: "ne" (observations per environment), "kd_meta" (decomposition settings, so kernel_model() can validate the cache), and "obs" (index vectors and level names).

Model structures

GE is a block-diagonal (within-environment) interaction; GxW is the reaction-norm interaction built from the environmental kernel.

Why the factored engine is exact

Every kernel here is a congruence transform of a small matrix. With ig mapping observations to genotypes and ie to environments,

$K_G = Z_g K_g Z_g' = K_g[ig, ig], \quad K_E = K_e[ie, ie]$

so the eigendecomposition of the n x n kernel follows from the small one in O(p^3) rather than O(n^3). Writing $K = A A'$ with $A = V[ig,]\sqrt{D}$, the nonzero eigenpairs come from the small Gram matrix $A'A$, giving an orthonormal U and exact eigenvalues.

The interaction kernel is an elementwise product of two expansions, which is the expansion of a Kronecker product: with $j = (ie-1) n_g + ig$, $K_{GE} = (K_e \otimes K_g)[j, j]$. On a complete balanced grid its eigenpairs are closed-form (eigenvalues $d_e \otimes d_g$, eigenvectors a Khatri-Rao product) so no large eigendecomposition is ever required.

Original Version

Costa-Neto et al (2021)

EnvRtype v.0.1.3, Sep 2026

Reusing kernels

Kernels depend only on genotypes and environments, never on the phenotype. Masking y for a cross-validation fold does not invalidate them, so build once and reuse across folds, chains and models.

Author(s)

Germano Costa Neto. Refactored version.

References

Costa-Neto G, Galli G, Carvalho HF, Crossa J, Fritsche-Neto R (2021). EnvRtype: a software to interplay enviromics and quantitative genomics in agriculture. G3, 11(4), jkab040.

See Also

kernel_model, decompose_kernels

Examples

data("maizeYield"); data("maizeG"); data("maizeWTH")

ECs <- W_matrix(env.data = maizeWTH[maizeWTH$daysFromStart < 100, ],
                var.id = c("FRUE", "PETP", "SRAD", "T2M_MAX"),
                statistic = "mean")
KE <- list(W = env_kernel(env.data = ECs)[[2]])
KG <- list(G = maizeG)

## factored engine (default): no n x n matrix is ever allocated
K <- get_kernel(K_G = KG, K_E = KE, data = maizeYield, model = "RNMM",
                env = "env", gid = "gid", y = "value", verbose = TRUE)

fit <- kernel_model(y = "value", data = maizeYield, random = K,
                    env = "env", gid = "gid",
                    iterations = 5000, burnin = 1000, verbose = TRUE)

## trim the eigenvalue tail: usually a large speedup, tiny accuracy cost
K99 <- get_kernel(K_G = KG, K_E = KE, data = maizeYield, model = "RNMM",
                  keep_var = 0.99)

## separate marker-predictable from residual genetic variance
Ki <- get_kernel(K_G = KG, K_E = KE, data = maizeYield, model = "RNMM",
                 intercept.random = TRUE)

Description

Queries the ISRIC SoilGrids REST API for soil properties at given coordinates.

Usage

get_soil(
  env.id = NULL,
  lat = NULL,
  lon = NULL,
  variables.names = "clay",
  depths = c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm"),
  stat = c("mean", "q_5", "q_50", "q_95", "uncertainty"),
  wide = TRUE,
  sleep = 12,
  tries = 3L,
  timeout = 60,
  verbose = TRUE
)

Arguments

Details

Accepted properties: bdod, cec, cfvo, clay, nitrogen, ocd, ocs, phh2o, sand, silt, soc, wv0010, wv0033, wv1500.

SoilGrids returns values in mapped units (e.g. clay in g/kg x 10). No unit conversion is applied; the unit column records what was served.

Value

A data.frame of soil features.

Author(s)

Brandon Monier, modified by Germano Costa Neto

References

Poggio, L. et al. (2021). SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. SOIL 7, 217-240.

See Also

get_weather, get_spatial

Examples

## Not run: 
## Requires the live SoilGrids service.
get_soil(env.id = "NM", lat = -13.05, lon = -56.05,
         variables.names = c("clay", "nitrogen"))

## End(Not run)

Description

Downloads SoilGrids properties one environment at a time with the same checkpointing scheme as get_weather_resumable: each result is written to disk immediately and recorded in an append-only log, so a crashed run resumes with restart_from_log.

SoilGrids is rate-limited and slower than NASA POWER, which makes long soil jobs the most frequent victims of a mid-run failure.

Usage

get_soil_resumable(
  env.id = NULL,
  lat = NULL,
  lon = NULL,
  variables.names = "clay",
  depths = c("0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm"),
  stat = c("mean", "q_5", "q_50", "q_95", "uncertainty"),
  dir.path = NULL,
  run.id = "soil_run",
  sleep = 20,
  tries = 3L,
  resume = TRUE,
  verbose = TRUE
)

Arguments

Details

Results are stored long on disk, one row per property-depth, and reshaped to wide only at assembly. Storing long is what makes resumption safe: parts fetched at different times still stack correctly, whereas pre-widened parts could have mismatched columns.

See get_weather_resumable for the full checkpoint design.

Value

A data.frame of soil properties (wide by default), with attributes run.dir, n_ok, n_failed and failed.

See Also

restart_from_log, read_progress_log, get_soil

Examples

## Not run: 
## Long-running download job requiring access to SoilGrids.
grid <- data.frame(env = paste0("P", 1:200),
                   lat = runif(200, -13, -12),
                   lon = runif(200, -56, -55))

soil <- get_soil_resumable(
  env.id = grid$env, lat = grid$lat, lon = grid$lon,
  variables.names = c("clay", "sand", "silt", "soc", "phh2o"),
  run.id = "zoning_grid")

## After a crash, in a new session:
soil <- restart_from_log(run.id = "zoning_grid")

## Points in the ocean never resolve -- keep the successes
soil <- restart_from_log(run.id = "zoning_grid", retry.failed = FALSE)

## End(Not run)

Description

Extracts point values from a raster at the coordinates of each environment. This function supersedes extract_GIS(), which relied on the retired raster and sp packages.

Usage

get_spatial(
  digital.raster = NULL,
  which.raster.number = NULL,
  env.dataframe = NULL,
  env.id = NULL,
  lat = NULL,
  lng = NULL,
  name.feature = NULL,
  merge = TRUE,
  crs = 4326,
  verbose = TRUE
)

Arguments

Details

Coordinates are validated before extraction: latitude must lie in [-90, 90] and longitude in [-180, 180]. This catches the common lat/lon transposition, which otherwise yields silent NAs or values from the wrong hemisphere.