如何在R中利用RasterBrick土地覆盖图实现面积加权统计

Hi Andy, let's work through your questions step by step—you're already off to a great start with importing and cropping your land cover data! Below are targeted solutions for each of your core needs:

1. Convert Cropped RasterBrick to High-Resolution Grid Shapefile #

If you want to turn your raster into a grid of polygon cells (one per raster pixel), use raster::rasterToPolygons (or the faster terra::as.polygons for large datasets). Note: High-resolution grids can get very large, so only do this if you explicitly need vector polygons (for most analysis, sticking to raster format is more efficient).

# First, extract the single land cover band from your RasterBrick
germany_land_single <- germany_land[[1]]

# Convert raster to polygon grid (dissolve=FALSE keeps every pixel as a separate polygon)
land_grid <- rasterToPolygons(germany_land_single, dissolve = FALSE)

# Save as Shapefile
writeOGR(
  obj = land_grid,
  dsn = "~/germany_land_grid",  # Output folder
  layer = "land_cover_grid",    # Shapefile name prefix
  driver = "ESRI Shapefile",
  overwrite_layer = TRUE
)

Pro Tip: #

For better performance with large rasters, switch to the terra package (the modern replacement for raster):

library(terra)
germany_land_terra <- rast(germany_land_single)
land_grid_terra <- as.polygons(germany_land_terra, dissolve = FALSE)
writeVector(land_grid_terra, "~/germany_land_grid/land_cover_grid.shp")

2. Efficiently Extract Target Land Cover Data for Germany NUTS-3 Regions #

The most efficient way to link raster land cover to vector NUTS-3 regions is using the exactextractr package—it handles partial raster-poly overlaps accurately (unlike the base raster::extract which only uses pixel centers).

Step 1: Prep Your Data #

First, ensure your NUTS-3 shapefile matches the CRS of your land cover raster:

library(exactextractr)
library(dplyr)

# Align CRS between NUTS-3 and land cover
Germany_NUTS3 <- spTransform(Germany_NUTS3, CRSobj = crs(germany_land))

# Define your target land cover classes (use Corine CLC codes)
# Example: Urban (111-122) + Agricultural (211-243)
target_classes <- c(111, 112, 121, 122, 211, 212, 221, 222, 231, 241, 242, 243)

Step 2: Extract & Summarize Land Cover by NUTS-3 #

# Extract area of each target class per NUTS-3 region
land_cover_summary <- exact_extract(
  x = germany_land[[1]],
  y = Germany_NUTS3,
  fun = function(value, coverage_fraction) {
    tibble(
      land_class = value,
      area_m2 = coverage_fraction * prod(res(germany_land[[1]]))  # Calculate pixel area
    )
  }
) %>%
  bind_rows(.id = "nuts3_row_id") %>%  # Link to NUTS-3 row index
  filter(land_class %in% target_classes) %>%
  group_by(nuts3_row_id, land_class) %>%
  summarise(total_area_m2 = sum(area_m2)) %>%
  ungroup()

# Merge results back to NUTS-3 attribute table
Germany_NUTS3@data <- Germany_NUTS3@data %>%
  mutate(nuts3_row_id = row_number()) %>%
  left_join(land_cover_summary, by = "nuts3_row_id")

Alternative (Base Raster Package): #

If you can't install exactextractr, use raster::extract—but note it's less accurate for partial overlaps:

# Extract all land cover pixels per NUTS-3 region
extracted_classes <- extract(germany_land[[1]], Germany_NUTS3)

# Calculate total urban/agricultural area per region
Germany_NUTS3@data$total_urban_m2 <- sapply(extracted_classes, function(x) {
  sum(x %in% c(111,112,121,122)) * prod(res(germany_land[[1]]))
})
Germany_NUTS3@data$total_agri_m2 <- sapply(extracted_classes, function(x) {
  sum(x %in% c(211,212,221,222,231,241,242,243)) * prod(res(germany_land[[1]]))
})

3. Calculate Area-Weighted Annual Mean Temperature #

To compute this, you'll need a temperature raster (e.g., WorldClim, ERA5) aligned with your land cover and NUTS-3 data. Here's how to combine all layers:

Step 1: Prep Temperature Raster #

# Load and align your temperature raster (example path)
temp_raster <- raster("~/germany_annual_temp.tif")

# Resample to match land cover raster resolution/CRS
temp_raster <- resample(temp_raster, germany_land[[1]], method = "bilinear")

Step 2: Compute Area-Weighted Temp by Land Cover Class #

# Extract temperature and land cover data per NUTS-3 region
weighted_temp_data <- exact_extract(
  x = temp_raster,
  y = Germany_NUTS3,
  fun = function(value, coverage_fraction, land_cover) {
    tibble(
      temp = value,
      area_m2 = coverage_fraction * prod(res(temp_raster)),
      land_class = land_cover
    )
  },
  append_cols = "land_cover",  # Attach land cover values to extraction
  data = germany_land[[1]]
) %>%
  bind_rows(.id = "nuts3_row_id") %>%
  filter(land_class %in% target_classes) %>%

# Calculate weighted mean temp per class, then overall area-weighted temp
  group_by(nuts3_row_id, land_class) %>%
  summarise(
    class_mean_temp = weighted.mean(temp, w = area_m2),
    class_total_area = sum(area_m2)
  ) %>%
  ungroup() %>%
  group_by(nuts3_row_id) %>%
  summarise(
    area_weighted_temp = sum(class_mean_temp * (class_total_area / sum(class_total_area)))
  ) %>%
  ungroup()

# Merge results to NUTS-3 shapefile
Germany_NUTS3@data <- Germany_NUTS3@data %>%
  mutate(nuts3_row_id = row_number()) %>%
  left_join(weighted_temp_data, by = "nuts3_row_id")

Simplified Option (No Land Cover Breakdown): #

If you just need the area-weighted temp for the entire NUTS-3 region (ignoring land cover):

Germany_NUTS3@data$nuts3_weighted_temp <- exact_extract(
  temp_raster, Germany_NUTS3, fun = "weighted.mean"
)

Key Notes: #

• CRS Consistency: Always verify all rasters/vectors share the same coordinate reference system before analysis.
• Memory Management: For large datasets, use terra instead of raster—it's faster and uses memory more efficiently.
• Corine Codes: Double-check your target land cover codes against the official Corine Land Cover classification to avoid errors.

内容的提问来源于stack exchange，提问作者Andy.Jian