## Examples of some functions in "R" for exploratory 
## data analysis for spatial data

## some libraries needed for functions below
library(geoR)
library(scatterplot3d)
library(fields)
library(lattice)
library(akima)#need this for interp
#library(geoRglm)
## Note: image and contour are available in the R base package
## Note: the coordinates (x,y) need to be provided in ascending
## order for image and contour to work


## read swiss rain data
swissrain = read.table("/Users/pcaragea/teaching/Stat406/2008/Rstuff/swissrain.txt",header=TRUE,sep=" ")
attach(swissrain)#id,xccord,ycoord,rain

## interp takes above data and puts it into matrix form appropriate for image and contour plots
## image and contour plots want input=list with x=vector of x-values, y=vector of y-values and
## z=matrix of data values at each combination of (x,y) values
swissrain.interp = interp(xcoord,ycoord,rain)
image(swissrain.interp)
contour(swissrain.interp, add=T)
## 3D- PLOTTING
## dropline plot
scatterplot3d(x=xcoord,y=ycoord,z=rain,xlab="X coordinate",ylab="Y coordinate",zlab="Rain",type='h',cex.axis=1.3)

## other 3d plotting functions
cloud(ycoord ~ xcoord * rain,xlab="X coordinate",ylab="Y coordinate",zlab="Rain", type = 'p')

plot.surface(swissrain.interp, type="p") # perspective/drape plot
plot.surface(swissrain.interp, type="C")  # contours with a legend

## empirical variograms 
variog.rain=variog(coords=cbind(xcoord,ycoord),data=rain)
plot(variog)


#Alternatively, put data into "geodata" format
#(need the input in the form X,Y,Z)
rain.g=as.geodata(cbind(xcoord,ycoord,rain))
summary(rain.g)
plot(rain.g)
points(rain.g)
variog.rain.g=variog(rain.g)
plot(variog.rain.g)
variog.rain.cloud<-variog(rain.g,option="cloud",max.dist=150000)
plot(variog.rain.cloud)
#cloud with boxplot---MOM
variog.rain.boxplot<-variog(rain.g,bin.cloud=T)
plot(variog.rain.boxplot,pch=20,bin.cloud=T);title(main="Binned variog(Class) with boxplots")

##compute the robust variogram---Cressie and Hawkins
variog.rain.rob=variog(coords=cbind(xcoord,ycoord),data=rain,estimator.type="modulus",uvec=20)
plot(variog.rain.rob)