#' Function to generate Thiessen polygons from SILO sites
#'
#'@param SILOdata - data loaded from SILO based on site list
#'@param path - file path to save Thiessen polygon shapefile
#'@param shpname - filename to save ESRI shapefile (no extension)
#'@param boundary - optional either a filename(including path) of a boundary, e.g. catchment boundary, to apply, or the boundary as a sfc_MULTIPOLYGON object
#'
#'@return A simple feature geometry (sf::sfc object) of the polgyons created. Shape file saved to path \\ shpname
#'
#' If boundary is specified weights are written to the attribute table of the polygon return, which can be extracted with
#' \code{st_drop_geometry(returnedfeature[c("Station","weights")])}
#'
#'@examples
#' \dontrun{
#'X<-SILOLoad(c("24001","24002","24003"),path="./SWTools/extdata")
#'p<-SILOThiessenShp(X,tempdir(),"Theissens")
#'a<-SILOSiteSummary(X)
#'ggplot(p)+geom_sf(aes(fill=AnnualRainfall))+
#'geom_point(data=a,aes(Longitude,Latitude))+
#'geom_text(data=a,aes(Longitude,Latitude,label=Site),nudge_y = 0.02)
#'}
#'@export
SILOThiessenShp<- function(SILOdata,path,shpname,boundary=NULL){
SiteTable <- SILOSiteSummary(SILOdata) #table summarising SILO sites
SILOLocs <- sf::st_as_sf(SiteTable, coords = c("Longitude","Latitude"))
if(is.character(boundary))
{
area<-sf::st_read(boundary,stringsAsFactors=FALSE)
area<-area$geometry
}else if(is.null(boundary))
{
area<-sf::st_polygon() #empty polygon used for boundary of Thiessen
}
if(!inherits(area,"sfc_MULTIPOLYGON")) stop(paste("Boundary not created, possibly couldn't read in",boundary))
if(length(area)>1) stop(paste(boundary,"has",nrow(area),"features, only 1 supported"))
. <- NULL #address global variable . warning. bit of a hack.
# compute Voronoi (Thiessen) polygons
TPoly <- SILOLocs %>%
sf::st_geometry() %>%
do.call(c, .) %>%
sf::st_voronoi(area) %>%
sf::st_collection_extract()
# match polygons to points for combining with site data table
SILOLocs$polys <- sf::st_intersects(SILOLocs, TPoly) %>% unlist %>% TPoly[.]
TPoly_id <- sf::st_set_geometry(SILOLocs, "polys")#["Station"] #Attaching station ID to polygons for export
TPolyCRS <- sf::st_set_crs(TPoly_id, "+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs") #Defining projection - GDA94
if(!is.null(boundary))
{
area<-sf::st_transform(area,sf::st_crs(TPolyCRS))
TpolyCRS<-suppressMessages(sf::st_intersection(TPolyCRS,area))
TpolyCRS$AlbersArea<-sf::st_area(TpolyCRS)
TpolyCRS$weights<-as.numeric(TpolyCRS$AlbersArea/sum(TpolyCRS$AlbersArea))
}
return(TpolyCRS)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.