# In this document we provide code for a function to be loaded in R. 
# After loading the function type "grid.plus()" in the R console to run it with default parameters.
# The function will create a grid sampling design with additional random sampling units placed onto grid lines. 
# The user can modify several input variables specifying the sampling design. The details are explained below.  
# Most importantly, the user can specify input coordinates for their sampling area. 
# The function will then provide coordinates of sampling units within this area 
# according to a grid sampling design with additional random sampling units.
# In case no input is specified, the ouput will be an example sampling design.  
# By default the function will export a csv file with the sampling coordinates to the working directory
# This code requires the PBSmapping and sp packages.

grid.plus<-function(X=c(5.275,5.28,5.278,5.27), Y=c(53.27, 53.266,53.266,53.27), isd=25, frp=.25, frcell=0.5, plotit=T, save=T, output.filename="gridplus", zone=31){

	#input arguments
	# X,Y - LL coordinates bordering the sampling area (in decimal degrees)
	# isd - inter sample distance (in meters)
	# frp - fraction of random sampling units 
	# frcell - specifies the border around the sampling area as fraction of isd 
	# plotit - logical specifying if graphical output is required
	# save - logical specifying if exporting the coords to a csv file is required
	# output.filename - filenames of csv export file 
	# zone - utm zone used by the function convUL (from the library PBSmapping) to transform coordinates between LL and UTM 

	#output values
	# X,Y - LL coordinates of sampling units acording to grid sampling with additional random sampling untis placed onto grid lines 
	# point - type of sampling point (1: grid sampling station, and 2: random sampling station)
	# ID - unique identifier of each sampling station

	#load libraries
	require(PBSmapping,quietly=T) 
	require(sp,quietly=T)

	#calculate buffer around sampling area 
	buffer<-frcell*isd

	#create dataframe of coordinates
	d<-data.frame(
		X=X,
		Y=Y
		)

	#stop function if the input contains less than 3 sets of coordinates 
	if(dim(d)[1]<3){
		stop("specify at least 3 sets of XY-coordinates bordering the sampling area")
	}

	#transform LL coordinates to UTM
	attr(d, "projection")<- "LL"
	attr(d, "zone")<- zone
	d<-convUL(d, km=F)

	#create square bounding box 
	bbox<-rbind(
		c(min(d$X),min(d$Y)), 	#bottom left 
		c(min(d$X),max(d$Y)),	#top left
		c(max(d$X),max(d$Y)),	#top right
		c(max(d$X),min(d$Y))	#bottom right
		)

	#create convex hull 
	c.hull<-calcConvexHull(d)	

	### create grid sampling design ###
	x<-seq((min(bbox[,1])-buffer),(max(bbox[,1])+buffer),isd)
	y<-seq((min(bbox[,2])-buffer),(max(bbox[,2])+buffer),isd)
	GRID<-expand.grid(x,y)

	#remove sampling stations that fall outside of convex hull (sampling area) 
	idx<-point.in.polygon(GRID[,1],GRID[,2],c.hull$X,c.hull$Y)
	GRID<-subset(GRID,idx==1)
	GRID$point<-"1"
	names(GRID)<-c("X","Y","point")

	#subsample a fraction (i.e. frp) of grid sampling stations for adding random sampling stations
	n<-round(dim(GRID)[1]*frp) 
	subsample<-sample(c(1:dim(GRID)[1]),n,replace=F)
	subsample.coords<-GRID[,1:2][subsample,]

	#calculate distance offset (in X and Y) for the random sampling stations 
	#between half the inter sample distance around the grid sampling station 
	subsample.dist<-cbind(runif(n, min=-0.5*isd,0.5*isd),runif(n, min=-0.5*isd,0.5*isd))

	#define on which grid line (i.e. N-S, or E-W) the random sampling station will be located and set the other offset to 0
	subsample.gridline<-rbinom(n=n,size=1,prob=0.5)
	subsample.dist<-subsample.dist*(cbind(subsample.gridline,1-subsample.gridline)) 
	
	#finalize coordinates by adding the distance offset to the grid sampling station's coordinates
	subsample.coords<-subsample.coords+subsample.dist
	
	#create output
	subsample<-as.data.frame(cbind(subsample.coords[,1:2],rep("2",n)))
	names(subsample)<-c("X","Y","point")
	out<-rbind(GRID,subsample)
	out$ID<-1:dim(out)[1]

	#transform UTM coords back to LL
	attr(out, "projection")<- "UTM"
	attr(out, "zone")<- zone
	out<-convUL(out, km=F)

	#plot the sampling design
	if(plotit==T){
		plot(out$X,out$Y,pch=19,col=out$point, main=paste(dim(out)[1]," samples", sep=""),
		xlab="lattitude",ylab="longitude")
		}
	
	#save sampling design in 
	if(save==T){
		if(is.null(output.filename)){
			stop("No output filename specified")
			}
		write.table(out,paste(getwd(),"/",output.filename,".csv",sep=""),row.names = F,sep=",")
		}
	return(out)
	}