Analysis of Spatial Coexistence and Neighbor Relationships in Forest Using Planar Graph Model
Longxiao
2021-04-21
In the study of organisms that live in a sessile manner (such as various plants, corals, and sedentary mollusks), direct observations of their interaction behaviors are difficult. Therefore, the interaction relationships are often inferred from their distribution or spatial coexistence. For example, it is generally inferred that species exhibiting low co-occurrence tend to have strong competitive interactions, whereas those with frequent co-occurrence are likely to experience weak competition or engage in mutualistic interactions. In the case of woody plants, spatial coexistence typically considers the neighbor relationships between individual plants. However, previous studies have generally described the neighbors of an individual as all other individuals within a certain spatial range (Hubbell et al., 2001; Peters, 2003; Stoll and Newbery, 2005; Comita et al., 2010). This method is computationally simple and requires fewer resources, but its definition of neighbors is vague, and the selection of range is often based on experience rather than a rigorous mathematical foundation.
To better define and analyze the spatial coexistence relationships of woody plants, we developed a method based on planar graphs to analyze neighbor relationships.
#packages
library(ggplot2)
library(pheatmap)
library(dplyr)
library(deldir)1 Defining Neighbor Relationships in a Planar Graph
Consider a set of points in a finite plane, which represent the individuals in forest.
# generate example: poisson distributed random points in 1*1 plane, set density as \rho.
example_generator<-function(rho = 100, X = 1, Y = 1, seed = 123){
set.seed(seed)
# possion distributed points
n_poisson <- rpois(1, lambda = rho * X * Y)
poisson_points <- data.frame(
x = runif(n_poisson, 0, X),
y = runif(n_poisson, 0, Y),
type = "sample"
)
return(poisson_points)
}rho = 100
X = 1
Y = 1
seed = 123
tree_sample <- example_generator()
#visualization
p1 <- ggplot(tree_sample, aes(x, y)) +
geom_point(color = "blue", alpha = 0.6) +
coord_fixed() +
ggtitle("simulated tree individuals in a spot")+
theme_minimal()
p1
We propose the following model: each individual is associated with a neighborhood centered on itself, with neighborhoods of different individuals being non-overlapping. Furthermore, the union of all these neighborhoods completely covers the sampling area. In other words, these neighborhoods form a partition of the sampling region.
Neighbor relationships are defined by the spatial proximity of the neighborhoods. Specifically, if the neighborhoods of two points are adjacent, these two points are considered neighbors. In this model, the neighbor relationship can be represented as a set of edges in the plane, where these edges do not intersect. In other words, the total neighbor relationships form a planar graph.
2 Complete Occupation in Subtropical Forests
Considering that subtropical forests often have a high canopy density, we may assume that the occupation of the plane by the set of points is sufficiently complete. From this, we infer that the number of neighbor relationships between points should be maximized, ensuring that there are no disconnected edges on the plane that would meet the required conditions. This leads to the generation of a complete planar graph.
3 Minimizing the Spatial Distance Between Neighboring Points
Since we are focusing on spatial neighbor relationships, the goal is to find the nearest set of neighboring points such that the spatial distance between neighbors is minimized. Additionally, for any given point, its neighbors must form a polygon containing only that point. This necessitates the creation of a minimum complete planar graph.
4 Mathematical Constraints for the Planar Graph \(G(V, E)\)
As state above, the graph \(G(V, E)\) must satisfy the following mathematical constraints:
Non-intersecting Edges: \[ \forall e_i, e_j \in E, e_i \cap e_j = \emptyset \] This means that the edges \(e_i\) and \(e_j\) should not intersect.
Minimization of Edge Length: \[ \min \sum e_i \quad \text{for} \quad e_i \in E \] This minimizes the total length of the edges.
Maximization of Vertex Degree: \[ \max \sum \text{degree}(v_i) \quad \text{for} \quad v_i \in V \] This maximizes the degree of the vertices (i.e., the number of edges connected to a vertex).
Where: - \(V\) and \(E\) represent the set of vertices and edges of the graph \(G\), - \(\text{degree}(v_i)\) denotes the degree of the vertex \(v_i\), which is the number of edges connected to that vertex.
5 The Algorithm Used to Implement the Model
We implement this process in real-world woody plant distribution data (here illustrated with the example generated above) using a global search approach. The process begins by generating a set of boundary points that define the sampling plot boundaries. This avoids incorrect neighbor relationships that might be formed due to boundary constraints.
boundary_generator <- function(rho = 100, X = 1, Y = 1){
# average distance scale
avg_dist <- sqrt(X * Y / rho)
# boundary area
x_seq <- seq(-avg_dist, X + avg_dist, by = avg_dist)
y_seq <- seq(-avg_dist, Y + avg_dist, by = avg_dist)
boundary_x <- expand.grid(
x = x_seq, y = c(-avg_dist, Y + avg_dist))
boundary_y <- expand.grid(
x = c(-avg_dist, X + avg_dist), y = y_seq)
boundary_points <- rbind(boundary_x, boundary_y)
boundary_points <- unique(boundary_points)
boundary_points$type <- "boundary"
return(boundary_points)
}boundary <- boundary_generator()
sample_data <- rbind(tree_sample,boundary)
p2 <- ggplot(sample_data, aes(x, y, color = type)) +
geom_point(alpha = 0.6) +
scale_color_manual(values = c("boundary" = "red", "sample" = "blue"))+
coord_fixed() +
ggtitle("Sample Points and Boundary Points") +
theme_minimal()
p2
Using the individual location data, we compute the pairwise distances between all individuals in the sampling plot, thus constructing a distance matrix. This matrix represents a complete graph for the set of points, where the edges include all potential connections.
#calculate distance matrix
dist_generator <- function(all_points){
#distance matrix
dist_matrix <- as.matrix(dist(all_points[, c("x", "y")]))
#lower part
lower_idx <- which(lower.tri(dist_matrix), arr.ind = TRUE)
#translate to a dataframe
distance_df <- data.frame(
i = lower_idx[, 1],
j = lower_idx[, 2],
distance = dist_matrix[lower_idx]
)
return(distance_df)
}distance_df <- dist_generator(sample_data)
ggplot(distance_df, aes(x = i, y = j, fill = distance)) +
geom_tile() +
scale_fill_viridis_c(option = "magma") +
coord_fixed() +
theme_minimal() +
ggtitle("Distance Matrix Heatmap")
Next, we create a set of neighbor relationships. To reduce computational complexity, we filter the edges within a specific distance range (for example, 3 or 4 times of the average distance) and sort them by distance.
# sort and preprocess distance dataframe
dist_sort <- function(distance_df, max_dist){
distance_df <- distance_df %>%
filter(distance < max_dist) %>%
arrange(distance)
return(distance_df)
}avg_dist <- sqrt(X * Y / rho)
distance_df_sort <- dist_sort(distance_df, 4*avg_dist)
ggplot(distance_df_sort, aes(x = i, y = j, fill = distance)) +
geom_tile() +
scale_fill_viridis_c(option = "magma") +
coord_fixed() +
theme_minimal() +
ggtitle("Distance Matrix Heatmap")
We then check each edge from smallest to largest to see if it intersects with any existing edges in the current neighbor relationship set. If it does not intersect, it is added to the set. This process continues until the complete set of neighbor relationships is generated.
# examine intersect with both end of the edge
is_intersect <- function(p1, p2, q1, q2) {
cross <- function(p, q, r) {
(q$x - p$x) * (r$y - p$y) - (q$y - p$y) * (r$x - p$x)
}
(cross(p1, p2, q1) * cross(p1, p2, q2) < 0) &&
(cross(q1, q2, p1) * cross(q1, q2, p2) < 0)
}neighbor_generator <- function(all_points, distance_df_sort){
selected_edges <- data.frame(i = integer(), j = integer(), distance = double())
for (k in 1:nrow(distance_df_sort)) {
i <- distance_df_sort$i[k]
j <- distance_df_sort$j[k]
d <- distance_df_sort$distance[k]
p1 <- all_points[i, c("x", "y")]
p2 <- all_points[j, c("x", "y")]
has_intersection <- FALSE
if (nrow(selected_edges) > 0) {
for (l in 1:nrow(selected_edges)) {
i_exist <- selected_edges$i[l]
j_exist <- selected_edges$j[l]
q1 <- all_points[i_exist, c("x", "y")]
q2 <- all_points[j_exist, c("x", "y")]
if (is_intersect(p1, p2, q1, q2)) {
has_intersection <- TRUE
break
}
}
}
if (d == 0 || !has_intersection) {
selected_edges <- rbind(selected_edges, data.frame(i = i, j = j, distance = d))
}
}
return(selected_edges)
}edges_data <- neighbor_generator(sample_data,distance_df_sort)
head(edges_data)## i j distance
## 1 72 45 0.009909599
## 2 30 19 0.012091432
## 3 86 18 0.012683174
## 4 63 49 0.014607308
## 5 69 20 0.019251597
## 6 85 5 0.020385177# both ends of the edges
edges_ends <- edges_data %>%
rowwise() %>%
mutate(x1 = sample_data$x[i], y1 = sample_data$y[i], x2 = sample_data$x[j], y2 = sample_data$y[j])
p3<-ggplot() +
geom_segment(data = edges_ends, aes(x = x1, y = y1, xend = x2, yend = y2), color = "darkgray", linewidth = 0.5, alpha = 0.7) +
geom_point(data = sample_data, aes(x = x, y = y, color = type), alpha = 0.6) +
scale_color_manual(values = c("boundary" = "red", "sample" = "blue"))+
ggtitle("Final Network of Selected Segments") +
theme_minimal() +
coord_fixed()
p3
Then, we delete boundary points and edges related to boundary, the final result contains all coexistence pairs.
# delete boundary points
edges_no_boundary <- edges_data %>%
filter(i %in% rownames(tree_sample) & j %in% rownames(tree_sample))
edges_end_no_boundary <- edges_no_boundary %>%
rowwise() %>%
mutate(x1 = tree_sample$x[as.numeric(i)], y1 = tree_sample$y[as.numeric(i)],
x2 = tree_sample$x[as.numeric(j)], y2 = tree_sample$y[as.numeric(j)])
p4<-ggplot() +
geom_segment(data = edges_end_no_boundary, aes(x = x1, y = y1, xend = x2, yend = y2), color = "darkgray", linewidth = 0.5, alpha = 0.7) +
geom_point(data = tree_sample, aes(x = x, y = y), color = "blue", size = 2, alpha = 0.7) +
ggtitle("Network") +
theme_minimal() +
coord_fixed()
p4
And we could inferring the interaction strength as well as interaction networks accordingly.
6 Discussion
The algorithm we implemented here follows a greedy approach. However, when the number of individuals increases significantly, this global algorithm becomes inefficient and time-consuming. To enhance its performance, several improvements can be considered:
Divide the area into smaller patches: This can help by reducing the problem size in each patch, making the computation more manageable.
Utilize spatial indexing techniques: When the area is extremely large, spatial index methods can speed up the process by quickly identifying relevant regions.
Limit the number of final edges: By setting a constraint on the maximum number of edges, we can avoid considering distant edges that are unnecessary. In a complete planar graph, it can be proven that the number of edges is bounded by \(3n−3\), where \(n\) is the number of points.
Additionally, another more efficient approach is Delaunay triangulation, which is significantly faster than the method we used here. Delaunay triangulation represents the dual graph of a Voronoi diagram and provides a more computationally efficient way to connect the points.
delaunay_result <- deldir(sample_data$x, sample_data$y)
delaunay_edges <- delaunay_result$delsgs
ggplot() +
geom_segment(data = delaunay_edges,
aes(x = x1, y = y1, xend = x2, yend = y2),
color = "darkgray", linewidth = 0.5) +
geom_point(data = sample_data, aes(x = x, y = y, color = type), alpha = 0.6) +
scale_color_manual(values = c("boundary" = "red", "sample" = "blue"))+
coord_fixed() +
theme_minimal() +
ggtitle("Delaunay Triangulation")
I became acquainted with the Voronoi graph only after completing this work; otherwise, I would have opted to use it from the outset. Delaunay triangulation, however, doesn’t seem to perfectly align with our model in some respects. While I have not yet fully clarified the exact differences between them, it is clear that they sometimes produce different edges.
Additionally, (to make this work appear at least somewhat meaningful), the algorithm we developed can be easily extended to accommodate various other conditions. For instance, if we assume that the area of an individual’s neighborhood is related to its biomass (which could be measured by height or DBH), we can adapt our method accordingly. In such cases, the only modification needed would be to adjust how we calculate the distance, such as by dividing it by the height.
7 Previous Code (by 2022)
install.packages("dplyr")
install.packages("ggplot2")
install.packages("data.table")# Construct neighbor web of trees
# By Chang Longxiao, 2019
# Read data, including real data and edges that are chosen
data.final <- read.table("data.final.csv", header = TRUE, sep = ",", encoding = "UTF-8")
data.edge <- read.table("edge.final.csv", header = TRUE, sep = ",", encoding = "UTF-8")
names(data.final) <- c("I", "sp", "x", "y", "h", "DBH")
names(data.edge) <- c("I", "sp", "x", "y", "h", "DBH")
n <- length(table(data.final$I))
# Divide data by site
data.list <- list()
for (i in 1:n) {
data.sub <- subset(data.final, I == i)
data.edge$h <- mean(data.sub$h)
data.edge$DBH <- mean(data.sub$DBH)
data.list[[i]] <- rbind(data.sub, data.edge)
}
# Generate the web and store by site
result.list <- list()
for (i in 1:n) {
a <- as.data.frame(data.list[i])
result.list[[i]] <- neibor.calculate(a, 5, weight = F)
}
# Function to generate neighbor network
neibor.calculate <- function(tree.plot, gapl, weight) {
num <- length(tree.plot[, 1])
num.count <- length(tree.plot[tree.plot$I != 0, 1])
m.dis <- matrix(0, num, num)
for (i in 1:num) {
for (j in 1:num) {
m.dis[i, j] <- sqrt((tree.plot$x[i] - tree.plot$x[j])^2 + (tree.plot$y[i] - tree.plot$y[j])^2)
if (weight)
m.dis[i, j] <- m.dis[i, j] / (tree.plot$DBH[i] + tree.plot$DBH[j])
}
}
# Structure distance matrix
exis <- matrix(0, num, num)
a <- 0
for (i in 1:(num - 1)) {
for (j in (i + 1):num) {
if (m.dis[i, j] <= gapl) a <- a + 1
else exis[i, j] <- -1
}
}
# Distances that match the gap
while (a > 0) {
p <- 1
q <- 2
l <- gapl
for (i in 1:(num - 1)) {
for (j in (i + 1):num) {
if (exis[i, j] == 0 & m.dis[i, j] < l) {
p <- i
q <- j
l <- m.dis[i, j]
}
}
}
print(l)
# Min distance now
if (l == 0) exis[p, q] <- 1
if (l > 0) {
for (i in 1:(num - 1)) {
if (i == p) next
for (j in (i + 1):num) {
if (j == q) next
if (exis[i, j] == 1 & m.dis[i, j] > 0) {
ij1 <- (tree.plot$x[j] - tree.plot$x[i]) * (tree.plot$y[p] - tree.plot$y[i]) -
(tree.plot$y[j] - tree.plot$y[i]) * (tree.plot$x[p] - tree.plot$x[i])
ij2 <- (tree.plot$x[j] - tree.plot$x[i]) * (tree.plot$y[q] - tree.plot$y[i]) -
(tree.plot$y[j] - tree.plot$y[i]) * (tree.plot$x[q] - tree.plot$x[i])
pq1 <- (tree.plot$x[q] - tree.plot$x[p]) * (tree.plot$y[i] - tree.plot$y[p]) -
(tree.plot$y[q] - tree.plot$y[p]) * (tree.plot$x[i] - tree.plot$x[p])
pq2 <- (tree.plot$x[q] - tree.plot$x[p]) * (tree.plot$y[j] - tree.plot$y[p]) -
(tree.plot$y[q] - tree.plot$y[p]) * (tree.plot$x[j] - tree.plot$x[p])
if (is.na(ij1 * ij2 < 0 & pq1 * pq2 < 0)) {
exis[p, q] <- -1
break
}
if (ij1 * ij2 < 0 & pq1 * pq2 < 0) {
exis[p, q] <- -1
break
}
}
}
if (exis[p, q] == -1) break
}
if (exis[p, q] == 0) exis[p, q] <- 1
}
a <- a - 1
print(a)
}
# Find edges of the graph
for (i in 2:num) {
for (j in 1:(i - 1)) {
if (exis[j, i] == -1) exis[j, i] <- 0
exis[i, j] <- exis[j, i]
}
}
# Fill the matrix
result <- list(m.dis = m.dis[1:num.count, 1:num.count], exis = exis[1:num.count, 1:num.count])
return(result)
}
# Translate the neighborhood matrix into co-occurrence matrix
inter.list <- list()
for (i in 1:n) {
inter.list[[i]] <- inter.frame(data.list[[i]], result.list[[i]]$exis, 45)
}
# Function to create the co-occurrence matrix
inter.frame <- function(sp.data, exis, sp.num) {
inter.num <- sum(exis) / 2
i.frame <- as.data.frame(matrix(0, inter.num, sp.num))
n <- length(sp.data[, 1]) - 88
flag <- 1
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
if (exis[i, j] == 1) {
i.frame[flag, sp.data$sp[i]] <- 1
i.frame[flag, sp.data$sp[j]] <- 1
flag <- flag + 1
}
}
}
return(i.frame)
}