r语言 - 将sf转换为未标记的ppp

Question

我希望将sf对象转换为未标记的ppp。(根据这篇文章，现在支持从sf到ppp的转换。)

library(sf)
#Initialise sf object
pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS["Korea 2000 / Unified CS",n    BASEGEOGCRS["Korea 2000",n        DATUM["Geocentric datum of Korea",n            ELLIPSOID["GRS 1980",6378137,298.257222101,n                LENGTHUNIT["metre",1]]],n        PRIMEM["Greenwich",0,n            ANGLEUNIT["degree",0.0174532925199433]],n        ID["EPSG",4737]],n    CONVERSION["Korea Unified Belt",n        METHOD["Transverse Mercator",n            ID["EPSG",9807]],n        PARAMETER["Latitude of natural origin",38,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8801]],n        PARAMETER["Longitude of natural origin",127.5,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8802]],n        PARAMETER["Scale factor at natural origin",0.9996,n            SCALEUNIT["unity",1],n            ID["EPSG",8805]],n        PARAMETER["False easting",1000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8806]],n        PARAMETER["False northing",2000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8807]]],n    CS[Cartesian,2],n        AXIS["northing (X)",north,n            ORDER[1],n            LENGTHUNIT["metre",1]],n        AXIS["easting (Y)",east,n            ORDER[2],n            LENGTHUNIT["metre",1]],n    USAGE[n        SCOPE["unknown"],n        AREA["Korea, Republic of (South Korea)"],n        BBOX[28.6,122.71,40.27,134.28]],n    ID["EPSG",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
library(spatstat)
#Convert to ppp format: method 1
pp1 <- as.ppp(pp)
#Warning message:
#In as.ppp.sf(pp) : only first attribute column is used for marks
pp1 <- unmark(pp1)
#Convert to ppp format: method 2
pp2 <- as.ppp(st_coordinates(pp), st_convex_hull(st_union(pp)))
#Warning message:
#data contain duplicated points
identical(pp1, pp2)
[1] FALSE

第一种方法需要标记，我还没有能够找出如何关闭(指定marks=NULL不起作用)。我删除标记之后，但这是一个解决方案。

考虑到警告，第二个方法可能是错误的，尽管我基于这个解决方案。两个输出不相同

这些方法中有正确的吗?为什么方法2会产生重复的点?有没有更直接的转换方法，没有变通方法?

Answer 1

我不确定下面的答案适用于spatstat的每个版本，但我认为，如果你想生成一个未标记的ppp对象，你应该运行

# packages
suppressPackageStartupMessages({
library(sf)
library(spatstat)
})
# Initialise sf object
pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS["Korea 2000 / Unified CS",n    BASEGEOGCRS["Korea 2000",n        DATUM["Geocentric datum of Korea",n            ELLIPSOID["GRS 1980",6378137,298.257222101,n                LENGTHUNIT["metre",1]]],n        PRIMEM["Greenwich",0,n            ANGLEUNIT["degree",0.0174532925199433]],n        ID["EPSG",4737]],n    CONVERSION["Korea Unified Belt",n        METHOD["Transverse Mercator",n            ID["EPSG",9807]],n        PARAMETER["Latitude of natural origin",38,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8801]],n        PARAMETER["Longitude of natural origin",127.5,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8802]],n        PARAMETER["Scale factor at natural origin",0.9996,n            SCALEUNIT["unity",1],n            ID["EPSG",8805]],n        PARAMETER["False easting",1000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8806]],n        PARAMETER["False northing",2000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8807]]],n    CS[Cartesian,2],n        AXIS["northing (X)",north,n            ORDER[1],n            LENGTHUNIT["metre",1]],n        AXIS["easting (Y)",east,n            ORDER[2],n            LENGTHUNIT["metre",1]],n    USAGE[n        SCOPE["unknown"],n        AREA["Korea, Republic of (South Korea)"],n        BBOX[28.6,122.71,40.27,134.28]],n    ID["EPSG",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
# Convert to ppp format: method 1
(pp1 <- as.ppp(st_geometry(pp)))
#> Planar point pattern: 10 points
#> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units

^{由reprex包(v0.3.0)创建于2021-03-01}

您提供的两个方法是不相等的，因为在sf包中定义的as.ppp.sfc函数使用一个矩形owin对象:

# packages
suppressPackageStartupMessages({
library(sf)
library(spatstat)
})
#Initialise sf object
pp <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS["Korea 2000 / Unified CS",n    BASEGEOGCRS["Korea 2000",n        DATUM["Geocentric datum of Korea",n            ELLIPSOID["GRS 1980",6378137,298.257222101,n                LENGTHUNIT["metre",1]]],n        PRIMEM["Greenwich",0,n            ANGLEUNIT["degree",0.0174532925199433]],n        ID["EPSG",4737]],n    CONVERSION["Korea Unified Belt",n        METHOD["Transverse Mercator",n            ID["EPSG",9807]],n        PARAMETER["Latitude of natural origin",38,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8801]],n        PARAMETER["Longitude of natural origin",127.5,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8802]],n        PARAMETER["Scale factor at natural origin",0.9996,n            SCALEUNIT["unity",1],n            ID["EPSG",8805]],n        PARAMETER["False easting",1000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8806]],n        PARAMETER["False northing",2000000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8807]]],n    CS[Cartesian,2],n        AXIS["northing (X)",north,n            ORDER[1],n            LENGTHUNIT["metre",1]],n        AXIS["easting (Y)",east,n            ORDER[2],n            LENGTHUNIT["metre",1]],n    USAGE[n        SCOPE["unknown"],n        AREA["Korea, Republic of (South Korea)"],n        BBOX[28.6,122.71,40.27,134.28]],n    ID["EPSG",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
# Convert to ppp format: method 1
(pp1 <- as.ppp(st_geometry(pp)))
#> Planar point pattern: 10 points
#> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units
# Method 2
(pp2 <- as.ppp(
X = st_coordinates(pp), 
W = as.owin(st_bbox(pp))
))
#> Warning: data contain duplicated points
#> Planar point pattern: 10 points
#> window: rectangle = [950492.1, 960500] x [1939982.2, 1955543.8] units
identical(pp1, pp2)
#> [1] TRUE

而且，第一个方法没有返回任何警告消息，因为作者设置了check = FALSE。

anyDuplicated(pp1)
#> [1] TRUE

^{由reprex包(v0.3.0)创建于2021-03-01}

详情见此处

Answer 2

这两个命令应该给出相同的结果。得到不同的结果是spatstat包最近更改引起的临时问题。

(这些更改打破了sf,maptools和sp包所做的一些假设，导致第二个命令无法找到as.owin的正确方法，因此它默认为矩形窗口)

第一种方法没有问题。数据pp包含标记值，因此as.ppp处理标记，并警告它只使用第一列。如果您不想要这些标记，只需在后面使用unmark即可。

在第二种方法中得到的警告告诉您一些数据点位于相同的位置。在第一种方法中，您不会得到此警告，因为同一位置的点具有不同的标记值。

警告不是错误。除了as.owin的调度存在临时问题外，这两种方法都是可以的。

我建议你采用第一种方法。

如果您确实需要使用第二种方法，则需要更新所有这些包，并且在spatstat的情况下，安装github存储库中可用的开发版本。如果这仍然造成麻烦，请等待一周左右，直到所有四个软件包的完整软件包更新在CRAN上发布。