Why vectorize the data when you can get at this, much more efficiently, in a raster data structure. It is always easier, and less computationally expensive, to overlay cells rather than vector topology. This also removes the spatial uncertainty of converting an n-resolution raster to a vector and then intersecting the results. You could always convert the final raster overlay results to polygons.
First, lets create some data emulating what you have described. I am using "raster.gaussian.smooth" to create a correlated "patchy" landscape and then stretching back to a 0-10 scale thus, the extra steps.
library(raster)
library(spatialEco)
s = 3 # Gaussian sigma parameter for raster.gaussian.smooth function
w = 21 # Matrix (window) size for raster.gaussian.smooth function
s1.summer <- raster(nrows=500, ncols=500)
s1.summer[] <- round(runif(ncell(s1.summer), 0,10),0)
s1.summer <- raster.gaussian.smooth(s1.summer, sigma = s, n = w)
s1.summer <- raster.transformation(s1.summer, "stretch", 0, 10)
s2.summer <- raster(nrows=500, ncols=500)
s2.summer[] <- round(runif(ncell(s2.summer), 0,10),0)
s2.summer <- raster.gaussian.smooth(s2.summer, sigma = s, n = w)
s2.summer <- raster.transformation(s2.summer, "stretch", 0, 10)
s1.winter <- raster(nrows=500, ncols=500)
s1.winter[] <- round(runif(ncell(s1.winter), 0,10),0)
s1.winter <- raster.gaussian.smooth(s1.winter, sigma = s, n = w)
s1.winter <- raster.transformation(s1.winter, "stretch", 0, 10)
s2.winter <- raster(nrows=500, ncols=500)
s2.winter[] <- round(runif(ncell(s2.winter), 0,10),0)
s2.winter <- raster.gaussian.smooth(s2.winter, sigma = s, n = w)
s2.winter <- raster.transformation(s2.winter, "stretch", 0, 10)
par(mfrow=c(2,2))
plot(s1.summer)
plot(s2.summer)
plot(s1.winter)
plot(s2.winter)
We can now reclassify each raster to [0,1] based on the occupancy threshold and then sum the results thus, illustrating seasonal overlap.
occ <- 5 # occupancy threshold
s1.summer <- setValues(s1.summer, ifelse(getValues(s1.summer) >= occ, 1, 0) )
s2.summer <- setValues(s2.summer, ifelse(getValues(s2.summer) >= occ, 1, 0) )
s1.winter <- setValues(s1.winter, ifelse(getValues(s1.winter) >= occ, 1, 0) )
s2.winter <- setValues(s2.winter, ifelse(getValues(s2.winter) >= occ, 1, 0) )
summer.overlap <- s1.summer + s2.summer
winter.overlap <- s1.winter + s2.winter
par(mfrow=c(2,1))
plot(summer.overlap, main="summer species overlap")
plot(winter.overlap, main="winter species overlap")
The resulting summer and winter overlap rasters could then be vectorized based on the summed values or, by reclassifying to [0,1]. However, even in treating the overlap as a binomial process, I would not expect a single polygon.