HW1_P2_i_July

#### Problem 2 - Homework 1 for CVEN6833 - Spatial map July
#### Author: Alvaro Ossandon

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

#data required: 
Month=7 # 1: Monthly precipitation of January; 7: Monthly precipitation of July
save=FALSE # True for saving figure as pdf file

#### SOURCE LIBRARIES (suppress package loading messages) AND SET WORKING DIRECTORY
mainDir="C:/Users/DELL/Dropbox/Boulder/Courses/Advance data analysis/HW/HW1"
setwd(mainDir)
suppressPackageStartupMessages(source("Lib_HW1.R"))
source("Lib_HW1.R")

#### create a folder to save the figures
if (Month==1){
  diraux=paste(mainDir,"/Montly_Precep_January",sep = "")  
  name_fig="Colorado Topography and Station Locations"
  name_fig1="Colorado Mean January Precipitation (1987-2007)"
}else{
  diraux=paste(mainDir,"/Montly_Precep_July",sep = "")
  name_fig="Colorado Topography and Station Locations"
  name_fig1="Colorado Mean July Precipitation (1987-2007)"
}
setwd(diraux)

# Change projection to geographical coordinates for the location of stations
proj4string <- "+proj=utm +zone=13 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
xx=read.table("http://cires1.colorado.edu/~aslater/CVEN_6833/colo_utm_precip.dat")
xy=cbind(xx[,2],xx[,1])
pj=project(xy, proj4string, inv=TRUE)
Est= data.frame(lat=pj[,2], lon=pj[,1])

#read the monthly precipitation
x = read.table(paste("http://cires1.colorado.edu/~aslater/CVEN_6833/colo_monthly_precip_0",Month,".dat",sep=""))
names(x)=c("Lat","Lon","Elevation","Precipitation")
a<- which(x[,4]>-999.999)
x2<-x[a,]

#  read the topography map
x1 = read.table("http://cires1.colorado.edu/~aslater/CVEN_6833/colo_dem.dat")
names(x1)=c("Lat","Lon","Elevation")

# interpolate the topography map and monthly precipitation in order to create a grid for each one
fld <- with(x1, interp(x =Lon , y = Lat, z = Elevation, extrap=TRUE, duplicate = "strip"))
fld2 <- with(x2, interp(x =Lon , y = Lat, z = Precipitation, extrap=TRUE, duplicate = "strip"))

# get the colorado state limits
states <- map_data("state")
west_coast <- subset(states, region %in% "colorado")

# prepare data in long format
df <- melt(fld$z, na.rm = TRUE)
names(df) <- c("x", "y", "Elevation")
df$Lon <- fld$x[df$x]
df$Lat <- fld$y[df$y]  
df2 <- melt(fld2$z, na.rm = TRUE)
names(df2) <- c("x", "y", "Precipitation")
df2$Lon <- fld2$x[df2$x]
df2$Lat <- fld2$y[df2$y]          
   

g1=ggplot()+
  geom_tile(data = df, aes(x = Lon, y = Lat, fill=Elevation))+
  ggtitle(name_fig) +
  xlab("Longitude") +
  ylab("Latitude")+
  xlim(range(x1$Lon,x2$Lon,Est$lon))+
  ylim(range(x1$Lat,x2$Lat,Est$lat))+
  scale_fill_gradientn(name=" (m)", colours = terrain.colors(10),
                       breaks = as.integer(c(seq(800,round(max(df$Elevation),digits = 0),by=200))), limits=c(800,max(df$Elevation)), oob=squish)+
  theme(legend.position="right")+
   guides(fill = guide_colorbar(barwidth =1, barheight = 20))+
  theme(panel.background = element_rect(fill = "white", colour = "grey50"),
         panel.grid.major = element_line(colour = "grey50", linetype = "dashed"),
         plot.title = element_text(size = 15, face = "bold", hjust = 0.5),
         legend.title = element_text(size = 10, vjust = 0.8),
         axis.text = element_text(size = 12),
         axis.title.x = element_text(size = 15, vjust = -0.5),
         axis.title.y = element_text(size = 15, vjust = -0.1),
         legend.text = element_text(size = 10))
    g2=g1+
  #stat_contour(data = df2, aes(x = Lon, y = Lat, z=Precipitation, colour =..level..),geom = "contour",show.legend=TRUE ,breaks =c(seq(as.integer(range(df2$Precipitation))[1],as.integer(range(df2$Precipitation))[2], by=20)))+
  #g3= direct.label(g1, "bottom.pieces")+
   geom_polygon(data = west_coast, aes(x = long, y = lat, group = group),alpha = 0, fill = "white", color = "gray30",size = 1) +
      geom_point(data=Est, aes(x=lon, y=lat), colour ="red", size=1)
    if (save==TRUE){
      ggsave(paste(name_fig,".pdf",sep = ""), width = 9, height = 6) # save figure
      dev.off()
    }
 show(g2)

if (save==TRUE){
  pdf(paste(name_fig1,".pdf",sep = ""), width = 9, height = 6) # save figure
}
surface(fld2,xlab="Longitude",ylab ="Latitude",main=name_fig1
        ,horizontal = FALSE,legend.shrink = 0.9, zlim=range(x2$Precipitation), xlim=range(x1$Lon,x2$Lon,Est$lon),ylim=range(x1$Lat,x2$Lat,Est$lat))
grid(col= "grey50", lty = "dashed")
US( add=TRUE, col="gray30", lwd=2)

if (save==TRUE){
  dev.off()
}