R/allocPoly.r

In BradHubley/SpatialHub: Survey Design and Spatial Analysis of Fisheries Data

#' @title allocPoly
#' @description ALLOCATION by STRATIFIED RANDOM DESIGN based on bottom type (or another polygon based non-continuous stratifying variable)
#' @param poly.lst = list containing PBSmapping::PolySet and PBSmapping::PolyData that describe strata polygons, PolyData may contain allocation to specifiy the number of station to be selected for each strata and repeats to indicate how many stations should be selected from repeated.tows
#' @param bounding.poly = PBSmapping::PolySet describing the area to be surveyed
#' @param ntows = total number of stations to be selected
#' @param mindist = minimum distance (km) or buffer between stations
#' @param pool.size = multiplyer to n.tows. This function first generates a pool of random locations, then randomly selects from the pool stations that fulfill the criteria (allocation by strata, buffer distance, etc.)
#' @param repeated.tows = PBSmapping::EventData of previous stations to be selected from if repeats exist in poly.lst[[2]] (PolyData)
#' @param map = preset location to pass to bioMap
#' @param lplace = placement of legend if !is.null(map)
#' @param show.pool = logical, if TRUE the initial pool of random stations is plotted on the map
#' @param UTMzone = PBSmapping::PolySet attribute, auto-generated if missing
#' @importFrom graphics legend
#' @importFrom stats na.omit
#' @importFrom PBSmapping calcArea
#' @importFrom PBSmapping combinePolys
#' @importFrom PBSmapping convUL
#' @importFrom spatstat.geom nndist
#' @importFrom spatstat.geom as.ppp
#' @importFrom spatstat.geom owin
#' @author Brad Hubley
#' @export


allocPoly<-function(poly.lst,bounding.poly,ntows,mindist=1,pool.size=4,repeated.tows=NULL, map=NULL,lplace='bottomleft',show.pool=F,UTMzone){

	options(warn=-1)
	# create pool of random points
	if(missing(ntows))ntows<-sum(poly.lst[[2]]$allocation)
	npool=ntows*pool.size

	if(!missing(bounding.poly))surveyed.polys<-joinPolys(poly.lst[[1]],bounding.poly,operation="INT")

	if(missing(bounding.poly)){
	surveyed.polys<-poly.lst[[1]]
	bounding.poly<-poly.lst[[1]][chull(poly.lst[[1]]$X,poly.lst[[1]]$Y),]
		bounding.poly$POS<-1:nrow(bounding.poly)
		bounding.poly$SID<-1
	}
	attr(bounding.poly,"projection")<-"LL"
	if(!missing(UTMzone))attr(bounding.poly,"zone")<-UTMzone

	# start with a pool of random stations
	pool.EventData<-genran(npool,bounding.poly,mindist=mindist)


	Poly.ID<-unique(poly.lst[[2]]$PID)
	strata<-as.character(unique(poly.lst[[2]]$PName))
	strataTows.lst<-list(NULL)

	# allocation provided
	if("allocation"%in%names(poly.lst[[2]])){
		towsi<-with(poly.lst[[2]],tapply(allocation,PName,unique))
		strataPolys.dat<-merge(surveyed.polys,subset(poly.lst[[2]],select=c("PID","PName")))
		attr(strataPolys.dat,"projection")<-"LL"
		if(!missing(UTMzone))attr(strataPolys.dat,"zone")<-UTMzone
		strataArea<-calcArea(strataPolys.dat,1)
	}
	else{
		# calculate area and proportional allocation
		strataPolys.lst<-list(NULL)
		strataArea<-c()
		towsi<-c()

		for(i in 1:length(strata)){
			strataPIDS<-poly.lst[[2]]$PID[poly.lst[[2]]$PName==strata[i]]
			tmp<-subset(surveyed.polys,PID%in%strataPIDS)
			if(nrow(tmp)>0){
				tmp$PName<-strata[i]
				strataPolys.lst[[i]]<-combinePolys(tmp)
				attr(strataPolys.lst[[i]],"projection")<-"LL"
				if(!missing(UTMzone))attr(strataPolys.lst[[i]],"zone")<-UTMzone
				strataArea[i]<-calcArea(strataPolys.lst[[i]],1)$area
				names(strataArea)[i]<-strata[i]
				print(strata[i])
			}
		}
		strataArea<-na.omit(strataArea)
		strataPolys.dat<-do.call("rbind",strataPolys.lst)
		towsi<-round(strataArea/sum(strataArea)*ntows)
		towsi<-towsi[towsi>0]
		towsi[1]<-ntows-sum(towsi[-1])
		strata<-names(towsi)
	}
	browser()
	for(i in 1:length(strata)){
		LocSet<-findPolys(pool.EventData,subset(strataPolys.dat,PName==strata[i]))
		strataTows.lst[[i]]<-data.frame(subset(pool.EventData,EID%in%LocSet$EID),Poly.ID=Poly.ID[i],STRATA=strata[i])

	}

	Tows<-do.call("rbind",strataTows.lst)
	Tows$EID<-1:nrow(Tows)
	rownames(Tows)<-1:nrow(Tows)
	attr(Tows,"projection")<-"LL"
	if(!missing(UTMzone))attr(Tows,"zone")<-UTMzone

	# randomly selects stations from last years	survey (repeated.tows)
	if(!is.null(repeated.tows)){

		names(repeated.tows)<-c("EID","X","Y","Poly.ID")
		repeated.tows$EID<-repeated.tows$EID+1000
		repeat.str<-poly.lst[[2]][!is.na(poly.lst[[2]]$repeats),]
		repeated.lst<-list(NULL)
		#browser()
		for(i in 1:length(repeat.str$PID)){
			str.tows<-subset(repeated.tows,Poly.ID==repeat.str$PID[i])
			repeated.lst[[i]]<-str.tows[sample(1:nrow(str.tows),repeat.str$repeats[repeat.str$PID==repeat.str$PID[i]]),]
			repeated.lst[[i]]$STRATA<-repeat.str$PName[repeat.str$PID==repeat.str$PID[i]]
		}
		repeated.tows<-do.call("rbind",repeated.lst)
		tmp<-rbind(Tows[,-4],repeated.tows)
		attr(tmp,"projection")<-"LL"
		if(!missing(UTMzone))attr(tmp,"zone")<-UTMzone
		tmp$nndist<-nndist(as.ppp(subset(convUL(tmp),select=c('X','Y')),with(convUL(bounding.poly),owin(range(X),range(Y)))))
		Tows<-subset(tmp,nndist>mindist&EID<1000)

	}

	Tows.lst<-list(NULL)
	strata<-as.character(unique(poly.lst[[2]]$PName))
	#browser()
	for(i in 1:length(strata)){
		Tows.lst[[i]]<-subset(Tows,STRATA==strata[i])[1:towsi[strata[i]],]

	}
	Tows<-do.call("rbind",Tows.lst)
	Tows$EID<-1:nrow(Tows)
	rownames(Tows)<-1:nrow(Tows)
	attr(Tows,"projection")<-"LL"
	if(!missing(UTMzone))attr(Tows,"zone")<-UTMzone

	if(!is.null(repeated.tows))Tows<-list(new.tows=Tows, repeated.tows=repeated.tows)

	if(!is.null(map)){
		bioMap(map,poly.lst=list(surveyed.polys,poly.lst[[2]]))
		bg.col<-tapply(poly.lst[[2]]$col,poly.lst[[2]]$PName,unique)
		if(is.null(repeated.tows))addPoints(Tows,pch=21, cex=1,bg=bg.col[as.character(Tows$STRATA)])
		if(!is.null(repeated.tows)){
			addPoints(Tows$new.tows,pch=21, cex=1,bg=bg.col[as.character(Tows$new.tows$STRATA)])
			addPoints(Tows$repeated.tows,pch=24, cex=1,bg=bg.col[as.character(Tows$repeated.tows$STRATA[order(Tows$repeated.tows$EID)])])
		}
		if(show.pool)addPoints(pool.EventData,pch=4,cex=0.4)
#		browser()
		if(!is.null(repeated.tows))legend(lplace,legend=names(bg.col[unique(as.character(Tows$new.tows$STRATA))]),pch=21,pt.bg=bg.col[unique(as.character(Tows$new.tows$STRATA))],bty='n',cex=1, inset = .02)
		if(is.null(repeated.tows))legend(lplace,legend=names(bg.col[unique(as.character(Tows$STRATA))]),pch=21,pt.bg=bg.col[unique(as.character(Tows$STRATA))],bty='n',cex=1, inset = .02)

	}

	options(warn=0)
	return(list(Tows=Tows,Areas=strataArea))

}