将 sf 转换为未标记的 ppp

问题描述

我希望将 sf 对象转换为未标记ppp。 (根据 this post,现在支持sf 转换为 ppp。)

library(sf)

#Initialise sf object
pp <- structure(list(X = c(959207.877070254,959660.734838225,951483.685462513,951527.767554883,958310.673042469,950492.05212104,959207.877070254,960500.020456073,959660.734838225),Y = c(1944457.42827898,1955543.76027363,1939982.16629396,1940216.55143212,1954704.68186897,1951434.68524296,1944457.42827898,1955292.64874361,1955543.76027363),geometry = structure(list(structure(c(959207.877070254,1944457.42827898),class = c("XY","POINT","sfg")),structure(c(959660.734838225,structure(c(951483.685462513,1939982.16629396),structure(c(951527.767554883,1940216.55143212),structure(c(958310.673042469,1954704.68186897),structure(c(950492.05212104,1951434.68524296),structure(c(959207.877070254,structure(c(960500.020456073,1955292.64874361),"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\",\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,0.0174532925199433],8801]],\n        ParaMETER[\"Longitude of natural origin\",127.5,8802]],\n        ParaMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],8805]],\n        ParaMETER[\"False easting\",1000000,\n            LENGTHUNIT[\"metre\",8806]],\n        ParaMETER[\"False northing\",2000000,8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\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 不起作用)。之后我会删除标记,但这是一种解决方法

考虑到警告,第二种方法可能是错误的,尽管我基于 this solution。 2 个输出不相同。

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

解决方法

我不确定以下答案是否适用于 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,959207.877070254,960500.020456073,959660.734838225),Y = c(1944457.42827898,1955543.76027363,1939982.16629396,1940216.55143212,1954704.68186897,1951434.68524296,1944457.42827898,1955292.64874361,1955543.76027363),geometry = structure(list(structure(c(959207.877070254,1944457.42827898),class = c("XY","POINT","sfg")),structure(c(959660.734838225,structure(c(951483.685462513,1939982.16629396),structure(c(951527.767554883,1940216.55143212),structure(c(958310.673042469,1954704.68186897),structure(c(950492.05212104,1951434.68524296),structure(c(959207.877070254,structure(c(960500.020456073,1955292.64874361),"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\",\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,0.0174532925199433],8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],8805]],\n        PARAMETER[\"False easting\",1000000,\n            LENGTHUNIT[\"metre\",8806]],\n        PARAMETER[\"False northing\",2000000,8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\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 package (v0.3.0) 于 2021 年 3 月 1 日创建

您提供的两种方法并不等效,因为 as.ppp.sfc 包中定义的 sf 函数使用矩形 owin 对象:

# packages
suppressPackageStartupMessages({
  library(sf)
  library(spatstat)
})

#Initialise sf object
pp <- structure(list(X = c(959207.877070254,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,1955543.8] units

identical(pp1,pp2)
#> [1] TRUE

此外,由于作者设置了 check = FALSE,第一个方法没有返回警告消息。

anyDuplicated(pp1)
#> [1] TRUE

reprex package (v0.3.0) 于 2021 年 3 月 1 日创建

有关详细信息,请参阅 here

,

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

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

第一种方法没有任何问题。数据 pp 包含标记值,因此 as.ppp 处理标记,并警告它仅使用第一列。如果您不想要标记,请在之后使用 unmark

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

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

我建议您使用第一种方法。

如果您确实需要使用第二种方法,则需要更新所有这些软件包,对于 spatstat,请安装 the github repository 提供的开发版本。 如果仍然导致问题,请等待一周左右,直到所有四个包的完整包更新在 CRAN 上发布。