Skip to content

Instantly share code, notes, and snippets.

@SwampThingPaul
Last active January 7, 2019 14:37
Show Gist options
  • Save SwampThingPaul/0398073092ba337bf86862a6e1d6d45d to your computer and use it in GitHub Desktop.
Save SwampThingPaul/0398073092ba337bf86862a6e1d6d45d to your computer and use it in GitHub Desktop.
Geospatial analysis - Thessian Polygon in 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)
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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment