从坐标列表中创建 Owin 对象
make Owin object from list of coordinates
我正在尝试构建用于 R 程序的物种分布多边形 rase。该程序需要一个 owin 对象，但示例数据还包括一个 SpatialPolygonDataFrame。您可以自己获取数据：data(rase_data, package = 'rase')
我从坐标列表开始（lat/long 每个物种）。感谢这个答案 ，我已经能够为列表的每个元素（每个物种）制作一个多边形。我需要找到一个 owin 对象。这是一些测试数据的输入，然后是我用来获取当前位置的代码。
#dput(specieslist)
specieslist <- structure(list(Species = c("A", "A", "A", "A", "A", "M", "A", "M", "A", "A", "A", "A", "A", "A", "M", "M", "A", "M", "A", "A", "A", "M", "M", "M", "A", "A", "A", "A", "A", "A", "A", "M", "A", "A", "M", "M", "A", "M", "M", "A"), lat = c(37.407002, 35.65242, 33.33891, 37.65105, 38.90657, 39.06893, 34.53998, 38.18311, 37.40006, 35.65242, 34.53998, 33.33891, 37.65105, 38.90657, 38.18311, 39.06893, 36.252183, 40.32555, 39.575983, 39.73548, 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, -120.77931, -119.85962, -116.91402, -119.88513, -121.05138, -120.37691, -120.86546, -121.775867, -121.91209, -121.554167, -121.70126, -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))
通过围绕点创建外壳来根据 species/points 制作多边形：
#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
#或
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", "A", "A", "A", "M", "A", "M", "A", "A", "A", "A", "A", "A", "M", "M", "A", "M", "A", "A", "A", "M", "M", "M", "A", "A", "A", "A", "A", "A", "A", "M", "A", "A", "M", "M", "A", "M", "M", "A"), lat = c(37.407002, 35.65242, 33.33891, 37.65105, 38.90657, 39.06893, 34.53998, 38.18311, 37.40006, 35.65242, 34.53998, 33.33891, 37.65105, 38.90657, 38.18311, 39.06893, 36.252183, 40.32555, 39.575983, 39.73548, 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, -120.77931, -119.85962, -116.91402, -119.88513, -121.05138, -120.37691, -120.86546, -121.775867, -121.91209, -121.554167, -121.70126, -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
我正在尝试构建用于 R 程序的物种分布多边形 rase。该程序需要一个 owin 对象，但示例数据还包括一个 SpatialPolygonDataFrame。您可以自己获取数据：data(rase_data, package = 'rase')
我从坐标列表开始（lat/long 每个物种）。感谢这个答案
#dput(specieslist)
specieslist <- structure(list(Species = c("A", "A", "A", "A", "A", "M", "A", "M", "A", "A", "A", "A", "A", "A", "M", "M", "A", "M", "A", "A", "A", "M", "M", "M", "A", "A", "A", "A", "A", "A", "A", "M", "A", "A", "M", "M", "A", "M", "M", "A"), lat = c(37.407002, 35.65242, 33.33891, 37.65105, 38.90657, 39.06893, 34.53998, 38.18311, 37.40006, 35.65242, 34.53998, 33.33891, 37.65105, 38.90657, 38.18311, 39.06893, 36.252183, 40.32555, 39.575983, 39.73548, 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, -120.77931, -119.85962, -116.91402, -119.88513, -121.05138, -120.37691, -120.86546, -121.775867, -121.91209, -121.554167, -121.70126, -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))
通过围绕点创建外壳来根据 species/points 制作多边形：
#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
#或
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", "A", "A", "A", "M", "A", "M", "A", "A", "A", "A", "A", "A", "M", "M", "A", "M", "A", "A", "A", "M", "M", "M", "A", "A", "A", "A", "A", "A", "A", "M", "A", "A", "M", "M", "A", "M", "M", "A"), lat = c(37.407002, 35.65242, 33.33891, 37.65105, 38.90657, 39.06893, 34.53998, 38.18311, 37.40006, 35.65242, 34.53998, 33.33891, 37.65105, 38.90657, 38.18311, 39.06893, 36.252183, 40.32555, 39.575983, 39.73548, 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, -120.77931, -119.85962, -116.91402, -119.88513, -121.05138, -120.37691, -120.86546, -121.775867, -121.91209, -121.554167, -121.70126, -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
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