R中为Polygon分配空间参考报错：找不到匹配方法求助

问题描述 #

尝试生成随机大小的地块以覆盖栅格图像的85%，但运行代码时持续报错：

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘proj4string<-’ for signature ‘"Polygon", "CRS"’

原代码：

library(sp)
library(ggplot2)
library(rgeos)
library(sf)
library(rgdal)
library(raster)

# image
site <- raster("C:/Users/ichth/OneDrive/Documents/Restoration_Implementation/Design_Project/Site_Raster.tif")

# plot size min/max in pixels
plot_size <- c(10, 25)

# find raster total area
total_area <- cellStats(site, sum)

# find area needed to cover 85%
target_area <- 0.85 * total_area

# create empty data frame 
plots <- SpatialPolygonsDataFrame(
  SpatialPolygons(list()), 
  data.frame(id = numeric(), area = numeric(), row.names = character()),
  match.ID = FALSE
)

# loop until target reached
while (sum(plots$area) < target_area) {
  
  # next plot random location and size
  x <- runif(1, xmin(site), xmax(site))
  y <- runif(1, ymin(site), ymax(site))
  size <- runif(1, plot_size[1], plot_size[2])
  
  # set coordinate system for plots
  crs(site) <- "+proj=tmerc +lat_0=24.3333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
  
  # create overall polygon
  pol <- Polygon(
    matrix(
      c(x, x + size, x + size, x, x, 
        y, y, y + size, y + size, y),
      ncol = 2, byrow = TRUE))
  
  # set coordinate system for plots
  proj4string(pol) <- CRS(as.character(crs(site)))
  
  # add to plots data frame
  spatial_poly <- SpatialPolygons(list(Polygons(list(pol), ID = as.character(length(plots) + 1))))
  plots <- rbind(plots, data.frame(id = length(plots) + 1, area = gArea(pol), row.names = as.character(length(plots) + 1)))
  coordinates(plots) <- c("x", "y")
  
  
  # check if target area has been reached
  if (sum(plots$area) >= target_area) {
    break
  }
  
  # plot on raster
  ggplot() +
    geom_raster(data = as.data.frame(site), aes(x = x, y = y, fill = value)) +
    geom_polygon(data = as.data.frame(plots), aes(x = x, y = y, fill = area), color = "white", alpha = 0.7) +
    coord_equal() +
    theme_void() +
    scale_fill_gradient(low = "white", high = "red")
}

错误原因 #

proj4string<-函数是为SpatialPolygons（空间多边形集合）对象设计的，不能直接给单个Polygon（仅几何图形）对象设置投影——Polygon类本身不具备存储空间参考的属性，因此触发方法不匹配的报错。

此外代码还有其他问题：

• 循环内重复设置栅格的CRS，完全没必要
• 合并SpatialPolygonsDataFrame的方式错误，直接rbind普通数据框会丢失空间属性
• 用gArea(pol)计算面积无效，gArea需要接收Spatial类对象
• 对已有的SpatialPolygonsDataFrame再次调用coordinates(plots) <- ...会破坏对象结构

修正后的代码 #

library(sp)
library(ggplot2)
library(rgeos)
library(raster)

# 加载栅格并设置CRS（只执行一次）
site <- raster("C:/Users/ichth/OneDrive/Documents/Restoration_Implementation/Design_Project/Site_Raster.tif")
crs(site) <- "+proj=tmerc +lat_0=24.3333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"

# 地块大小范围（单位与栅格坐标一致）
plot_size <- c(10, 25)

# 计算目标覆盖面积
total_area <- cellStats(site, sum)
target_area <- 0.85 * total_area

# 创建空的SpatialPolygonsDataFrame
plots <- SpatialPolygonsDataFrame(
  SpatialPolygons(list()), 
  data.frame(id = integer(), area = numeric(), row.names = character()),
  match.ID = FALSE
)

# 循环生成地块直到达到目标面积
while (sum(plots$area, na.rm = TRUE) < target_area) {
  # 随机生成位置和大小
  x <- runif(1, xmin(site), xmax(site))
  y <- runif(1, ymin(site), ymax(site))
  size <- runif(1, plot_size[1], plot_size[2])
  
  # 创建单个Polygon并包装为SpatialPolygons
  pol <- Polygon(matrix(
    c(x, x + size, x + size, x, x, 
      y, y, y + size, y + size, y),
    ncol = 2, byrow = TRUE
  ))
  pols <- Polygons(list(pol), ID = as.character(nrow(plots) + 1))
  spatial_poly <- SpatialPolygons(list(pols))
  
  # 给SpatialPolygons设置投影
  proj4string(spatial_poly) <- crs(site)
  
  # 计算当前地块面积
  poly_area <- gArea(spatial_poly)
  
  # 转换为SpatialPolygonsDataFrame并合并到plots
  new_plot <- SpatialPolygonsDataFrame(
    spatial_poly,
    data.frame(id = nrow(plots) + 1, area = poly_area, row.names = as.character(nrow(plots) + 1))
  )
  plots <- rbind(plots, new_plot)
  
  # 提前终止循环
  if (sum(plots$area) >= target_area) break
  
  # 可视化当前进度
  site_df <- as.data.frame(site, xy = TRUE)
  plots_df <- fortify(plots) %>% left_join(plots@data, by = "id")
  
  ggplot() +
    geom_raster(data = site_df, aes(x = x, y = y, fill = value)) +
    geom_polygon(data = plots_df, aes(x = long, y = lat, group = group, fill = area), color = "white", alpha = 0.7) +
    coord_equal() +
    theme_void() +
    scale_fill_gradient(low = "white", high = "red")
}

关键修复点 #

• 把CRS设置移到循环外，避免重复执行
• 先将Polygon包装为SpatialPolygons，再给这个对象设置投影
• 使用SpatialPolygonsDataFrame的正确合并方式，保证空间属性不丢失
• 用gArea(spatial_poly)计算有效面积
• 用fortify()转换空间对象适配ggplot的可视化要求

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