########################################################################
#                                                                      #
# MATH574 - Geostatistics                                              #
# Prof. Dr. Donald Myers                                               #
#                                                                      #
# Term-project: Anaylizing climatological data from Sao Paulo State,   #
# Brazil:                                                              #
#    - Annual precipitation (mm); and                                  #
#    - Annual mean temperature (oC).                                   #
#                                                                      #
# * elevation data is also provided                                    #
#                                                                      #
# Author: Rafael Rosolem (rafael@hwr.arizona.edu)                      #
# Department: Hydrology and Water Resources                            #
# The University of Arizona                                            #
#                                                                      #
# Date: 02.14.2007                                                     #
# Last update: 02.19.2007                                              #
#                                                                      #
########################################################################

# Opening libraries                                                       
########################################################################
library(geoR) # geoR package
library(gstat) # gstat package
library(scatterplot3d) # scatterplot3d package

# Setting the path and file name                                       
########################################################################
path <- "/home/rafael/math574/term-project/"
input <- "dataset_rafael.dat"

# Opening data file                                       
########################################################################
data <- read.table(paste(path,input,sep=""),header = TRUE)

# Converting data object to geodata class
########################################################################
data_geoR <- as.geodata(data,coords.col = c(1,2),data.col = c(3,4,5))

# Basics statistics
########################################################################
data_summary <- summary(data_geoR) # Data summary
data_stdev <- sd(data_geoR$data) # Standard deviation
data_cor <- cor(data_geoR$data) # Data correlation

# Histograms
########################################################################
X11(display = "",width = 9,height = 4,pointsize = 12,bg = "transparent")
#png(file=paste(path,"data_hist.png",sep=""),height=576,width=576,
#    res=600,bg="transparent")
   par(mfrow = c(1,3)) # 1x3 graph
   temp_hist <-  hist(data_geoR$data[,"temp"],breaks = c(17:26),
                      col = "red",border = "black",
		      main = NULL,
		     xlim = c(min(data_geoR$data[,"temp"])-1,
		     max(data_geoR$data[,"temp"])+1),ylim = NULL,
		     xlab = "Temperature (oC)",ylab = "Freq") # Temperature
   prec_hist <-  hist(data_geoR$data[,"prec"],breaks = c(250*0:14),
		      col = "darkblue",border = "black",
		     main = NULL,
		     xlim = c(min(data_geoR$data[,"prec"])-500,
		     max(data_geoR$data[,"prec"])+500),ylim = NULL,
		     xlab = "Precipitation (mm)",ylab = "Freq") # Precipitation
   elev_hist <-  hist(data_geoR$data[,"elev"],breaks = c(100*0:12),
		      col = "gray",border = "black",
		     main = NULL,
		     xlim = c(min(data_geoR$data[,"elev"]),
		     max(data_geoR$data[,"elev"])+200),ylim = NULL,
		     xlab = "Elevation (m)",ylab = "Freq") # Elevation
#dev.off()

# Scatterplots
########################################################################
#X11(display = "",width = 9,height = 9,pointsize = 12,bg = "transparent")
#png(file=paste(path,"data_pairs.png",sep=""),height=576,width=576,
#    res=300,bg="transparent")
   data_pairs <- pairs(data_geoR,
                       labels=c("Eastings (km)","Northings (km)","Temp(oC)","Prec(mm)","Elev(m)"))
#dev.off()

# 3d plots
########################################################################
X11(display = "",width = 12,height = 11,pointsize = 12,bg = "transparent")
#png(file=paste(path,"data_3d.png",sep=""),height=576,width=576,
#    res=600,bg="transparent")
   par(mfrow = c(2,2)) # 2x2 graph
   temp_3d <- scatterplot3d(data_geoR$coords[,"X"],
                            data_geoR$coords[,"Y"],
			    data_geoR$data[,"temp"],
			    xlab="Eastings (km)",ylab="Northings (km)",
			    zlab="Temp(oC)")
   temp_3d.lm <- lm(data_geoR$data[,"temp"] ~ 
                    data_geoR$coords[,"X"] + data_geoR$coords[,"Y"])
   temp_3d$plane3d(temp_3d.lm)

   prec_3d <- scatterplot3d(data_geoR$coords[,"X"],
                            data_geoR$coords[,"Y"],
			    data_geoR$data[,"prec"],
			    xlab="Eastings (km)",ylab="Northings (km)",
			    zlab="Prec(mm)")
   prec_3d.lm <- lm(data_geoR$data[,"prec"] ~
                    data_geoR$coords[,"X"] + data_geoR$coords[,"Y"])
   prec_3d$plane3d(prec_3d.lm)

   elev_3d <- scatterplot3d(data_geoR$coords[,"X"],
                            data_geoR$coords[,"Y"],
			    data_geoR$data[,"elev"],
			    xlab="Eastings (km)",ylab="Northings (km)",
			    zlab="Elev(m)")
   elev_3d.lm <- lm(data_geoR$data[,"elev"] ~
                    data_geoR$coords[,"X"] + data_geoR$coords[,"Y"])
   elev_3d$plane3d(elev_3d.lm)

   all_3d <- scatterplot3d(data_geoR$data,
			   xlab="Temp(oC)",ylab="Prec(mm)",
			   zlab="Elev(m)")
   all_3d.lm <- lm(data_geoR$data[,"elev"] ~
                   data_geoR$data[,"temp"] + data_geoR$data[,"prec"])
   all_3d$plane3d(all_3d.lm)

#dev.off()
