R语言raster包zonal函数计算结果出现重复值求助

Hi Diego, sorry to hear you're stuck with duplicate values in your zonal stats—let's break down what's likely going wrong and fix it step by step.

Key Observations #

Your results show that bands 31-36 have identical ET sums to bands 1-6, which shouldn't happen (your colleague's results show mostly zeros here). This suggests either:

• Your reclassified zone raster (edb) has incorrect zone assignments (e.g., pixels that should be in zone 31 are incorrectly mapped to zone 1, or vice versa)
• Your input rasters (ET and ED) aren't properly aligned, leading to mismatched pixel calculations
• The reclassification logic is missing edge cases or misdefining intervals

Step-by-Step Fixes #

1. Verify Raster Alignment #

First, make sure your ET and ED rasters have the same projection, resolution, and extent. Mismatched rasters can cause zonal stats to miscompute:

# Check if rasters are aligned
if (!compareRaster(ET, ED)) {
  # Resample ET to match ED's grid (use appropriate method for your data)
  ET <- resample(ET, ED, method = "bilinear")
}

2. Fix Reclassification Logic #

Your current code removes the last row of your reclassification matrix, which may leave some ED values unclassified or misclassified. Here's a more robust way to create your zones:

# Get ED's maximum value (with NA handling)
ed_max <- cellStats(ED, "max", na.rm = TRUE)

# Create intervals that cover all possible values
lowbound <- seq(0, ed_max, by = 1000)
upbound <- c(lowbound[-1], ed_max + 1000) # Ensure last interval covers edge values
band <- seq_along(lowbound) # Use a simple sequence for bands to avoid confusion

# Build reclassification matrix (no need to remove the last row!)
rclmat <- as.matrix(cbind(lowbound, upbound, band))

# Reclassify with explicit interval behavior (left-closed, right-open)
edb <- reclassify(ED, rclmat, right = FALSE)

3. Validate the Classified Zone Raster #

Check if edb has unexpected zone values before running zonal stats:

# List unique zones in edb
unique_zones <- unique(edb[])
print(unique_zones)

# Check how many pixels fall into zone 31 (should be 0 if your colleague's results are correct)
pixel_count_zone31 <- sum(edb[] == 31, na.rm = TRUE)
cat("Number of pixels in zone 31:", pixel_count_zone31, "\n")

4. Run Zonal Stats with Explicit NA Handling #

Make sure you're explicitly handling NA values in your zonal calculations:

# Calculate zonal sums with na.rm=TRUE to ignore NA pixels
sum_ET <- zonal(ET, edb, "sum", na.rm = TRUE)
sum_SH <- zonal(SH, edb, "sum", na.rm = TRUE)
sum_SC <- zonal(SC, edb, "sum", na.rm = TRUE)
sum_WH <- zonal(WH, edb, "sum", na.rm = TRUE)

5. Check for Duplicate Zones #

If you still see duplicates, verify that your edb raster doesn't have duplicate zone assignments for the same pixels:

# Check if there are any unexpected zone mappings
table(edb[], ED[], useNA = "always")

Why This Works #

• Raster Alignment: Ensures every pixel in ET maps exactly to the corresponding pixel in ED, eliminating cross-pixel calculation errors.
• Robust Reclassification: Covers all possible ED values (including edge cases) and uses explicit interval logic to avoid misassigning pixels to wrong zones.
• Validation Steps: Lets you confirm that your zone raster is correct before running compute-heavy zonal stats.

Give these steps a try, and your results should match your colleague's expected zeros instead of duplicate sums!

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