在R中基于像素计数法实现土地覆盖栅格的多边形面积估算

Hey there! I’ve been right where you are—staring at a large Corine raster, frustrated with slow extraction times, and wondering how to turn pixel values into usable area stats for dozens of polygons. Let’s break down a faster, more straightforward approach using modern R tools.

Why raster::extract is Slow (and What to Use Instead) #

The old raster package works, but it’s not optimized for large datasets or bulk polygon extractions. The terra package—built as a replacement for raster—uses faster underlying code and has cleaner syntax for exactly this kind of task. It’ll cut down your processing time significantly.

Step-by-Step Implementation #

1. Load Your Data with terra #

First, load your clipped Corine raster and projected polygons:

# Install terra if you haven't already
# install.packages("terra")
library(terra)

# Load clipped land cover raster
lc_raster <- rast("path/to/your/clipped_corine.tif")

# Load your projected polygons (supports shapefiles, GeoPackage, etc.)
polygons <- vect("path/to/your/projected_polygons.shp")

2. Extract Pixel Counts Efficiently #

Use terra::extract with the fun=table parameter to directly count land cover classes per polygon. This avoids loading all pixel values into memory at once:

# Extract land cover class counts for each polygon (ignore NA values)
lc_counts <- extract(lc_raster, polygons, fun=table, na.rm=TRUE)

• fun=table tells the function to return a frequency table of land cover classes for each polygon.
• na.rm=TRUE skips any no-data pixels in your raster.

3. Convert Counts to Area #

Next, calculate the area for each land cover class. First get the area of a single pixel (make sure your raster uses a meter-based projection, like UTM, for accurate area calculations):

# Calculate area of one pixel (square meters)
pixel_area <- res(lc_raster)[1] * res(lc_raster)[2]

# Convert the count list to a tidy data frame with area
lc_area_stats <- do.call(rbind, lapply(lc_counts, function(poly_counts) {
  # Handle polygons with no overlapping pixels (if any)
  if (length(poly_counts) == 0) {
    data.frame(
      land_cover_code = NA,
      pixel_count = 0,
      area_sqm = 0
    )
  } else {
    data.frame(
      land_cover_code = names(poly_counts),
      pixel_count = as.integer(poly_counts),
      area_sqm = as.integer(poly_counts) * pixel_area
    )
  }
}))

# Add polygon IDs to link stats back to your original polygons
# Replace "poly_id" with the actual ID column name in your vector data
lc_area_stats$poly_id <- rep(polygons$poly_id, sapply(lc_counts, length))

4. (Optional) Add Land Cover Class Names #

If your Corine raster uses numeric codes, you can map them to human-readable names using Corine’s official classification:

# Example Corine classification lookup table (expand with all classes you need)
corine_lookup <- data.frame(
  land_cover_code = c("111", "112", "211", "311"),
  land_cover_name = c(
    "Continuous Urban Fabric",
    "Discontinuous Urban Fabric",
    "Non-Irrigated Arable Land",
    "Broad-Leaved Forest"
  )
)

# Merge to add class names
lc_area_stats <- merge(lc_area_stats, corine_lookup, by="land_cover_code", all.x=TRUE)

If You Must Stick with raster #

If you can’t switch to terra right now, here’s an optimized version using raster and dplyr to avoid excessive memory usage:

library(raster)
library(dplyr)

lc_raster <- raster("path/to/your/clipped_corine.tif")
polygons <- shapefile("path/to/your/projected_polygons.shp")

# Extract pixel values with polygon IDs
extracted_pixels <- extract(lc_raster, polygons, df=TRUE, na.rm=TRUE)

# Calculate area stats
pixel_area <- res(lc_raster)[1] * res(lc_raster)[2]
lc_area_stats <- extracted_pixels %>%
  group_by(ID, layer) %>%
  summarise(
    pixel_count = n(),
    area_sqm = n() * pixel_area
  ) %>%
  rename(poly_id = ID, land_cover_code = layer)

Note: This will still be slower than terra for large datasets, but it’s better than the default extract workflow.

Key Tips for Speed #

• Ensure projections match: You mentioned you already did this, but double-check—mismatched projections will slow things down and cause errors.
• Use compressed raster formats: If your Corine raster isn’t already in a compressed format (like COG or GeoTIFF with LZW compression), converting it will reduce load times.
• Avoid loading full rasters into memory: terra handles on-disk rasters efficiently, so you don’t need to read the entire 92MB file into RAM at once.

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