问题描述
我希望将 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
包的更改引起的临时问题。
(这些更改破坏了包 sf
、maptools
和 sp
所做的一些假设,导致第二个命令无法找到 as.owin
的正确方法,因此它默认为矩形窗口)。
第一种方法没有任何问题。数据 pp
包含标记值,因此 as.ppp
处理标记,并警告它仅使用第一列。如果您不想要标记,请在之后使用 unmark
。
您在第二种方法中得到的警告告诉您某些数据点位于同一位置。在第一种方法中您不会收到此警告,因为同一位置的点具有不同的标记值。
警告不是错误。除了 as.owin
的调度存在临时问题外,这两种方法都可以。
我建议您使用第一种方法。
如果您确实需要使用第二种方法,则需要更新所有这些软件包,对于 spatstat
,请安装 the github repository 提供的开发版本。
如果仍然导致问题,请等待一周左右,直到所有四个包的完整包更新在 CRAN 上发布。