3

我是 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)

if (w.area < 0) stop(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()。所以我尝试加载一个马萨诸塞州的形状文件（此时我完全在猜测）。：

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, .. [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] ...[检查具有 2449 条边的多边形...] .. [etd 4:49] ....70 [ etd 4:27] ..... [etd 4:07] ....80 [etd 3:50] ..... [etd 3:36] ....90 [etd 3:22] 。 .... [等 3:11] ....100 [等。

我收到消息抱怨顶点相交等，但它未能构建多边形。

关于问题的一些背景：我正在尝试使用 spatstat 中的函数进行空间相对风险计算，即病例密度与对照的空间比率。为此，我需要一个观察窗口和该窗口内的点图。我可以作弊，使观察窗口围绕马萨诸塞州成为一个矩形，但这可能会扭曲海岸附近的值。在任何情况下，我都想学习如何在我以后使用这个包所做的任何工作中正确地做到这一点。感谢您的任何帮助，您可以提供。

4

3 回答 3

2

编辑 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://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)

在此处输入图像描述

于 2014-03-18T10:27:40.863 回答
1

最初的问题是在 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 章，您可以从此处免费下载。

于 2021-02-22T03:19:30.420 回答
0

我很熟悉你遇到的问题。如果我理解正确，您希望获得 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)
于 2014-01-17T15:40:53.960 回答