从坐标列表制作 Owin 对象

如何解决从坐标列表制作 Owin 对象

我正在尝试构建用于 R 程序 rase 的物种分布多边形。该程序需要一个 owin 对象，但示例数据还包括一个 SpatialpolygonDataFrame。您可以自己获取数据：data(rase_data,package = 'rase')

我从坐标列表开始（每个物种的纬度/经度）。感谢这个答案 here，我已经能够为列表的每个元素（每个物种）制作一个多边形。我需要找到一个 owin 对象。这是一些测试数据的 dput，然后是我用来获取当前位置的代码。

#dput(specieslist)
specieslist <- structure(list(Species = c("A","A","M","A"),lat = c(37.407002,35.65242,33.33891,37.65105,38.90657,39.06893,34.53998,38.18311,37.40006,36.252183,40.32555,39.575983,39.73548,37.82096,39.71557,38.7222,35.58556,36.3075,36.208,33.967875,35.375,39.73589,38.75774,36.61058,37.63605,36.586111,40.63344,39.80565,39.72601,39.70529,40.50957,37.81238),long = c(-122.232016,-120.77931,-116.91402,-119.88513,-121.05138,-120.86546,-119.85962,-120.37691,-122.23219,-121.775867,-121.91209,-121.554167,-121.70126,-120.14661,-121.61181,-120.98745,-120.9122,-121.4806,-121.816,-120.097752,-120.6456,-121.70175,-120.8443,-119.05645,-119.8728,-121.914722,-121.87499,-121.71465,-121.76862,-121.53125,-122.10229,-120.42828)),class = "data.frame",row.names = c(NA,-40L))

通过围绕点创建外壳来制作每个物种/点的多边形：

#create simple feature
library(sf)
df.sf <- specieslist %>%
  st_as_sf( coords = c("long","lat" ),crs = 4326 )
# perform fast visual check using mapview-package
#mapview::mapview( df.sf )

#group and summarise by species,and draw hulls
hulls <- df.sf %>%
  group_by( Species ) %>%
  summarise( geometry = st_combine( geometry ) ) %>%
  st_convex_hull()
##result
#mapview::mapview( list(df.sf,hulls ) )

现在我认为 df.sf（sf 点对象）成为 SpatialpolygonDataFrame 并且 hulls（sf 多边形对象）成为一个 owin 对象：

as(df.sf,"Spatial") -> df.sf_SPDF #this formats incorrectly though.

distribution <- st_transform(hulls,crs = 6345)
dist_owin <- as.owin(as_Spatial(distribution))

#Error: Only projected coordinates may be converted to spatstat class objects

#OR

as.owin(distribution)
#Error: 'owin' is not an exported object from 'namespace:spatstat'
maptools::as.owin(distribution)
#Error: 'as.owin' is not an exported object from 'namespace:maptools'

问题是：df.sf_SPDF 似乎格式不正确和 dist_owin 错误。我发现 R 中的所有这些空间工作都非常令人困惑。我已经为此工作了好几天了。

更新：如果我尝试另一种方法 - 将几何图形转换为多边形，然后制作 owin。这会产生一个错误：

hulls_poly <- st_cast(distribution$geometry,"polyGON") #. 
dist_owin <- as.owin(as_Spatial(hulls_poly))
#ERROR:  no method or default for coercing “sfc_polyGON” to “owin”

解决方法

我对 sf 的了解不足以解决这个问题，所以我通过 terra 展示它，但重要的部分是操作顺序。如果您愿意，您可以在 sf 中再次实现它。应该不需要恢复到旧的 Spatial* 对象

您的数据

specieslist <- structure(list(Species = c("A","A","M","A"),lat = c(37.407002,35.65242,33.33891,37.65105,38.90657,39.06893,34.53998,38.18311,37.40006,36.252183,40.32555,39.575983,39.73548,37.82096,39.71557,38.7222,35.58556,36.3075,36.208,33.967875,35.375,39.73589,38.75774,36.61058,37.63605,36.586111,40.63344,39.80565,39.72601,39.70529,40.50957,37.81238),long = c(-122.232016,-120.77931,-116.91402,-119.88513,-121.05138,-120.86546,-119.85962,-120.37691,-122.23219,-121.775867,-121.91209,-121.554167,-121.70126,-120.14661,-121.61181,-120.98745,-120.9122,-121.4806,-121.816,-120.097752,-120.6456,-121.70175,-120.8443,-119.05645,-119.8728,-121.914722,-121.87499,-121.71465,-121.76862,-121.53125,-122.10229,-120.42828)),class = "data.frame",row.names = c(NA,-40L))

首先，我制作了一个空间对象，在本例中为 SpatVector，然后我将其转换为平面 CRS --- 以消除这种影响。

您选择的 epsg:6345，即 +proj=utm +zone=16 不适合您的数据。 16 区是阿拉巴马州的经度。加利福尼亚州涵盖两个 UTM 区域，因此您不能使用它。而是使用例如“Teale Albers”，如果您的所有数据都仅限于金州。

 library(terra)
 #terra version 1.2.5
 v <- vect(specieslist,c("long","lat"),crs="epsg:4326")
 tacrs <- "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m"
 v <- project(v,tacrs)

为了简单起见，我展示了 1 个物种的工作流程

usp <- unique(v$Species)  
sp <- v[v$Species==usp[1]]

制作一个凸包，我想你会想要添加一个缓冲区。

ch <- terra::convexhull(sp)
bch <- buffer(ch,25000)

plot(bch)
points(sp)

现在通过 sf

创建 owin

library(sf)
library(spatstat)
sfobj <- st_as_sf(bch)
owin <- as.owin(sfobj)

你可以像这样在新的CRS中提取点

pxy <- terra::coords(sp)

现在创建一个 spatstat ppp 对象。

x <- ppp(pxy[,1],pxy[,2],window=owin)
#Warning message:
#data contain duplicated points

为避免上述警告，您可以在脚本开头使用 specieslist <- unique(specieslist)

x
#Planar point pattern: 27 points
#window: polygonal boundary
#enclosing rectangle: [-222286.97,312378.62] x [-539742.6,217425] units