# File: NELakesDesign.R
# Purpose: Example GRTS survey designs for New England lakes.
#					Lake locations are given by their centroid.
#					Data originally from National Hydrology Database
# Programmer: Tony Olsen
# Date: April 17, 2005
#
#######
# Load psurvey.design library
#######


#######
# Read in attribute table from shapefile saved as tab delimited file
att <- read.dbf("reg1_lakes")
names(att) <- tolower(names(att))  # I like lower case names
tail(att)

# Most column names are self-explanatory.
# lakres_alb_ - unique ID for each lake
# ftype - identifies as lake or reservoir
# sq_km - lake area in sq km
# x_coord - Albers projection x-coordinate
# y_coord - Albers projection y-coordinate
# level3 - Omernik Level 3 ecoregion
# level2 - Omernik Level 2 ecoregion
# level1 - Omernik Level 1 ecoregion
# change lake area to hectares
att$area.ha <- att$sq_km*100
summary(att$area.ha)
hist(att$area.ha, breaks=100)
hist(log10(att$area.ha), breaks=100)

##### This section is not necessary to do designs. Informational only.
#### Check lake type, spatial pattern and lake area
table(att$ftype)
table(att$st)
table(att$level1)
table(att$level2)
table(att$level3_nam)


#####  Four example designs given below.  Each can be run independently.

##### Equal Probability GRTS design for all lakes and reservoirs
Equaldsgn <- list(None=list(panel=c(PanelOne=300),
	            seltype='Equal'))

### select sample
sample(10000000,1) # run to get random seed
dsgntime <- proc.time()

Equalsites <- grts(design=Equaldsgn,
                   DesignID='EQUAL',
                   type.frame='finite',
                   src.frame='shapefile',
                   in.shape='reg1_lakes',
                   att.frame=att,
                   prjfilename='reg1_lakes',
                   out.shape='Equal Sites'
                   )
dsgntime <- ( proc.time() - dsgntime )/60  # time in minutes to complete design
dsgntime

# Write out design file
options(digits=11)   # make sure have sufficient precision for lake location
write.table(Equalsites,'EqualSites.tab',sep='\t',row.names=FALSE)
options(digits=7)


##### Stratified GRTS design for all lakes and reservoirs
# stratify by state
table(att$st)
Stratdsgn <- list(CT=list(panel=c(PanelOne=50),
			  seltype='Equal'),
		  MA=list(panel=c(PanelOne=50),
			  seltype='Equal'),
	  	  ME=list(panel=c(PanelOne=50),
			  seltype='Equal'),
		  NH=list(panel=c(PanelOne=50),
			  seltype='Equal'),
		  RI=list(panel=c(PanelOne=50),
			  seltype='Equal'),
		  VT=list(panel=c(PanelOne=50),
		          seltype='Equal')
		  	   )

### select sample
sample(10000000,1) # run to get random seed
dsgntime <- proc.time()

Stratsites <- grts(design=Stratdsgn,
                   DesignID='EQUAL',
                   type.frame='finite',
                   src.frame='shapefile',
                   in.shape='reg1_lakes',
                   att.frame=att,
                   stratum='st',
                   prjfilename='reg1_lakes',
                   out.shape='Stratified Sites'
                   )
dsgntime <- ( proc.time() - dsgntime )/60  # time in minutes to complete design
dsgntime

# Write out design file
options(digits=11)   # make sure have sufficient precision for lake location
write.table(Stratsites,'StratSites.tab',sep='\t',row.names=FALSE)
options(digits=7)

#### Site Summaries
addmargins( table(Stratsites$stratum,Stratsites$mdcaty) )




######  Unequal probability based on lake area and oversample
### Calculate lake area categories
att$lake.area.cat <- cut(att$area.ha, breaks=c(0,1,5,10,50,500,70000),
				           include.lowest=TRUE)
summary(att$lake.area.cat)

Areadsgn <- list(None=list(panel=c(PanelOne=300),
                 seltype='Unequal',
		     caty.n=c('[0,1]'=50, '(1,5]'=50, '(5,10]'=50, '(10,50]'=50,
		         '(50,500]'=50, '(500,7e+04]'=50 ),
			over=120)
		  	     )

### select sample
sample(10000000,1) # run to get random seed
dsgntime <- proc.time()

UnequalSites <- grts(design=Areadsgn,
                   DesignID='UNEQUAL',
                   type.frame='finite',
                   src.frame='shapefile',
                   in.shape='reg1_lakes',
                   att.frame=att,
                   mdcaty='lake.area.cat',
                   prjfilename='reg1_lakes',
                   out.shape='Unequal Sites'
                   )
dsgntime <- ( proc.time() - dsgntime )/60  # time in minutes to complete design
dsgntime

# Write out design file
options(digits=11)   # make sure have sufficient precision for lake location
write.table(Areasites,'AreaSites.tab',sep='\t',row.names=FALSE)
options(digits=7)

#### Site Summaries
addmargins( table(UnequalSites$mdcaty, UnequalSites$panel) )





######  Unequal probability based on lake area  with oversample
######  Add panel structure for survey over time
### Calculate lake area categories
att$lake.area.cat <- cut(att$area.ha, breaks=c(0,1,5,10,50,500,70000),
				           include.lowest=TRUE)
summary(att$lake.area.cat)

Paneldsgn <- list(None=list(panel=c(Annual=50, Year1=50, Year2=50, Year3=50,
				Year4=50, Year5=50),
                     seltype='Unequal',
			   caty.n=c('[0,1]'=50, '(1,5]'=50, '(5,10]'=50, '(10,50]'=50,
			   		'(50,500]'=50, '(500,7e+04]'=50 ),
			   over=120)
		  	     )

### select sample
sample(10000000,1) # run to get random seed
dsgntime <- proc.time()

Panelsites <- grts(design=Paneldsgn,
                   DesignID='UNEQUAL',
                   type.frame='finite',
                   src.frame='shapefile',
                   in.shape='reg1_lakes',
                   att.frame=att,
                   mdcaty='lake.area.cat',
                   prjfilename='reg1_lakes',
                   out.shape='Panel Sites'
                   )
dsgntime <- ( proc.time() - dsgntime )/60  # time in minutes to complete design
dsgntime

# Write out design file
options(digits=11)   # make sure have sufficient precision for lake location
write.table(Panelsites,'PanelSites.tab',sep='\t',row.names=FALSE)
options(digits=7)

#### Site Summaries
addmargins(table(Panelsites$mdcaty,Panelsites$panel ))