R语言按图层分组提取指定值范围的栅格单元坐标与值

Alright, let's figure out how to solve this problem. You've got 12 raster stacks each with 39 layers, and you need to pull out cells from specific layers (matched to your xydf's layerName) where the cell value falls within a fixed range around the corresponding T value. The extract() function alone isn't cutting it because it's designed for point extraction, not filtering entire layers by value range. Here's a step-by-step solution that works:

Step 1: Load Required Packages #

First, make sure you have these libraries installed and loaded—they'll handle raster operations, data grouping, and iterative processing:

library(raster)
library(dplyr)
library(purrr)

Step 2: Prepare Your Data (or Use Your Existing Data) #

I'll use simulated data to mirror your setup; replace this with your actual raster stacks and xydf dataset:

# Simulate a raster stack with 3 layers (match your 39-layer structure)
r1 <- raster(ncol=36, nrow=18, vals=1:(18*36))
r2 <- sqrt(r1)
r3 <- r1 * 0.1
r_stack <- stack(r1, r2, r3)
names(r_stack) <- c("layer.1", "layer.2", "layer.3") # Match layer names to xydf

# Simulate your xydf dataset
xydf <- data.frame(
  x = c(30, -40, 100),
  y = c(70, -60, -40),
  layerName = c("layer.1", "layer.1", "layer.2"),
  T = c(94.00, 555.00, 9.69)
)

# Define your fixed range (e.g., ±10; adjust this to your needs)
value_threshold <- 10

Step 3: Define a Custom Extraction Function #

This function will handle the heavy lifting: for each row in your grouped xydf, it targets the correct raster layer, filters cells within your T±threshold range, and pulls out the cell coordinates and values.

extract_matching_cells <- function(group_subset, raster_stack, threshold) {
  # Iterate over each row in the grouped subset
  pmap_dfr(group_subset, function(x, y, layerName, T, ...) {
    # Calculate the value range for this row's T value
    lower_bound <- T - threshold
    upper_bound <- T + threshold
    
    # Pull the target raster layer
    target_raster <- raster_stack[[layerName]]
    
    # Find all cells in the layer that fall within the value range
    matching_cells <- which(target_raster[] >= lower_bound & target_raster[] <= upper_bound)
    
    # If no cells match, skip this row
    if (length(matching_cells) == 0) return(NULL)
    
    # Extract coordinates and values for matching cells
    cell_coords <- xyFromCell(target_raster, matching_cells)
    cell_values <- extract(target_raster, matching_cells)
    
    # Format results into a tidy data frame
    data.frame(
      original_x = x,
      original_y = y,
      layerName = layerName,
      target_T = T,
      cell_x = cell_coords[, 1],
      cell_y = cell_coords[, 2],
      cell_value = cell_values
    )
  })
}

Step 4: Apply the Function Across Groups #

Use dplyr to group your xydf by layerName, then apply the custom function to each group:

# Group xydf by layerName, process each group, and combine results
final_results <- xydf %>%
  group_by(layerName) %>%
  group_split() %>%
  map_dfr(~ extract_matching_cells(.x, r_stack, value_threshold))

# View your extracted results
print(final_results)

Key Notes & Troubleshooting #

• Coordinate System Alignment: Double-check that your raster data and xydf share the same coordinate reference system (CRS). If not, use projectRaster() to reproject your rasters to match xydf's CRS.
• Large Rasters: If your rasters are extremely large, using target_raster[] might consume too much memory. For this case, consider using raster::calc() with a mask to filter cells, or process the raster in chunks.
• Variable Thresholds: If you need different thresholds for different rows, add a threshold column to xydf and pass it into the function instead of using a fixed value.

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