HW1_Prob2a

#### Anna Starodubtseva
#### CVEN 6833, HW 01
#### Problem 2(a) - Display the mean extreme seasonal precipitation along with the topography as a spatial map

# CLEAR MEMORY
rm(list=ls())
# Prevent Warnings
options(warn=-1)

#### Source Libraries and set working directory
mainDir="/Users/annastar/OneDrive - UCB-O365/CVEN 6833 - Advanced Data Analysis Techniques/Homework/HW 01/R"
setwd(mainDir)
suppressPackageStartupMessages(source("hw1_library.R"))
source("hw1_library.R")

## Read location data
loc=read.table(paste(mainDir, "/data/Precipitaton_data.txt", sep = ""), header = T)
prec = data.frame(lon=loc[,1],lat=loc[,2])

## Read precipitation data
# Summer
summ_prec = read.table(paste(mainDir, "/data/Summer_temporal/Max_Summer_Seas_Prec.txt", sep = ""), header = T)
# Omit extreme values
non_values <- which(summ_prec[,-1] > 1000, arr.ind = T)
summ_prec[non_values] = NaN
summ_prec = na.omit(summ_prec)
# Winter
wint_prec = read.table(paste(mainDir, "/data/Winter_temporal/Max_Winter_Seas_Prec.txt",sep = ""), header = T)
# Omit extreme values
non_values <- which(wint_prec[,-1] > 1000, arr.ind = T)
wint_prec[non_values] = NaN
wint_prec = na.omit(wint_prec)

## Calculate averages and combine with locations
summ_avg = colMeans(summ_prec[,-1]) # Mean of location columns
wint_avg = colMeans(wint_prec[,-1]) # Mean of location columns
prec$Summ = summ_avg
prec$Wint = wint_avg

wint_max = 1:ncol(wint_prec)-1
for (i in 2:ncol(wint_prec)) {
    wint_max[i] = max(wint_prec[,i])
}

## Plot spatial map
elev_grid=read.delim(paste(mainDir, "/data/Elevation_grid_1deg.txt", sep = ""), header = TRUE, sep = "\t", dec = ".")
r_pal=terrain.colors(10)
xlim1=c(-115,-101.3)
ylim1=c(29.5,42.5)
df <- melt(elev_grid$Elev, na.rm = TRUE)

df$Lon <- elev_grid$Lon
df$Lat <- elev_grid$Lat 
names(df) <- c("elev","Lon", "Lat")
zlime=range(elev_grid$Elev)
zlab=seq(0,3500, by=500)
siz1=14
siz2=10
MainStates <- map_data("state")
states_full=c("arizona","colorado","new mexico","utah")
station_prec = hcl.colors(10, palette = "Zissou 1")
name="Elevation"
g1=ggplot() +
  geom_tile(data = df, aes(x = Lon, y = Lat, fill=elev))+
  geom_polygon( data=MainStates, aes(x=long, y=lat, group=group), color="gray20", fill="transparent" )+
  ggtitle("") +
  xlab("Longitude") +
  ylab("Latitude")+
  scale_fill_gradientn(name="Elevation (m)", colours = r_pal,limits=zlime,breaks=zlab,labels=zlab,guide="colourbar")+
  theme_bw()+
  theme(legend.position="right")+
  theme(panel.grid.major = element_line(colour = "grey50", linetype = "dashed"),
        legend.text = element_text(size = 10),plot.margin=unit(c(0.0,0.01,0.01,0.01),"in"))+
  coord_sf(xlim = xlim1, ylim = ylim1,label_axes = list(bottom = "E", left = "N"), expand =FALSE)
  g2 = g1 + ggtitle("Summer Extreme Precipitation Average by Station") + geom_point(data=prec, aes(x=lon, y=lat, colour = Summ), size = 3) + scale_colour_gradientn(name = "Precipitation (mm)", colours = station_prec)
  g3 = g1 + ggtitle("Winter Extreme Precipitation Average by Station") + geom_point(data=prec, aes(x=lon, y=lat, colour = Wint), size = 3) + scale_colour_gradientn(name = "Precipitation (mm)", colours = station_prec)

show(g2)

show(g3)