如何在R的sf包中按年份分组计算空间点近邻距离

Hey there! Let's work through this problem using just the sf package plus some tidyverse tools (like dplyr and tidyr) since you want a straightforward, sf-native solution. I'll break down each step so it's easy to follow as an R beginner.

Step 1: Load Required Packages #

First, make sure you have these packages installed (run install.packages(c("sf", "dplyr", "tidyr", "readr", "stringr")) if you don't):

library(sf)
library(dplyr)
library(tidyr)
library(readr)
library(stringr)

Step 2: Prepare Your Data #

You already have your test data set up, but let's include it here for completeness:

trial <- read.table(text = "Country station_code lat_laea lon_laea year
Belize BF09-1 -2955950 1247610 2009
Belize BF09-10 -2953600 1248590 2009
Belize BF09-11 -2954620 1247900 2009
Belize BF11-13 -2958360 1244020 2011
Belize BF11-18 -2963740 1240290 2011
Belize BF11-19 -2963380 1242020 2011
Costa BraulioCarrilloNP-C16 -3640760 1821170 2011
Costa BraulioCarrilloNP-C17 -3640730 1823240 2011
Costa BraulioCarrilloNP-C18 -3642140 1817560 2011
Guatemala 40063 -3178260 1249780 2009
Guatemala 40596 -3183800 1246940 2009
Guatemala 43279 -3182640 1251560 2009", header = TRUE)

# Convert to sf object with the correct CRS (you already did this part!)
trial.sp <- st_as_sf(trial, coords = c("lat_laea", "lon_laea"), crs = 3035)

Step 3: Calculate & Filter Close Yearly Pairs #

This is the core part. We'll group the data by year, compute distances within each year, then filter for pairs under 1000 meters:

close_pairs <- trial.sp %>%
  # Group data by year so we only compare points within the same year
  group_by(year) %>%
  # Apply a custom function to each year's group
  group_modify(function(data, key) {
    # Calculate distance matrix for all points in the current year
    dist_matrix <- st_distance(data, data)
    
    # Convert the matrix to a long-format data frame for easy filtering
    dist_df <- as.data.frame(dist_matrix) %>%
      # Add an index for the first point in each pair
      mutate(point1 = 1:nrow(data)) %>%
      # Reshape from wide to long format
      pivot_longer(cols = -point1, names_to = "point2", values_to = "distance") %>%
      # Convert point2's column name (V1, V2, etc.) to a numeric index
      mutate(point2 = as.integer(str_remove(point2, "V"))) %>%
      # Keep only pairs where distance < 1000m, and avoid duplicate pairs (e.g., A-B and B-A)
      filter(distance < units::set_units(1000, "meters"), point1 < point2)
    
    # If there are valid pairs, merge in the station/country info for readability
    if(nrow(dist_df) > 0) {
      dist_df %>%
        left_join(data %>% select(station_code, Country) %>% mutate(point1 = 1:nrow(data)), by = "point1") %>%
        left_join(data %>% select(station_code, Country) %>% mutate(point2 = 1:nrow(data)), by = "point2") %>%
        rename(station_code1 = station_code.x, Country1 = Country.x,
               station_code2 = station_code.y, Country2 = Country.y) %>%
        # Add the year back in
        mutate(year = key$year) %>%
        # Reorder columns for clarity
        select(year, Country1, station_code1, Country2, station_code2, distance)
    } else {
      # If no valid pairs in this year, return an empty frame with the right columns
      data.frame(year = key$year, Country1 = character(), station_code1 = character(),
                 Country2 = character(), station_code2 = character(), distance = numeric())
    }
  }) %>%
  # Remove grouping to get a clean data frame
  ungroup()

# Check the result
print(close_pairs)

Step 4: Export to CSV #

Finally, export your filtered pairs to a CSV file:

write_csv(close_pairs, "close_pairs_under_1000m.csv")

Key Notes #

• We use group_by(year) to ensure we only compare points from the same year—no cross-year pairs included.
• The point1 < point2 filter prevents duplicate entries (we don't want both "BF09-1 & BF09-10" and "BF09-10 & BF09-1" in the output).
• All distance calculations use your existing projected CRS (3035), so distances are in meters as you need.

内容的提问来源于stack exchange，提问作者AnnK