在 spatstat 中使用 sf 多边形对象作为窗口

Question

你好所有潜在的帮助者，

我有一个SpatialPolygonDataFrame从tigris包中获得的对象，我想在创建ppp对象时将其用作多边形窗口。这是我尝试过的：

# Needed packages
library(spatstat)
library(sf)

# Download geospatial data for Lee county in Alabama (home of the great Auburn University by the way!)
county <- tigris::county_subdivisions(state = "Alabama", county = "Lee")

# The polygon of Lee county is subdivided, so I convert it to a single polygon after converting it to an sf object
county_sf <- st_as_sf(county)
county_one <- st_union(county_sf)

# A quick plot of the object outputs what I am expecting
plot(county_one)

# Now I create a planar point pattern and I use county_one as the window
p <- ppp(x = -85.4, y = 32.5, window = as.owin((county_one)))

# But the plot here shows that the window is just a rectangle and not the polygon :(
plot(p)

感谢您的帮助。

Answer 1

注意：我已编辑此答案以包含完整的详细信息。

作为@TimSalabim提到这是下的方式sf，但在那之前，你必须将旧办理sp类，如SpatialPolygons。使用类似 as_Spatial的sf，然后负载maptools和使用as.owin或 as(x, "owin")在上Spatial对象。

此外，您只能使用平面（投影）空间中的spatstat坐标， 而不能使用地球曲面上的坐标。您必须投影到相关的平面坐标系。也许 <epsg.io/6345> 在这种情况下可用。要投影到该坐标系，请使用 sf::st_transform(county_one, crs = 6345). 之后，您转换为 Spatial，然后owin。注意：选择相关投影是一门科学，我对它了解不多，所以如果你想确保你不会得到太扭曲的结果，那么做一些研究。

特别是使用原始示例，您可以执行以下操作：

# Needed packages
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: nlme
#> Loading required package: rpart
#> 
#> spatstat 1.62-2       (nickname: 'Shape-shifting lizard') 
#> For an introduction to spatstat, type 'beginner'
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: TRUE
library(tigris)
#> To enable 
#> caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
#> 
#> Attaching package: 'tigris'
#> The following object is masked from 'package:graphics':
#> 
#>     plot

county <- county_subdivisions(state = "Alabama", county = "Lee", class = "sf", progress_bar = FALSE)
county_one <- st_union(county)
plot(county_one)

county_flat <- st_transform(county_one, crs = 6345)
plot(county_flat)

county_owin <- as.owin(as_Spatial(county_flat))

县内100个随机点：

p <- runifpoint(100, win = county_owin)
plot(p)