diff --git a/NAMESPACE b/NAMESPACE
index 34840ccf..5b26743b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -45,8 +45,9 @@ export(listw2sn, sn2listw, read.gwt2nb,	write.sn2gwt,
         read.swmdbf2listw, read_swm_dbf, write.swmdbf, write_swm_dbf, write.sn2DBF,
         lm.LMtests, lm.RStests,
 	lm.morantest, localG, localG_perm, localmoran, localmoran_perm, moran,
-	moran.test, moran.mc, moran.plot, localmoran.sad, lm.morantest.sad,
-	nb2listw, nb2listwdist, nb2mat, listw2mat, mat2listw, nbdists, nblag,
+	moran.test, moran.mc, moran.plot, moran.plot.drop, moran.plot.seismogram,
+	localmoran.sad, lm.morantest.sad, nb2listw, nb2listwdist, nb2mat,
+	listw2mat, mat2listw, nbdists, nblag,
 	nblag_cumul, poly2nb, read.gal, write.nb.gal, read.geoda, relativeneigh,
 	soi.graph, sp.correlogram, sp.mantel.mc, set.spChkOption, chkIDs,
 	get.spChkOption, spNamedVec, tri2nb,
diff --git a/R/moran.plot.drop.R b/R/moran.plot.drop.R
new file mode 100644
index 00000000..57eac58c
--- /dev/null
+++ b/R/moran.plot.drop.R
@@ -0,0 +1,133 @@
+moran.plot.drop <- function(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL, significant = TRUE, xlab = NULL, ylab = NULL, return_df = TRUE, spChk = NULL, labels = NULL, zero.policy=attr(listw, "zero.policy")) {
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  if (!inherits(locmoran, "localmoran")) 
+    stop(paste(deparse(substitute(locmoran)), "is not a localmoran object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(significant))
+  stopifnot(is.logical(return_df))
+  stopifnot(is.numeric(alpha))
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x)) 
+    stop(paste(xname, "is not a numeric vector"))
+  if (anyNA(x)) 
+    stop("NA in X")
+  n <- length(listw$neighbours)
+  if (n != length(x)) 
+    stop("objects of different length")
+  usePadj <- FALSE
+  if (!is.null(adjusted_p)) {
+    if(n != length(adjusted_p))
+      stop("objects of different length")
+    if (anyNA(adjusted_p)) 
+      stop("NA in vector of adjusted p values")
+    usePadj <- TRUE
+  }
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(locmoran, listw)) 
+    stop("Check of data and weights ID integrity failed")
+  labs <- TRUE
+  if (is.logical(labels)) {
+    if(!labels)
+      labs <- FALSE
+    labels <- as.character(attr(listw, "region.id"))
+  } else if (!is.logical(labels) && !is.null(labels)) {
+    if(length(labels) != n) {
+      warning("Length of the labels vector does not match number of regions. region.id from the listw object is used instead.")
+      labels <- as.character(attr(listw, "region.id"))
+    }
+  } else if (is.null(labels)) {
+    labs <- FALSE
+    labels <- as.character(attr(listw, "region.id"))
+  }
+  if (is.null(xlab))
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged centred", xname)
+  
+  WX <- lag.listw(listw, scale(x, scale = F), zero.policy = zero.policy)
+  if (anyNA(WX)) warning("no-neighbour observation(s) found - use zero.policy=TRUE")
+  if(usePadj)
+    cv <- min(abs(locmoran[(which(adjusted_p <= alpha/2)), 4]))
+  else
+    cv <- abs(qnorm(1 - alpha, lower.tail = FALSE))
+
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))] <- b2[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))]
+  
+  if(significant) {
+    x_q1 <- x[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    y_q1 <- WX[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    b_q1 <- b[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    labels_q1 <- labels[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] >= cv)]
+    x_q2 <- x[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    y_q2 <- WX[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    b_q2 <- b[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    labels_q2 <- labels[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    x_q3 <- x[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    y_q3 <- WX[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    b_q3 <- b[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    labels_q3 <- labels[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] >= cv)]
+    x_q4 <- x[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    y_q4 <- WX[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    b_q4 <- b[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+    labels_q4 <- labels[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] <= (-1) * cv)]
+  } else {
+    x_q1 <- x[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    y_q1 <- WX[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    b_q1 <- b[which(x > mean(x) & WX > mean(WX) & locmoran[, 4] < cv)]
+    x_q2 <- x[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    y_q2 <- WX[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    b_q2 <- b[which(x > mean(x) & WX < mean(WX) & locmoran[, 4] > (-1) * cv)]
+    x_q3 <- x[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    y_q3 <- WX[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    b_q3 <- b[which(x < mean(x) & WX < mean(WX) & locmoran[, 4] < cv)]
+    x_q4 <- x[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+    y_q4 <- WX[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+    b_q4 <- b[which(x < mean(x) & WX > mean(WX) & locmoran[, 4] > (-1) * cv)]
+  }
+  
+  lw.lm <- lm(WX ~ x)
+  plot(x, WX, xlab=xlab, ylab=ylab, pch = 20, cex = 0.5, col = "gray70", xlim = c(min(x),max(x)), ylim = c(min(WX, b),max(WX, b)))
+  abline(h = mean(WX), lty = "dashed", col = "grey30")
+  abline(v = mean(x), lty = "dashed", col = "grey30")
+  abline(lw.lm, lty = "dotted", col = "grey40")
+  
+  for(i in 1:length(x_q1)) {
+    lines(c(x_q1[i], x_q1[i]), c(y_q1[i], b_q1[i]), type = "l", col = "lightpink", lwd = 1)
+    points(c(x_q1[i]), c(y_q1[i]), pch = 20, cex = 0.5, col = "firebrick")
+    if (labs) 
+      text(x_q1[i], y_q1[i], labels = labels_q1[i], pos = 2, cex = 0.5, col = "firebrick")
+  }
+
+  for(i in 1:length(x_q2)) {
+    lines(c(x_q2[i], x_q2[i]), c(y_q2[i], b_q2[i]), type = "l", col = "lightblue1", lwd = 1)
+    points(c(x_q2[i]), c(y_q2[i]), pch = 20, cex = 0.5, col = "royalblue")
+    if (labs) 
+      text(x_q2[i], y_q2[i], labels = labels_q2[i], pos = 2, cex = 0.5, col = "royalblue")
+  }
+      
+  for(i in 1:length(x_q3)) {
+    lines(c(x_q3[i], x_q3[i]), c(y_q3[i], b_q3[i]), type = "l", col = "lightpink", lwd = 1)
+    points(c(x_q3[i]), c(y_q3[i]), pch = 20, cex = 0.5, col = "firebrick")
+    if (labs) 
+      text(x_q3[i], y_q3[i], labels = labels_q3[i], pos = 2, cex = 0.5, col = "firebrick")
+  }
+    
+  for(i in 1:length(x_q4)) {
+    lines(c(x_q4[i], x_q4[i]), c(y_q4[i], b_q4[i]), type = "l", col = "lightblue1", lwd = 1)
+    points(c(x_q4[i]), c(y_q4[i]), pch = 20, cex = 0.5, col = "royalblue")
+    if (labs) 
+      text(x_q4[i], y_q4[i], labels = labels_q4[i], pos = 2, cex = 0.5, col = "royalblue")
+  }
+  
+  if(return_df) {
+    res <- data.frame(labels = labels, x = x, WX = WX, b = b, line_lengths = abs(WX - b))
+    invisible(res)
+  }
+}
diff --git a/R/moran.plot.seismogram.R b/R/moran.plot.seismogram.R
new file mode 100644
index 00000000..da83256a
--- /dev/null
+++ b/R/moran.plot.seismogram.R
@@ -0,0 +1,91 @@
+moran.plot.seismogram <- function(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL, xlab = NULL, ylab = NULL, return_df = TRUE, spChk = NULL, zero.policy=attr(listw, "zero.policy")) {
+  if (!inherits(listw, "listw")) 
+    stop(paste(deparse(substitute(listw)), "is not a listw object"))
+  if (!inherits(locmoran, "localmoran")) 
+    stop(paste(deparse(substitute(locmoran)), "is not a localmoran object"))
+  stopifnot(is.vector(x))
+  stopifnot(is.logical(return_df))
+  stopifnot(is.numeric(alpha))
+  if (is.null(zero.policy))
+    zero.policy <- get.ZeroPolicyOption()
+  stopifnot(is.logical(zero.policy))
+  xname <- deparse(substitute(x))
+  if (!is.numeric(x)) 
+    stop(paste(xname, "is not a numeric vector"))
+  if (any(is.na(x))) 
+    stop("NA in X")
+  n <- length(listw$neighbours)
+  if (n != length(x)) 
+    stop("objects of different length")
+  usePadj <- FALSE
+  if (!is.null(adjusted_p)) {
+    if(n != length(adjusted_p))
+      stop("objects of different length")
+    if (any(is.na(adjusted_p))) 
+      stop("NA in vector of adjusted p values")
+    usePadj <- TRUE
+  }
+  if (is.null(spChk)) 
+    spChk <- get.spChkOption()
+  if (spChk && !chkIDs(locmoran, listw)) 
+    stop("Check of locmoran and weights ID integrity failed")
+  if (is.null(xlab)) 
+    xlab <- xname
+  if (is.null(ylab)) 
+    ylab <- paste("spatially lagged centred", xname)
+
+  WX <- lag.listw(listw, scale(x, scale = F), zero.policy = zero.policy)
+  if (anyNA(WX)) warning("no-neighbour observation(s) found - use zero.policy=TRUE")
+  if(usePadj)
+    cv <- min(abs(locmoran[(which(adjusted_p <= alpha/2)),4]))
+  else
+    cv <- abs(qnorm(1 - alpha, lower.tail = FALSE))
+  b <- ((cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b2 <- ((-cv * sqrt(locmoran[, 3])) + locmoran[, 2]) * var(x) / (x - mean(x))
+  b[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))] <- b2[which((x < mean(x) & WX > mean(WX)) | (x > mean(x) & WX < mean(WX)))]
+  
+  x_q1 <- x[which(x > mean(x) & WX > mean(WX))]
+  y_q1 <- WX[which(x > mean(x) & WX > mean(WX))]
+  b_q1 <- b[which(x > mean(x) & WX > mean(WX))]
+  x_q2 <- x[which(x > mean(x) & WX < mean(WX))]
+  y_q2 <- WX[which(x > mean(x) & WX < mean(WX))]
+  b_q2 <- b[which(x > mean(x) & WX < mean(WX))]
+  x_q3 <- x[which(x < mean(x) & WX < mean(WX))]
+  y_q3 <- WX[which(x < mean(x) & WX < mean(WX))]
+  b_q3 <- b[which(x < mean(x) & WX < mean(WX))]
+  x_q4 <- x[which(x < mean(x) & WX > mean(WX))]
+  y_q4 <- WX[which(x < mean(x) & WX > mean(WX))]
+  b_q4 <- b[which(x < mean(x) & WX > mean(WX))]
+
+  lw.lm <- lm(WX ~ x)
+  plot(x, WX, xlab=xlab, ylab=ylab, pch = 20, cex = 0.5, col = "gray70", xlim = c(min(x),max(x)), ylim = c(min(WX, b),max(WX, b)))
+  abline(h = mean(WX), lty = "dashed", col = "grey30")
+  abline(v = mean(x), lty = "dashed", col = "grey30")
+  abline(lw.lm, lty = "dotted", col = "grey40")
+    
+  df_q1 <- data.frame(x_q1, b_q1)
+  df_q1 <- df_q1[order(x_q1),]
+  df_q2 <- data.frame(x_q2, b_q2)
+  df_q2 <- df_q2[order(x_q2),]
+  df_q3 <- data.frame(x_q3, b_q3)
+  df_q3 <- df_q3[order(x_q3),]
+  df_q4 <- data.frame(x_q4, b_q4)
+  df_q4 <- df_q4[order(x_q4),]
+
+  for(i in 1:(length(df_q1$x_q1) - 1))
+    lines(c(df_q1$x_q1[i], df_q1$x_q1[i+1]), c(df_q1$b_q1[i], df_q1$b_q1[i+1]), type = "l", col = "firebrick")
+  
+  for(i in 1:length(df_q2$x_q2)) 
+    lines(c(df_q2$x_q2[i], df_q2$x_q2[i+1]), c(df_q2$b_q2[i], df_q2$b_q2[i+1]), type = "l", col = "royalblue")
+        
+  for(i in 1:length(df_q3$x_q3))
+    lines(c(df_q3$x_q3[i], df_q3$x_q3[i+1]), c(df_q3$b_q3[i], df_q3$b_q3[i+1]), type = "l", col = "firebrick")
+    
+  for(i in 1:length(df_q4$x_q4))
+    lines(c(df_q4$x_q4[i], df_q4$x_q4[i+1]), c(df_q4$b_q4[i], df_q4$b_q4[i+1]), type = "l", col = "royalblue")
+
+  if(return_df) {
+    res <- data.frame(labels=as.character(attr(listw, "region.id")), x=x, wx=WX, b=b)
+    invisible(res)
+  }
+}
diff --git a/man/moran.plot.Rd b/man/moran.plot.Rd
index e64141a4..05a01e35 100644
--- a/man/moran.plot.Rd
+++ b/man/moran.plot.Rd
@@ -19,7 +19,7 @@ moran.plot(x, listw, y=NULL, zero.policy=attr(listw, "zero.policy"), spChk=NULL,
   \item{spChk}{should the data vector names be checked against the spatial objects for identity integrity, TRUE, or FALSE, default NULL to use \code{get.spChkOption()}}
   \item{labels}{character labels for points with high influence measures, if set to FALSE, no labels are plotted for points with large influence}
   \item{xlab}{label for x axis}
-  \item{ylab}{label for x axis}
+  \item{ylab}{label for y axis}
   \item{quiet}{default NULL, use !verbose global option value; if TRUE, output of summary of influence object suppressed}
   \item{plot}{default TRUE, if false, plotting is suppressed}
   \item{return_df}{default TRUE, invisibly return a data.frame object; if FALSE invisibly return an influence measures object}
diff --git a/man/moran.plot.drop.Rd b/man/moran.plot.drop.Rd
new file mode 100644
index 00000000..527100d8
--- /dev/null
+++ b/man/moran.plot.drop.Rd
@@ -0,0 +1,57 @@
+\name{moran.plot.drop}
+\alias{moran.plot.drop}
+
+\title{Moran drop plot}
+
+\description{
+A version of the Moran scatterplot (see \code{\link{moran.plot}}) supplemented by lines indicating \emph{p} values for visual inspection of statistical significance.
+}
+\usage{
+moran.plot.drop(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL,
+  significant = TRUE, xlab = NULL, ylab = NULL, return_df = TRUE,
+  spChk = NULL, labels = NULL, zero.policy=attr(listw, "zero.policy"))
+}
+
+\arguments{
+  \item{x}{a numerical vector holding the attribute of interest}
+  \item{listw}{a \code{listw} spatial weights object}
+  \item{locmoran}{a fitted object of type \code{localmoran}}
+  \item{alpha}{default 0.05; the desired significance level regarding local Moran's \emph{I}}
+  \item{adjusted_p}{default NULL; an optional vector of \emph{p} values adjusted to account for multiple testing as is returned by \code{\link{p.adjustSP}}; if NULL, standard normal approximation is used to determine critical values}
+  \item{significant}{default TRUE; a parameter indicating whether to display critical value distances of significant (default) or non-significant observations}
+  \item{xlab}{default NULL; an optional label for the x-axis}
+  \item{ylab}{default NULL; an optional label for the y-axis}
+  \item{return_df}{default TRUE; invisibly return a data.frame object, if FALSE do not return anything}
+  \item{spChk}{default NULL to use \code{get.spChkOption()}; should the locmoran names be checked against the listw spatial objects for identity integrity, TRUE, or FALSE}
+  \item{labels}{default NULL; no labels are plotted by default; region IDs from the \code{listw} object are used as labels of significant observations if set to TRUE; custom labels are used if a character vector is provided}
+  \item{zero.policy}{default option stored in the listw object; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors}
+}
+
+\details{
+The Moran drop plot is a version of the Moran scatterplot supplemented by visual indications of \emph{p} values. The standard Moran scatterplot provides information about the effect size but not about the level of confidence to determine whether the effects shown might be more than just random outcomes. The Moran drop plot marks significant points in red (positive) and blue (negative spatial autocorrelation) and adds so-called drop lines connecting the significant observations to the scatterplot positions of their associated critical values. The coordinates of the latter are determined either under the assumption of approximate standard normality of the z-scores of the local Moran's \emph{I} values, or based on \emph{p} values provided through \code{adjusted_p}. In the latter case, the critical value is approximated based on the observation with the highest adjusted \emph{p} value that still satisfies the selected significance level. The longer the lines, the lower the associated \emph{p} values. The visualisation thus enables visual inspection of statistical significance while maintaining the relationship to both attribute value ranges and scatterplot quadrants. It is also possible to invert the visualised relationship and display the distances of non-significant observations to their corresponding critical values (if significant is set to FALSE).
+}
+
+\value{
+When return_df is TRUE, a data frame object with the following members is returned:
+  \item{labels}{either the labels provided or the region ids, if not specified}
+  \item{x}{the attribute values}
+  \item{wx}{the spatial lags of the centred attribute values}
+  \item{b}{the y-coordinates (i.e. hypothetical spatial lags) of the critical values given x}
+  \item{line_lengths}{the absolute distances between b and x}
+}
+
+\references{Westerholt, R. (2024): Extending the Moran scatterplot by indications of critical values and \emph{p}-values: introducing the Moran seismogram and the drop plot. Proceedings of the 32nd Annual GIS Research UK Conference (GISRUK), Leeds, UK. \url{https://doi.org/10.5281/zenodo.10897792}}
+
+\author{Rene Westerholt \email{rene.westerholt@tu-dortmund.de}}
+
+\seealso{\code{\link{moran.plot}}}
+
+\examples{
+# Boston example (CMEDV; owner-occupied housing in USD)
+data(boston)
+boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1])
+boston.nb <- poly2nb(boston.tr)
+boston.listw <- nb2listw(boston.nb)
+moran.plot.drop(boston.c$CMEDV, boston.listw, localmoran(boston.c$CMEDV, boston.listw), 0.01,
+  significant = TRUE, labels = NULL)
+}
diff --git a/man/moran.plot.seismogram.Rd b/man/moran.plot.seismogram.Rd
new file mode 100644
index 00000000..605837e3
--- /dev/null
+++ b/man/moran.plot.seismogram.Rd
@@ -0,0 +1,52 @@
+\name{moran.plot.seismogram}
+\alias{moran.plot.seismogram}
+
+\title{Moran seismogram}
+
+\description{
+A variant of the Moran scatterplot (see \code{\link{moran.plot}}) supplemented by lines connecting location-wise critical value configurations of attribute values and spatial lags. The plot allows for visual inspection of potential spatial weights misspecifiation.
+}
+\usage{
+moran.plot.seismogram(x, listw, locmoran, alpha = 0.05, adjusted_p = NULL,
+  xlab = NULL, ylab = NULL, return_df = TRUE, spChk = NULL, zero.policy = attr(listw, "zero.policy"))
+}
+
+\arguments{
+  \item{x}{a numerical vector holding the attribute of interest}
+  \item{listw}{a \code{listw} spatial weights object}
+  \item{locmoran}{a fitted object of type localmoran}
+  \item{alpha}{default 0.05; the desired significance level regarding local Moran's \emph{I}}
+  \item{adjusted_p}{default NULL; an optional vector of \emph{p} values adjusted to account for multiple testing as is returned by \code{\link{p.adjustSP}}; if NULL, standard normal approximation is used to determine critical values}
+  \item{xlab}{default NULL; an optional label for the x-axis}
+  \item{ylab}{default NULL; an optional label for the y-axis}
+  \item{return_df}{default TRUE; invisibly return a data.frame object, if FALSE do not return anything}
+  \item{spChk}{default NULL to use \code{get.spChkOption()}; should the locmoran names be checked against the listw spatial objects for identity integrity, TRUE, or FALSE}
+  \item{zero.policy}{default option stored in the listw object; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors}
+}
+
+\details{
+The Moran seismogram is a version of the Moran scatterplot that is complemented by lines connecting the location-specific critical values of local Moran's \emph{I} in the plot. The y-coordinates associated with the critical value are determined either under the assumption of approximate standard normality of the z-scores of the local Moran's \emph{I} values, or based on \emph{p} values provided through \code{adjusted_p}. In the latter case, the critical value is approximated based on the observation with the highest adjusted \emph{p} value that still satisfies the selected significance level. Lines in quadrants with positive spatial autocorrelation are shown in red and lines in quadrants with negative spatial autocorrelation are shown in blue. The representation is similar to a seismogram for detecting earthquakes and thus reveals potentially suspicious local spatial weights configurations by visualising spikes. The latter are displayed in an integrated manner with their positions in the attribute value range and in connection with the types of the associated spatial patterns (by the quadrants of the scatterplot).
+}
+
+\value{
+When return_df is TRUE, a data frame object with the following members is returned:
+  \item{labels}{either the labels provided or the region ids, if not specified}
+  \item{x}{the attribute values}
+  \item{wx}{the spatial lags of the centred attribute values}
+  \item{b}{the y-coordinates (i.e. hypothetical spatial lags) of the critical values given x}
+}
+
+\references{Westerholt, R. (2024): Extending the Moran scatterplot by indications of critical values and \emph{p}-values: introducing the Moran seismogram and the drop plot. Proceedings of the 32nd Annual GIS Research UK Conference (GISRUK), Leeds, UK. \url{https://doi.org/10.5281/zenodo.10897792}}
+
+\author{Rene Westerholt \email{rene.westerholt@tu-dortmund.de}}
+
+\seealso{\code{\link{moran.plot}}}
+
+\examples{
+# Boston example (CMEDV; owner-occupied housing in USD)
+data(boston)
+boston.tr <- sf::st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1])
+boston.nb <- poly2nb(boston.tr)
+boston.listw <- nb2listw(boston.nb)
+moran.plot.seismogram(boston.c$CMEDV, boston.listw, localmoran(boston.c$CMEDV, boston.listw), 0.01, zero.policy = TRUE)
+}