如何确保在spatstat::owin(poly=<polygon>)函数中输入的多边形不具有“负面积”？

Question

我是R spatstat包的新用户，使用owin()函数创建多边形观测窗口时遇到了问题。以下是代码：

library("maps")
library ("sp")` 
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T) # This returns a data frame includding x and y components that form a polygon of massachusetts mainland`

mass.win <- owin(poly=data.frame(x=mass.map$x, y=mass.map$y)

如果 (w.area < 0) 出错了，停止(paste("Area of polygon is negative -", "maybe traversed in >wrong direction?")) : 需要TRUE/FALSE的地方缺失数值

我尝试了一些方法，比如反转多边形的顺序，但是仍然出现了同样的错误。

 mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))

多边形包含重复的顶点

多边形自相交 在owin(poly = data.frame(x = rev(mass.map$x), y = rev(mass.map$y)))中出现错误：多边形数据包含重复的顶点和自相交

后来我想也许map()返回的多边形不适合传递给owin()。所以我尝试加载马萨诸塞州的一个shape文件（在这一点上我完全是瞎猜的）：

x <- readShapePoly("../Geog/OUTLINE25K_POLY") ## The shape file for MASS, loaded from MassGIS website
mass.poly <- x <- readShapePoly("../Geog/OUTLINE25K_POLY", force_ring=T, delete_null_obj=T) ## I got following error whether or not I used force_ring

mass.owin <- as(mass.poly, "owin") 检查1006个多边形...1，多边形1包含重复顶点 [正在检查具有91844个边的多边形...] 2、3、.. [预计剩余时间1小时21分52秒] ....10 [预计剩余时间36分12秒] ..... [预计剩余时间23分10秒] ....20 [预计剩余时间16分59秒] ..... [预计剩余时间13分22秒] ....30 [预计剩余时间11分01秒] ..... [预计剩余时间9分21秒] ....40 [预计剩余时间8分06秒] ..... [预计剩余时间7分09秒] ....50 [预计剩余时间6分23秒] ..... [预计剩余时间5分46秒] ....60 [预计剩余时间5分15秒] ...[正在检查具有2449个边的多边形...] .. [预计剩余时间4分49秒] ....70 [预计剩余时间4分27秒] ..... [预计剩余时间4分07秒] ....80 [预计剩余时间3分50秒] ..... [预计剩余时间3分36秒] ....90 [预计剩余时间3分22秒] ..... [预计剩余时间3分11秒] ....100 [ 等等。

我收到了关于相交顶点等问题的消息，并且无法构建多边形。

一些问题的背景：我正在尝试使用spatstat中的函数进行空间相对风险计算，即病例与对照组密度的空间比率。为此，我需要一个观察窗口和该窗口内的点图。我可以欺骗并使观察窗口成为围绕马萨诸塞州的矩形，但这可能会扭曲沿海地区的值。无论如何，我想学习如何正确地完成任何将来我使用此软件包的工作。感谢您提供的任何帮助。

- user2793761

---

首先，第一个错误信息告诉你它根本无法获取值。像其他消息一样，这几乎肯定是由于你的多边形具有自相交和/或重复部分。听起来你需要一个更干净（可能不那么详细）的多边形来表示你的地理区域。 - Carl Witthoft

Answer 1

EDIT 2021-02-04: 我再次遇到了 spatstat 中的顺时针/逆时针问题，并发现只有我自己的答案才能解决这个问题。答案并没有提供很多帮助。我已经改进了它，使其对更广泛的应用程序更有用。

spatstat 包需要 owin 对象坐标按逆时针方向指定。来自版本 1.36-0 的文档中的一句引用：

单个多边形：如果 poly 是一个具有两列的矩阵或数据框，或者是具有相等长度的两个分量向量 x 和 y 的结构体，则这些值被解释为围绕窗口的多边形的顶点的笛卡尔坐标。顶点必须按逆时针方向列出。不应重复任何顶点（即不要重复第一个顶点）。

library("maps")
library ("sp")
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T)

首先，您需要确定您的多边形是顺时针还是逆时针。可以使用这里提供的方法来找到答案。将其表述为一个函数：

#' @title Check whether points for an owin are clockwise
#' @param x a dataframe with x coordinates in the first column and y coordinates in the second. 
#' @details Similarly to owin, the polygon should not be closed
#' @return A logical telling whether the polygon is arranged clockwise.
#' @author The idea has been scavenged from https://dev59.com/r3M_5IYBdhLWcg3w6X5e#1165943

clockwise <- function(x) {
    
  x.coords <- c(x[[1]], x[[1]][1])
  y.coords <- c(x[[2]], x[[2]][1])
  
  double.area <- sum(sapply(2:length(x.coords), function(i) {
    (x.coords[i] - x.coords[i-1])*(y.coords[i] + y.coords[i-1])
  }))
  
  double.area > 0
} 

现在，请检查mass.map中的坐标是否顺时针方向：

clockwise(data.frame(x=mass.map$x, y=mass.map$y))
#> TRUE

他们是。使用rev函数将坐标逆时针旋转，该函数可以颠倒x和y的向量：

mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))
plot(mass.win)

- Mikko

Answer 2

原始问题发布于7年前（2014年）。spatstat包自那时以来已经发生了巨大变化。

关于重复顶点和自相交多边形的错误消息不再出现，因为现在的代码会自动纠正这些多边形数据中的缺陷。

错误消息 Area of polygon is negative - maybe traversed in wrong direction? 仍然会出现。这是因为owin要求多边形顶点围绕外边界按逆时针顺序列出（围绕孔的边界则按顺时针顺序）。如果您的窗口是单个多边形而不是多个多边形的集合，则解决方法是使用rev反转坐标向量的顺序：

library(maps)
library(sp)
library(spatstat)
mass.map <- map("state", "massachusetts:main", fill=TRUE)
mass.win <- owin(poly=lapply(mass.map, rev))

我们无法自动执行此操作，因为这并不总是用户想要的。欲了解更多信息，请参阅 spatstat书第3章，您可以从此处免费下载。

- Adrian Baddeley

Answer 3

我非常熟悉你所遇到的问题。 如果我理解正确，你想要获得一个owin类的多边形对象。以下代码将帮助你获得你想要的。

#Loading the required packages
library (sp)
library(maps)

# Reading the dataset
mass.map <- map("state", "massachusetts:main", fill=T)

IDs <- sapply(strsplit(mass.map$names, ":"), function(x) x[1])
# converting can be done by loading the powerful package "maptools"
library(maptools)
# Converting into SpatialPolygons-object
mass.poly <- map2SpatialPolygons(mass.map, IDs = IDs)

# Converting by coercing into owin-object
mass.owin <- as.owin.SpatialPolygons(mass.poly)

- And_R