# File: OregonEstuary.R
# Purpose: Illustrate area sample for an estuarine resource
# Programmer: Tony Olsen
# Date: May 5, 2005

#######
Load SP library
#######
#######
# Load psurvey.design library
#######


# Read in dbf file as tab file
att <- read.dbf("Oregon2001Estuaries")
head(att)
tail(att)

# change calculated area of estuary polygons to sq km from sq meters
att$area_mdm <- att$area_mdm/1000/1000

# Deterine NCA size classes present
summary(att$CLASS)

# create an estuary size variable
att$EstSize <- cut(att$SQ_KILOMET, breaks=c(0, 1, 25,50), include.lowest=TRUE,
					labels=c('Small', 'Medium', 'Large') )

dsgn <- list("< 100 SQ.Kilometers"=list(
       panel=c(PanelOne=25, PanelTwo=25),
        				seltype="Unequal", 
        				caty.n=c(Small=17, Medium=17, Large=16), 
        				over=50) )

sample(100000000,1) # run to get random seed

dsgntime <- proc.time()
sites <- grts(design=dsgn, 
			DesignID="ORUNEQ04", 
			type.frame="area",
   			src.frame="shapefile",
   			in.shape='Oregon2001Estuaries',
   			att.frame=att,
   			stratum="CLASS", 
   			mdcaty="EstSize",  
   			shapefile=TRUE, 
   			prjfilename='Oregon2001Estuaries', 
   			out.shape="OREstuaryUnequalSites")

dsgntime <- proc.time() - dsgntime
dsgntime

#### end of design process.  Shapefile automatically written

# IF want to add Lat/Lon to sites file
# function to compute Deg:Min:Sec
dms <- function(dd) {
	x <- abs(dd)
	d <- trunc(x)
	m <- trunc((x - d) * 60)
	s <- (x - d - m/60) * 3600 
	d <- sign(dd) * d
	data.frame(deg=d, min=m, sec=s)
}
					
tmp <- albersgeod(sites$x, sites$y, sph='Clarke1866', 
				clon=-96, clat=23, sp1=29.5, sp2=45.5)
londms <- dms(tmp$lon)
latdms <- dms(tmp$lat)
nc <- ncol(sites)
sites <- data.frame(sites[,1:3],LonDD=tmp$lon,LatDD=tmp$lat,
			LonDeg=londms$deg, LonMin=londms$min, LonSec=londms$sec,
			LatDeg=latdms$deg, LatMin=latdms$min, LatSec=latdms$sec,
			sites[,4:ncol(sites) ] 
							)
# Write out design file
options(digits=11)
write.table(sites, 'OregonSites.csv', sep=',', row.names=FALSE)
options(digits=7)