Geospatial analysis - Thessian Polygon in R

thessian_clip.r

require(spstat)

point.dat;# a shapefile of points with associated data.
study.area;# a shapefile of the study area/sampling area.

# Generate Thessian polygon and assign CRS
th=as(dirichlet(as.ppp(point.dat)),"SpatialPolygons")
proj4string(th)=proj4string(point.dat)

# Join thessian polygon with actual data
th.z=over(th,point.dat)

# Convert to a spatial polygon
th.spdf=SpatialPolygonsDataFrame(th,th.z)

# Clip to study area
th.clp=raster::intersect(study.area,th.spdf)

thessian_clip_alt.r

require(dismo)
require(sp)
require(raster)
require(rgeos)
  
point.dat;# a shapefile of points with associated data.
study.area;# a shapefile of the study area/sampling area.

bbox.da=c(bbox(study.area)[1,1],bbox(study.area)[1,2],bbox(study.area)[2,1],bbox(study.area)[2,2])

th=dismo::voronoi(point.dat,ext=bbox.da)
th.z=sp::over(th,point.dat)
th.z.spdf=sp::SpatialPolygonsDataFrame(th,th.z)
th.clp=raster::intersect(study.area,th.z.spdf)
th.clp$area_sqkm=rgeos::gArea(th.clp,byid=T)*1e-6;#Assumes coordinate system in meters

## combining this all into a function
thessian_create=function(data,clip,plot=T,full.extent=T){
  require(dismo)
  require(sp)
  require(raster)
  require(rgeos)
  if(full.extent==T){bbox.da=c(bbox(clip)[1,1],bbox(clip)[1,2],bbox(clip)[2,1],bbox(clip)[2,2])}else{bbox.da=NULL}
  th=voronoi(data,ext=bbox.da)
  th.z=sp::over(th,data)
  th.z.spdf=sp::SpatialPolygonsDataFrame(th,th.z)
  th.clp=raster::intersect(clip,th.z.spdf)
  th.clp$area_sqkm=rgeos::gArea(th.clp,byid=T)*1e-6
  if(plot==T){plot(th.clp)}
  return(th.clp)
}