r - 在 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而不是地球曲面上的坐标。你有
投影到相关的平面坐标系。也许 是
在这种情况下可以使用。要投影到此坐标系，请使用 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)