【问题标题】:Use sf polygon object as window in spatstat在 spatstat 中使用 sf 多边形对象作为窗口
【发布时间】:2020-10-01 23:22:52
【问题描述】:

您好所有潜在的帮助者,

我有一个从tigris 包中获得的SpatialPolygonDataFrame 对象,我想在创建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)

感谢您的帮助。

【问题讨论】:

标签: r geospatial sf spatstat


【解决方案1】:

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

正如@TimSalabim 提到的,这正在sf 中进行,但在那之前你必须 遍历旧的sp 类,例如SpatialPolygons。使用类似的东西 as_Spatial in sf 然后加载 maptools 并使用 as.owinas(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)

【讨论】:

  • 感谢您的回答。 sf 函数是 as_Spatial。不幸的是,当我使用 as.owin: Error in maptools:::as.owin.SpatialPolygons(county_sp) : Only projected coordinates may be converted to spatstat class objects 时出现此错误。
  • 您只能使用spatstat 的平面(投影)空间中的坐标,而不能使用地球曲面上的坐标。您必须投影到相关的坐标系。也许epsg.io/6345 是相关的。要投影到此坐标系,请使用sf::st_transform(county_one, crs = 6345)。之后您转换为Spatial,然后转换为owin。请注意,选择相关投影是一门科学,我对此了解不多,因此如果您想确保不会得到太扭曲的结果,请进行一些研究。
【解决方案2】:

这里只想指出,sf 类的强制方法现在由 sf 包注册(如果这是正确的话)。我不完全理解找到方法的 R 魔法,但它确实有效。

library(sf)
library(spatstat)

> methods(as.owin)
 [1] as.owin.boxx                    as.owin.data.frame              as.owin.default                 as.owin.distfun                
 [5] as.owin.dppm                    as.owin.funxy                   as.owin.im                      as.owin.influence.ppm          
 [9] as.owin.kppm                    as.owin.layered                 as.owin.leverage.ppm            as.owin.linfun                 
[13] as.owin.linnet                  as.owin.lintess                 as.owin.lpp                     as.owin.lppm                   
[17] as.owin.msr                     as.owin.MULTIPOLYGON*           as.owin.nnfun                   as.owin.owin                   
[21] as.owin.POLYGON*                as.owin.ppm                     as.owin.ppp                     as.owin.psp                    
[25] as.owin.quad                    as.owin.quadratcount            as.owin.quadrattest             as.owin.rmhmodel               
[29] as.owin.sf*                     as.owin.sfc*                    as.owin.sfc_MULTIPOLYGON*       as.owin.sfc_POLYGON*           
[33] as.owin.SpatialGridDataFrame*   as.owin.SpatialPixelsDataFrame* as.owin.SpatialPolygons*        as.owin.tess                   
see '?methods' for accessing help and source code

因此,假设您已正确预测数据(如 @Ege Rubak 所述),直接调用 as.owin 应该可以:

library(tigris)

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

county <- st_union(county)

county <- st_transform(county, crs = 6345)

window <- as.owin(county)

plot(window)

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2022-01-17
    • 2020-10-03
    • 2021-01-11
    • 2020-11-14
    • 1970-01-01
    • 2022-07-26
    • 1970-01-01
    • 2019-11-22
    相关资源
    最近更新 更多