# Sample script to get agrimet data into Splus
# and do some silly stuff with it.
# Version 1.0 by Ashton, May 2000

agri_read.table("C:\\WINNT\\Profiles\\guest\\Desktop\\agxl.txt", header=T)
module(spatial)

##############
# First let's analyze this data as a point pattern.
##############

ag.pt_spp(agri)   # x and y are columns in agri

### The next lines are from Joel's assignment 1 source, slightly modified

# par(pty="s")    # produces square plotting region

# produce dot map of Agrimet station data

rx_range(ag.pt$x)
ry_range(ag.pt$y)

plot(ag.pt,xlim=rx,ylim=ry)
title(main="Agrimet Station Locations")

# calculate and plot kernel density estimates

k100_kern2d(ag.pt[,1],ag.pt[,2],bw=100,nx=50,ny=50,lims=c(rx,ry),
        kernfun=kernquart)
plot(ag.pt,xlim=rx,ylim=ry)
image(k100,add=T)
points(ag.pt)
title(main="Kernel Intensity Estimates of Agrimet Station Locations\nBW = 100")

# ... and you can go on from here
# for example... change the bandwidth to 200, 300, etc
# run fhat and ghat
# go for khat, simulations, etc

##################
# Now lets work with the elevation data ala assignment 3
##################

# image and perspective plots of coal ash (from hw3, problem 1)

ag.itp_interp(agri$x,agri$y,agri$elev)
image(ag.itp,xlab="X",ylab="Y")
points(agri$x, agri$y)
title(main="Agrimet region elevations")
image.legend(ag.itp,x= 650, y=400, size=c(1,2), hor=F)

persp(ag.itp)
title(main="Agrimet region elevations")

###########
# Lets fit a planar trend surface ala hwk 3
##########

ag.lm_lm(elev ~ x + y, data=agri)  # fit the model
print(summary(ag.lm))              # look at the summary stats!!

ag.lm.itp_interp(agri$x,agri$y,fitted(ag.lm))
contour(ag.lm.itp,nint=10)    # plot a contour map of the trend
points(agri$x, agri$y)
title("Linear Trend Model for Agrimet Region Elevations")
persp(ag.lm.itp)
title("Linear Trend Model for Agrimet Region Elevations")

rlm_resid(ag.lm)
ag.lmr.itp_interp(agri$x,agri$y,rlm)
contour(ag.lmr.itp,nint=10)   # plot a contour map of the residuals

########
# Variograms of residuals!
##########

v.rlm_variogram(rlm ~ loc(x,y), data=agri)
plot(v.rlm)
title("Omnidirectional Variogram Plot")

v4.rlm_variogram(rlm ~ loc(x,y), data=agri, azimuth=c(0,45,90,135), tol.azimuth=22.5)
plot(v4.rlm)
title("Directional Variogram Plots")


########
# So what if it's anisotropic; Fit a variogram model anyway!
#######

# remember not to run model.variogram embedded in your script!
# model.variogram(v.rlm, gauss.vgram, range=200, sill=1300000, nugget=100000)

# Here's my model estimates
range_230
sill_1300000
nugget_100000
your.cov_gauss.cov

ag.coef_c(range,sill,nugget)
names(ag.coef)_c("range","sill","nugget")


#########
# Krig that bad boy using glm!
##########

ag.kg_krige(elev ~ loc(x,y) + x + y, data=agri, covfun=your.cov, 
	range=ag.coef["range"], sill=ag.coef["sill"], nugget=ag.coef["nugget"])

# get the coeffients for the glm model to compare to OLS
xm_(ag.kg$xlim[2]+ag.kg$xlim[1])/2
xd_(ag.kg$xlim[2]-ag.kg$xlim[1])/2
ym_(ag.kg$ylim[2]+ag.kg$ylim[1])/2
yd_(ag.kg$ylim[2]-ag.kg$ylim[1])/2
co.kg_c(ag.kg$coef[1]-xm*ag.kg$coef[2]/xd-ym*ag.kg$coef[3]/yd,
        ag.kg$coef[2]/xd,ag.kg$coef[3]/yd)
print(co.kg)
print(ag.lm$coef)

# predict values for the whole area

ag.pk_predict.krige(ag.kg, grid=list(x=c(rx,40), y=c(ry,40)))
x_unique(ag.pk$x)
y_unique(ag.pk$y)
zh_matrix(ag.pk$fit,40,40)
zse_matrix(ag.pk$se,40,40)
par(pty="s")
contour(x,y,zh,nint=20)
points(ag.pt)
title(main="Kriging Interpolation of Agrimet Elevation")

persp(x,y,zh)
title(main="Agrimet Elevations from GSM/Kriging")

contour(x,y,zse,nint=10)
points(ag.pt)
title(main="Kriging Standard Errors")
persp(x,y,zse)
title(main="Standard Error Surface")