r - 如何确保 spatstat::owin(poly=) 中的多边形没有 "negative area"

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)

Error in if (w.area < 0) stop(paste("Area of polygon is negative -", "maybe traversed in >wrong direction?")) : missing value where TRUE/FALSE needed

我尝试了诸如反转多边形顺序之类的方法，但得到了同样的错误。

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

多边形包含重复的顶点

Polygon is self-intersecting Error in owin(poly = data.frame(x = rev(mass.map$x), y = rev(mass.map$y))) : Polygon data contain duplicated vertices and self-intersection

然后我想也许map()返回的多边形并不意味着要输入owin()。所以我尝试加载一个马萨诸塞形状文件(我现在完全在猜测)。:

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") Checking 1006 polygons...1, Polygon 1 contains duplicated vertices [Checking polygon with 91844 edges...] 2, 3, .. [etd 1:21:52] ....10 [etd 36:12] ..... [etd 23:10] ....20 [etd 16:59] ..... [etd 13:22] ....30 [etd 11:01] ..... [etd 9:21] ....40 [etd 8:06] ..... [etd 7:09] ....50 [etd 6:23] ..... [etd 5:46] ....60 [etd 5:15] ...[Checking polygon with 2449 edges...] .. [etd 4:49] ....70 [etd 4:27] ..... [etd 4:07] ....80 [etd 3:50] ..... [etd 3:36] ....90 [etd 3:22] ..... [etd 3:11] ....100 [ etc.

我收到提示顶点相交等的消息，但未能构建多边形。

有关问题的一些背景:我正在尝试使用 spatstat 中的函数进行空间相对风险计算，即病例密度与对照密度的空间比率。为此，我需要一个观察窗口和该窗口内的点图。我可以作弊，将观察窗口设置为马萨诸塞州周围的矩形，但这可能会扭曲海岸附近的值。无论如何，我想了解如何在我将来使用此包进行的任何工作中正确执行此操作。感谢您的任何帮助，您可以提供。

Answer 1

编辑 2021-02-04:我再次遇到了这个顺时针/逆时针问题，并发现我自己的答案是解决 spatstat 问题的唯一答案。答案不是很有帮助。对其进行了改进，使其可用于更广泛的应用程序。

spatstat 包希望 owin 对象坐标指定逆时针。引用1.36-0版本的文档:

single polygon: If poly is a matrix or data frame with two columns, ora structure with two component vectors x and y of equal length, thenthese values are interpreted as the cartesian coordinates of thevertices of a polygon circumscribing the window. The vertices must belisted anticlockwise. No vertex should be repeated (i.e. do not repeatthe first vertex).

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

首先，您需要确定多边形是顺时针还是逆时针运行。做法here可以用来找出答案。将其表示为函数:

#' @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://stackoverflow.com/a/1165943/1082004

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)