# File: adjustwts.txt
# Purpose: 1999 + 2000 survey design for Estuaries: adjust 
#  weights using Design File, frame information, and Site Evaluation file
# Programmer: Tony Olsen
# Date: January 2006
#
#Example updated for R 2.2.1 and psurvey.analysis 2.9
#
# Frame Area divided into two strata: Tidal Creeks and Open water
#  Tidal Creeks, identified as < 100 m wide on GIS cover
#  Open Water areas, identified as > 100 m wide on GIS cover, represent 
#  larger estuarine rivers and sounds
#   Area of estuaries based on GIS coverages of NWI maps.
#  The design documentation has  105.829522 sq km of tidal creeks;
#  628.509298 sq km of open water estuaries.
framearea <- c('Open Water'=628.509298, 'Tidal Creek'=105.829522 )

# Import Design File for Estuaries - Tidal Creeks + Open Water
DesignFile <- read.delim('original_data\\est_design.txt')
head(DesignFile)


# What are the stratum and panel present?
table(DesignFile$stratum)
table(DesignFile$panel)
table(DesignFile$stratum, DesignFile$Type) # can use Type as stratum

# Import Open Water Site Evaluation file
SiteEval <- read.delim('original_data\\siteeval.txt')
tail(SiteEval)

# Check codes for Sampling Status
table(SiteEval$EE.Status)


# Note on code meanings
# TS - target and sampled site
# NT - site was non-target (dry, too shallow,...)
# NN - site not needed (not evaluated)

# Create DesignStatus file
DesignStatus <- merge(DesignFile, SiteEval, 
						by.x='Site.id', by.y='EE.Site.ID')

# Add equal area x,y coordinates to be used in variance estimation
# Note that longitude must be expressed as a negative
tmp <- marinus(DesignStatus$lat.dd,DesignStatus$long.dd)
DesignStatus$xmarinus <- tmp[,'x']
DesignStatus$ymarinus <- tmp[,'y']

# check to see status codes & not use NN if present
table(DesignStatus$EE.Status)
tst <- DesignStatus$EE.Status == 'NN'
DesignStatus <- DesignStatus[!tst,]

# Adjust weights to match frame area
# Determine which sites to include in weight adjustment.
# Note that original design weights assumed sample size of 30; 
sites <- rep(TRUE,nrow(DesignStatus))
DesignStatus$final.wt <- adjwgt(sites,DesignStatus$Weight ,
			DesignStatus$Type,framearea)
# Check out weights
tapply(DesignStatus$final.wt, DesignStatus$EE.Type, sum)

# Write out Design Status File
write.table(DesignStatus, file="designstatus.csv",sep=",",row.names=FALSE)