I'm trying to build species distribution polygons for use in the R program rase. The program requires an owin object but sample data also includes a SpatialPolygonDataFrame. You can get the data yourself with: data(rase_data, package = 'rase')
I'm starting with a list of coordinates (lat/long per species). Thanks to this answer here, I've been able to make a polygon per element of the list (each species). I need to get to an owin object. Here's the dput of some test data and then code I've used to get where I'm at.
#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))
Make the polygon per species/points by creating hull around 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 ) )
Now I think df.sf (sf points object) becomes the SpatialPolygonDataFrame and hulls (sf polygon object) becomes an owin object:
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'
The problems are: df.sf_SPDF seems to be formatted incorrectly and Dist_owin errors.
I find all this spatial work in R very confusing. I've been working on this for several days now.
UPDATE: If I try another way- convert geometry to polygon and then make the owin. This produces an error:
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”
I do not know sf enough to fix this, so I show it via terra but the important part is the sequence of operations. You can implement that in sf again if you wish. There should be no need to revert to the old Spatial* objects
Your data
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))
First I make a spatial object, a SpatVector in this case, and I transform it to a planar CRS --- to get that out of the way.
Your choice of epsg:6345, that is +proj=utm +zone=16 is inappropriate for your data. Zone 16 is for the longitude of Alabama. California covers two UTM zones so you cannot use that. Instead use e.g. "Teale Albers" if all your data are confined to the Golden State.
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)
To simplify things, I show a workflow for 1 species
usp <- unique(v$Species)
sp <- v[v$Species==usp[1]]
Make a convex hull, and I think you would want to add a buffer.
ch <- terra::convexhull(sp)
bch <- buffer(ch, 25000)
plot(bch)
points(sp)
Now make the owin via sf
library(sf)
library(spatstat)
sfobj <- st_as_sf(bch)
owin <- as.owin(sfobj)
You can extract the points in new CRS like this
pxy <- terra::coords(sp)
And now create a spatstat ppp object.
x <- ppp(pxy[,1], pxy[,2], window=owin)
#Warning message:
#data contain duplicated points
To avoid the above warning, you could use, at the start of the script, specieslist <- unique(specieslist)
x
#Planar point pattern: 27 points
#window: polygonal boundary
#enclosing rectangle: [-222286.97, 312378.62] x [-539742.6, 217425] units