################################################################################
# Functions to:
#   1. extract estuary information and NWI for user points
#   2. create kml files to view data in GoogleEarth
#
#
# Created:  Lee McCoy, Melanie Frazier & Deborah A. Reusser
# Date:  01/30/2012
###############################################################################S

print("estuRe v. 1")
print("Loading necessary packages")

# ---------------------------------------------------------------------
## load required libraries

if(!("package:rgdal" %in% search())) require(rgdal)
if(!("package:maptools" %in% search())) require(maptools)

# ---------------------------------------------------------------------
## merge function that preserves order:
merge.with.order <- function(x,y, ..., sort = T, keep_order)
	{
		# this function works just like merge, only that it adds the option to return the merged data.frame ordered by x (1) or by y (2)
		add.id.column.to.data <- function(DATA)
		{
			data.frame(DATA, id... = seq_len(nrow(DATA)))
		}
		# add.id.column.to.data(data.frame(x = rnorm(5), x2 = rnorm(5)))
		order.by.id...and.remove.it <- function(DATA)
		{
			# gets in a data.frame with the "id..." column.  Orders by it and returns it
			if(!any(colnames(DATA)=="id...")) stop("The function order.by.id...and.remove.it only works with data.frame objects which includes the 'id...' order column")
 
			ss_r <- order(DATA$id...)
			ss_c <- colnames(DATA) != "id..."
			DATA[ss_r, ss_c]
		}
 
		# tmp <- function(x) x==1; 1	# why we must check what to do if it is missing or not...
		# tmp()
 
		if(!missing(keep_order))
		{
			if(keep_order == 1) return(order.by.id...and.remove.it(merge(x=add.id.column.to.data(x),y=y,..., sort = FALSE)))
			if(keep_order == 2) return(order.by.id...and.remove.it(merge(x=x,y=add.id.column.to.data(y),..., sort = FALSE)))
			# if you didn't get "return" by now - issue a warning.
			warning("The function merge.with.order only accepts NULL/1/2 values for the keep_order variable")
		} else {return(merge(x=x,y=y,..., sort = sort))}
	}

# ---------------------------------------------------------------------
## Function 1: Saves 2 or 3 data files and creates an object that can be
## used in subsequent manipulations
points2est <- function(x, lon, lat, climate=TRUE, FileName="NewData", saveDir){

#FileName = "CDA"
#setwd("C:\\Documents and Settings\\mfrazi03\\Desktop\\WestuRe_v1_master")
#x=read.csv("ExampleData.csv")
#lon="x"
#lat="y"
#climate=TRUE
#saveDir="C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReV5.12\\SampleOutput"


if(sum(x[,lon]>0)) stop("Longitude coordinates should be negative")
if(!is.numeric(x[,lon])) stop("Longitude contains non-numeric data")
if(!is.numeric(x[,lat])) stop("Latitude contains non-numeric data")

DataDim <- dim(x)[1]
x <- x[!(is.na(x[,lon]) | is.na(x[,lat])),]
MissingLatLon <- DataDim - dim(x)[1]
if(MissingLatLon>0) warning(paste(MissingLatLon, "data points removed due to missing latitude/longitude", sep=" "))

# loading data:
print("Wait for it: preparing your data")
salinity <- readOGR(dsn=".\\Data\\NOAAsalinity", layer="SalinityZonesNOAA", verbose=FALSE)
watershed <- readOGR(dsn=".\\Data\\Watershed", layer="Watershed", verbose=FALSE)
estuaryData <- read.csv(".\\Data\\EstuaryData.csv", header=TRUE)
watershedData <- read.csv(".\\Data\\WatershedData.csv", header=TRUE)
climateData <- read.csv(".\\Data\\ClimateMonthly.csv", header=TRUE)

  # Not all estuaries have a unique watershed (due to NWI updates).
  # at this point, I will remove these and then add a note in data
CDAests <- estuaryData$wsID[duplicated(estuaryData$wsID)]
estuaryData <- estuaryData[!duplicated(estuaryData$wsID),]  #keeps only the first instance of duplicated data, this prevents problems in subsequent merges.


# converting user data to spatial points dataframe:
  coordinates(x) <- as.formula(paste("~",lon,"+",lat,sep=""))
  proj4string(x) <-  CRS(proj4string(salinity))

  #-------------------------------------------------------------------------
# associate user sample points with NOAA salinity zone data
x@data$SalinityZone <- over(x, salinity)$SZ_NAME

# associate user sample points with watershed databases
x@data$wsID <- over(x, watershed)$wsID
dataEstWS <- estuaryData[estuaryData$wsID %in% x@data$wsID,]
# special cases
if(sum(dataEstWS$wsID %in% c("P090a", "P090w")>0)) {
    dataEstWS <- rbind(dataEstWS, estuaryData[estuaryData$wsID == "P090aXP090w",])
}
if(sum(dataEstWS$wsID %in% c("P130x02", "P130x02"))>0) {
    dataEstWS <- rbind(dataEstWS, estuaryData[estuaryData$wsID == "P130x01XP130x02",])
}

watershedData <- watershedData[, -which(names(watershedData)== "estCount")]
dataEstWS <- merge(dataEstWS, watershedData, by="wsID", all.x=TRUE)



# some notes for when a sample point is located in an estuary with a CDA for watershed
dataEstWS$wsNotes <- as.character(dataEstWS$wsNotes)
dataEstWS$wsNotes[dataEstWS$wsType=="CDA"] <- "Est. in CDA, watershed data may not be meaningful for estuary"
dataEstWS$estLatitude[dataEstWS$wsID %in% CDAests] <- NA
dataEstWS$estLongitude[dataEstWS$wsID %in% CDAests] <- NA

#####################################################
### Obtaining NWI data
#####################################################

# downloading the appropriate NWI shape files
if(sum(dataEstWS$estState=="CA", na.rm=TRUE)>0) {
    CA_NWI <- readOGR(dsn=".\\Data\\NWI_CA", layer="NWI_CA", verbose=FALSE)
    CA_NWI <- spChFIDs(CA_NWI, paste("CA", sapply(slot(CA_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

if((sum(dataEstWS$estState=="OR", na.rm=TRUE)>0)|(sum(dataEstWS$estState=="WA/OR"))) {
OR_NWI <- readOGR(dsn=".\\Data\\NWI_OR", layer="NWI_OR", verbose=FALSE)
OR_NWI <- spChFIDs(OR_NWI, paste("OR", sapply(slot(OR_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

if(sum(dataEstWS$estState=="WA", na.rm=TRUE)>0) {
    WA_NWI <- readOGR(dsn=".\\Data\\NWI_WA", layer="NWI_WA", verbose=FALSE)
    WA_NWI <- spChFIDs(WA_NWI, paste("WA", sapply(slot(WA_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

NWIfiles <- ls()[grep("_NWI", ls())]
if(length(NWIfiles)==0) warning("No data is located in NWI estuaries")
if(length(NWIfiles)==1) NWI <- get(NWIfiles[1])
if(length(NWIfiles)==2) NWI <- rbind(get(NWIfiles[1]), get(NWIfiles[2]))
if(length(NWIfiles)==3) NWI <- rbind(get(NWIfiles[1]), get(NWIfiles[2]), get(NWIfiles[3]))

# sapply(slot(WA..NWI, "polygons"), function(x) slot(x, "ID"))
x@data$NWIcode <- over(x, NWI)$NWIcode
x@data$NWIshort <- substring(x@data$NWIcode, 1, 2)
x@data <- merge.with.order(x@data, dataEstWS[,1:2], by="wsID", all.x=TRUE, sort=FALSE, keep_order=1)

  # Determining whether points are located outside of watershed data
  # This provides a warning.
  # NOTE: modified code from raster package
    bb <- bbox(watershed)
    p <- rbind(c(bb[1,1], bb[2,1]), c(bb[1,1], bb[2,2]), c(bb[1,2],
        bb[2,2]), c(bb[1,2], bb[2,1]), c(bb[1,1], bb[2,1]))
    bbPoly <- SpatialPolygons(list(Polygons(list(Polygon(p)), 1)))
    proj4string(bbPoly) <- CRS(proj4string(watershed))
    inPoly <- over(x, bbPoly)
    bbPolyNA <- sum(is.na(inPoly))
    if(bbPolyNA>0) warning(paste(bbPolyNA, "data points appear to be outside the geographic extent, you may want to check this!", sep=" "))


#####################################################
### Obtaining climate data
#####################################################

climateData <- climateData[(climateData$wsID %in% dataEstWS$wsID), ] 
climateData$wsID <- as.factor(as.character(climateData$wsID))

#PPT data
wsPPT <- aggregate(wsPPT ~ wsID, data=climateData, FUN=sum, na.rm=TRUE)
names(wsPPT)[2] <- "wsPPT_cumulative"
dataEstWS <- merge(dataEstWS, wsPPT, by="wsID", all.x=TRUE)
dataEstWS$AreaNormFreshFlow <- round((dataEstWS$wsArea*dataEstWS$wsPPT_cumulative * 0.001)/dataEstWS$estArea, 1)

#wsTemp data
wsTemp <- aggregate(wsTair_avg ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
wsTemp$wsTair_avg <- round(wsTemp$wsTair_avg, 2)
names(wsTemp)[2] <- "wsTair_YearMean"
dataEstWS <- merge(dataEstWS, wsTemp, by="wsID", all.x=TRUE)

#estTemp data
estTemp <- aggregate(estTair_avg ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
estTemp$estTair_avg <- round(estTemp$estTair_avg, 2)
names(estTemp)[2] <- "estTair_YearMean"
dataEstWS <- merge(dataEstWS, estTemp, by="wsID", all.x=TRUE)

#SSTemp data
SSTemp <- aggregate(estSST ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
SSTemp$estSST <- round(SSTemp$estSST, 2)
names(SSTemp)[2] <- "estSST_YearMean"
dataEstWS <- merge(dataEstWS, SSTemp, by="wsID", all.x=TRUE)


if(climate){
climateData <- merge(dataEstWS[,1:2], climateData, by="wsID", all.y=TRUE)
write.csv(climateData, paste(saveDir, paste(FileName, "_ClimateMonthly", ".csv", sep=""),  sep="\\"), row.names=FALSE)
}


# save csv files
write.csv(data.frame(x), paste(saveDir, paste(FileName, "_Samples", ".csv", sep=""),  sep="\\"), row.names=FALSE)
write.csv(dataEstWS, paste(saveDir, paste(FileName, "_Estuary", ".csv", sep=""),  sep="\\"), row.names=FALSE)

# function object
watershed <- watershed[watershed$wsID %in% dataEstWS$wsID,]
NWI <- NWI[NWI$wsID %in% dataEstWS$wsID,]
NWI <- NWI[!(is.na(NWI$wsID)),]
list(data=x, ests=dataEstWS, climate=climateData, watershedPoly = watershed, NWIPoly = NWI)
}

##testing function
#setwd("C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan302012")
#fakedata <- read.csv("ExampleData.csv")
#trial <- points2est(fakedata, lon="x", lat="y", climate=TRUE, FileName="Test",
#saveDir="C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan302012\\SampleOutput")

#summary(trial)

# ---------------------------------------------------------------------
## Function 1.5: Same function as points2est, but uses estuary name to 
##                obtain data.
## used in subsequent manipulations

names2est <- function(x, climate=TRUE, FileName="NewData", saveDir){

#x <- c("Big Creek Estuary, N Yachats", "Yaquina Estuary")  #testing
#climate=TRUE
#FileName="NewData"
#saveDir="C:\\Documents and Settings\\mfrazi03\\Desktop\\WestuRe_v1\\SampleOutput"

estuaryData <- read.csv(".\\Data\\EstuaryData.csv", header=TRUE)
ests <- toupper(x)
location <- c()
for(i in 1:length(ests)){
location <- c(location, grep(ests[i], toupper(estuaryData$estName), fixed=TRUE))
}
estuaryData <- estuaryData[location,]
lon <- "estLongitude"
lat <- "estLatitude"

# loading data:
print("Wait for it: preparing your data")
watershed <- readOGR(dsn=".\\Data\\Watershed", layer="Watershed", verbose=FALSE)
watershedData <- read.csv(".\\Data\\WatershedData.csv", header=TRUE)
climateData <- read.csv(".\\Data\\ClimateMonthly.csv", header=TRUE)

# converting user data to spatial points dataframe:
  coordinates(estuaryData) <- as.formula(paste("~",lon,"+",lat,sep=""))
  proj4string(estuaryData) <-  CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

dataEstWS <- merge(estuaryData, watershedData, by="wsID", all.x=TRUE)


#####################################################
### Obtaining NWI data
#####################################################

# downloading the appropriate NWI shape files
if(sum(dataEstWS$estState=="CA", na.rm=TRUE)>0) {
    CA_NWI <- readOGR(dsn=".\\Data\\NWI_CA", layer="NWI_CA", verbose=FALSE)
    CA_NWI <- spChFIDs(CA_NWI, paste("CA", sapply(slot(CA_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

if((sum(dataEstWS$estState=="OR", na.rm=TRUE)>0)|(sum(dataEstWS$estState=="WA/OR"))) {
OR_NWI <- readOGR(dsn=".\\Data\\NWI_OR", layer="NWI_OR", verbose=FALSE)
OR_NWI <- spChFIDs(OR_NWI, paste("OR", sapply(slot(OR_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

if(sum(dataEstWS$estState=="WA", na.rm=TRUE)>0) {
    WA_NWI <- readOGR(dsn=".\\Data\\NWI_WA", layer="NWI_WA", verbose=FALSE)
    WA_NWI <- spChFIDs(WA_NWI, paste("WA", sapply(slot(WA_NWI, "polygons"), function(x) slot(x, "ID")), sep="_"))
}

NWIfiles <- ls()[grep("_NWI", ls())]
if(length(NWIfiles)==0) warning("No data is located in NWI estuaries")
if(length(NWIfiles)==1) NWI <- get(NWIfiles[1])
if(length(NWIfiles)==2) NWI <- rbind(get(NWIfiles[1]), get(NWIfiles[2]))
if(length(NWIfiles)==3) NWI <- rbind(get(NWIfiles[1]), get(NWIfiles[2]), get(NWIfiles[3]))

#####################################################
### Obtaining climate data
#####################################################
climateData <- climateData[(climateData$wsID %in% dataEstWS$wsID), ]
climateData$wsID <- as.factor(as.character(climateData$wsID))

#PPT data
wsPPT <- aggregate(wsPPT ~ wsID, data=climateData, FUN=sum, na.rm=TRUE)
names(wsPPT)[2] <- "wsPPT_cumulative"
dataEstWS <- merge(dataEstWS, wsPPT, by="wsID", all.x=TRUE)
dataEstWS$AreaNormFreshFlow <- round((dataEstWS$wsArea * dataEstWS$wsPPT_cumulative * 0.001)/dataEstWS$estArea, 1) #multiplied by 0.001 to convert cm to m

#wsTemp data
wsTemp <- aggregate(wsTair_avg ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
wsTemp$wsTair_avg <- round(wsTemp$wsTair_avg, 2)
names(wsTemp)[2] <- "wsTair_YearMean"
dataEstWS <- merge(dataEstWS, wsTemp, by="wsID", all.x=TRUE)

#estTemp data
estTemp <- aggregate(estTair_avg ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
estTemp$estTair_avg <- round(estTemp$estTair_avg, 2)
names(estTemp)[2] <- "estTair_YearMean"
dataEstWS <- merge(dataEstWS, estTemp, by="wsID", all.x=TRUE)

#SSTemp data
SSTemp <- aggregate(estSST ~ wsID, data=climateData, FUN=mean, na.rm=TRUE)
SSTemp$estSST <- round(SSTemp$estSST, 2)
names(SSTemp)[2] <- "estSST_YearMean"
dataEstWS <- merge(dataEstWS, SSTemp, by="wsID", all.x=TRUE)


if(climate){
climateData <- merge(dataEstWS[,1:2], climateData, by="wsID", all.y=TRUE)
write.csv(climateData, paste(saveDir, paste(FileName, "_ClimateMonthly", ".csv", sep=""),  sep="\\"), row.names=FALSE)
}


# save csv files
write.csv(dataEstWS, paste(saveDir, paste(FileName, "_Estuary", ".csv", sep=""),  sep="\\"), row.names=FALSE)

# function object
watershed <- watershed[watershed$wsID %in% dataEstWS$wsID,]
NWI <- NWI[NWI$wsID %in% dataEstWS$wsID,]
NWI <- NWI[!(is.na(NWI$wsID)),]
list(data=estuaryData, ests=dataEstWS, climate=climateData, watershedPoly = watershed, NWIPoly = NWI)
}

##testing function
#setwd("C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan302012")
#x <- c("Beaver Creek Estuary")  #testing
#trial <- names2est(x, climate=TRUE, FileName="Test",
#saveDir="C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan302012\\SampleOutput")

#summary(trial)



# ---------------------------------------------------------------------
## Function 2: Creates kml maps of data

points2kml <- function(x, kmlWS=TRUE, kmlNWI=TRUE, IDvar="", saveDir, FileName="NewData"){
#x <- trial #output from points2est function
#kmlWS <- TRUE #whether to create kml file of watershed
#kmlNWI <- TRUE #whether to create kml file of NWI
#saveDir <- "C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan302012\\SampleOutput" #where to save the kml file output
#FileName <- "Test" #Preface title for kml output
#IDvar="estName" # Name of ID variable from orginal dataset.  If used, this ID variable will be printed next to sample points in KML file
print("Wait for it: preparing your maps")
print("This may take a while. Be patient.")
save2 <- paste(saveDir, "\\", FileName, sep="")

data <- x[["data"]]
ests <- x[["ests"]]
if(kmlWS) watershed <- x[["watershedPoly"]]
if(kmlNWI) NWI <- x[["NWIPoly"]]


## Point data kml file
filename <- file(paste(save2, "_SamplePoints.kml", sep=""), "w")
#write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>", filename)
write("<kml xmlns=\"http://www.opengis.net/kml/2.2\">", filename, append = TRUE)
write("<Document>", filename, append = TRUE)
write(paste("<name>", "Sample Points", "</name>", sep=" "), filename, append = TRUE)
write("<Style id=\"pnt1\">", filename, append = TRUE)
     write("<LabelStyle>", filename, append = TRUE)
       write("<scale>0.7</scale>", filename, append = TRUE)
     write("</LabelStyle>", filename, append = TRUE)
     write("<IconStyle>", filename, append = TRUE)
       write("<color>ff0000ff</color>", filename, append = TRUE)
       write("<scale>1</scale>", filename, append = TRUE)
       write("<Icon>", filename, append = TRUE)
       write(" <href>http://maps.google.com/mapfiles/kml/shapes/donut.png</href>", filename, append = TRUE)
       write("</Icon>", filename, append = TRUE)
     write("</IconStyle>", filename, append = TRUE)
   write("</Style>", filename, append = TRUE)

for(i in 1:dim(data@data)[1]){
                   write("<Placemark>", filename, append = TRUE)
                   write(paste("<name>", ifelse(is.na(IDvar), (1:dim(data@data)[1])[i], as.character(data@data[i,which(names(data@data)==IDvar)])),"</name>", sep=""), filename, append = TRUE)
                   write("<styleUrl>#pnt1</styleUrl>", filename, append = TRUE)
                   write("<Point>", filename, append = TRUE)
                   write(paste("<coordinates>", data@coords[i,1], ",", data@coords[i,2], "</coordinates>"), filename, append = TRUE)
                   write("</Point>", filename, append = TRUE)
                   write("<ExtendedData>", filename, append = TRUE)
                   for(j in 1:length(names(data@data))){
                       write(paste("<Data name=", "\"", names(data@data)[j],"\"", ">", sep=""), filename, append=TRUE)
                       write(paste("<value>", data@data[i,j], "</value>", sep=""), filename, append=TRUE)
                       write("</Data>", filename, append=TRUE)
                   }
                   write("</ExtendedData>", filename, append = TRUE)
                   write("</Placemark>", filename, append = TRUE)
               }
write("</Document>", filename, append = TRUE)
write("</kml>", filename, append = TRUE)
close(filename)



## watershed data kml file
if(kmlWS){
filename <- file(paste(save2, "_Watersheds.kml", sep=""), "w")

out <- sapply(slot(watershed, "polygons"), function(x) { kmlPolygon(x,
    lwd=2, border='darkorange',
    name="wsID",
    description = as(watershed, "data.frame")[slot(x, "ID"), "wsID"],
    #name = as(watershed, "data.frame")[slot(x, "ID"), "wsID"]
                                                                       ) })
cat(kmlPolygon(kmlname="Watersheds")$header,
    file=filename, sep="\n")
cat(unlist(out["style",]), file=filename, sep="\n")
cat(unlist(out["content",]), file=filename, sep="\n")
cat(kmlPolygon()$footer, file=filename, sep="\n")
close(filename)
}
#############################################################
## NWI data kml file
## Colors and code based on code by Tom Philippi
## More great information here: http://science.nature.nps.gov/im/monitor/stats/Rcourse/index.cfm

if(kmlNWI){
filename <- paste(save2, "_Estuaries.kml", sep="")

MapColors <- read.csv(".\\RTools\\MapColors.csv",header=TRUE)
MapColors <- MapColors[1:length(unique(NWI@data$NWIcode)),]
ColMap <- with(MapColors,rgb(R,G,B,maxColorValue=255))

# transparency
ColMapAlpha <- paste(ColMap,"95",sep="")

#ColMapAlpha <- rep(paste(ColMap,"95",sep=""), length(unique(NWI@data$NWIcode)))
names(ColMapAlpha) <- unique(NWI@data$NWIcode)


kmlPolygons <- function(PGdv,kmlFileName,Pname,Pdescription,Kname,Kdescription,ColorTheme,ColorMap) {

    kmlFile <- file(kmlFileName, "w")
    Xout <- sapply(slot(get(PGdv), "polygons"), function(x) {
                             kmlPolygon(x,
                                                #name=as(get(PGdv), "data.frame")[slot(x, "ID"), Pname], #more complex names from data
                                                name="National Wetlands Inventory",
                                                col=get(ColorMap)[as.character(as(get(PGdv), "data.frame")[slot(x, "ID"), ColorTheme])],
                                                lwd=1.5, border='white', visibility=TRUE,
                                               description=as(get(PGdv), "data.frame")[slot(x, "ID"), Pdescription] )                                                               
                         })

#cat(kmlPolygon(kmlname=Kname, #use these two lines if a header is desired
#       kmldescription=Kdescription)$header,
cat(kmlPolygon(kmlname=Kname)$header,
       file=kmlFile, sep="\n")
cat(unlist(Xout["style",]), file=kmlFile, sep="\n")
cat(unlist(Xout["content",]), file=kmlFile, sep="\n")
cat(kmlPolygon()$footer, file=kmlFile, sep="\n")
close(kmlFile)
}

# Call function:
kmlPolygons(
                      PGdv="NWI",
                      kmlFileName=filename,
                     #Pname="NWIcode",
                     Pdescription="NWIcode",
                     Kname="NWI Codes",
                     #Kdescription="June 2011 Draft CABR Vegetation Map",
                     ColorTheme="NWIcode",
                     ColorMap="ColMapAlpha"
                     )
}
}
#points2kml(trial, kmlWS=TRUE, kmlNWI=TRUE, IDvar="ID", saveDir="C:\\Documents and Settings\\mfrazi03\\Desktop\\EstuReJan132012", FileName="Test1")
