## Do this after the code in oak1.s
## x1 is a 14 x 10 matrix of poison oak sighting totals


## 1) Smooth the image, onto an n1 x n1 grid, ranging from [0,50] x [0,50], 
## using bandwidth h
gridlon = rep(seq(0,50,length=14),times=rep(10,14))
gridlat = rep(seq(0,50,length=10),14)
x2 = c(t(x1))
## Now gridlon, gridlat, x2 correspond to the data, pixel by pixel.
h = 2
n1 = 7
mygrid1 = seq(1,50,length=n1)
mygrid2 = seq(1,50,length=n1)
a1 = matrix(0,ncol=n1,nrow=n1)
for(i in 1:n1){
  for(j in 1:n1){
    a1[i,j] = sum(x2 * dnorm(
	    sqrt((mygrid1[i]-gridlon)^2 + (mygrid2[j]-gridlat)^2),sd=h)) / 
    sum(dnorm(sqrt((mygrid1[i]-gridlon)^2 + (mygrid2[j]-gridlat)^2),sd=h))
    }
  }

par(mfrow=c(1,2))
image(mygrid1,mygrid2,a1,xlab="x",ylab="y",zlim=range(x1),
    col=grey(c(64:20)/64))
mtext(s=3,paste("smoothed totals, bw = ",as.character(h)))

## legend
x = a1
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,
    zlim=range(x),col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"sightings",srt=-90)

## try different values of h above

## 2) Fit the spatial trend to x1 by regression
a2 = lm(x2 ~ gridlon + gridlat)
a3 = a2$coef
## plot the results on a grid of the same dimension as the data
a4 = matrix(0,ncol=10,nrow=14)
mygrid1 = seq(1,50,length=14)
mygrid2 = seq(1,50,length=10)
for(i in 1:14){
    for(j in 1:10){
      a4[i,j] = a3[1] + a3[2]*mygrid1[i] + a3[3]*mygrid2[j]
      }
    }
image(mygrid1,mygrid2,a4,xlab="x",ylab="y",
    zlim=range(a4),col=grey(c(64:20)/64))
mtext(s=3, "spatial trend fit by regression")

## legend
x = a4
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,
    zlim=range(x),col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"sightings",srt=-90)

## 3) Show residuals from regression

a5 = a2$resid
a6 = matrix(a5,ncol=10,byrow=T)
image(mygrid1,mygrid2,a6,xlab="x",ylab="y",zlim=range(a5),
    col=grey(c(64:20)/64))
mtext(s=3,"Residuals from linear regression")
######### LEGEND:
x = a6
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"oak sightings",srt=-90)
sqrt(mean(x2^2)) ## RMS of the data
## compare with 
sqrt(mean(a6^2)) ## RMS of the regression residuals

## 4) Median Polish
med.polish = function(data,k){
## Kevin Nichols' code
n1 = length(data[,1])
n2 = length(data[1,])
col1 = apply(data,1,median)
data1 = data - col1
row1 = apply(data1,2,median)
data2 = data1 - matrix(rep(row1,n1),ncol=n2,byrow=T)
overall2 = median(col1) 
col2 = col1 - median(col1)
col3 = col2
row3 = row1
for(i in 2:k){
col1 = apply(data2,1,median)
data1 = data2 - col1
row1 = row3 - median(row3)
overall1 = overall2 + median(row3)
col2 = col3 + col1 
row2 = apply(data1,2,median)
data2 = data1 - matrix(rep(row2,n1),ncol=n2,byrow=T)
row3 = row1 + row2
overall2 = overall1 + median(col2)
col3 = col2 - median(col2)
}
return(list(data2,col3,row3,overall2))
}

par(mfrow=c(2,2))
image(mygrid1, mygrid2,x1,zlim=range(x1),col=grey(c(64:20)/64),
    xlab="",ylab="",main="oak sightings") 
mp1 = med.polish(x1,1)
## legend
x = x1
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,
    zlim=range(x2),col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"oak sightings",srt=-90)
image(mygrid1,mygrid2,mp1[[1]],zlim=range(x1),col=gray((64:20)/64), 
    xlab="",ylab="",main="residuals after \n median polishing")
plot(c(1,max(length(mp1[[2]]),length(mp1[[3]]))),range(c(mp1[[2]],mp1[[3]])), 
    xlab="row/column number",ylab="value",type="n") 
lines(mp1[[2]],lty=2) ##  row effect
lines(mp1[[3]],lty=1) ## column effect
legend(7,.5,legend=c("row effect","column effect"),lty=c(1,2))
mp1[[4]] ## overall mean


## 5) Median polish of smoothed data on a 7 x 7 grid from part 1.
par(mfrow=c(2,2))
n1 = 7
mygrid1 = seq(1,50,length=n1)
mygrid2 = seq(1,50,length=n1)
image(mygrid1,mygrid2,a1,xlab="x",ylab="y",zlim=range(x1),
    col=grey(c(64:20)/64))
mtext(s=3,paste("smoothed sighting totals, bw = ",as.character(h)))
## legend
x = a1
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,
    zlim=range(x2),col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"oak sightings",srt=-90)
mp1 = med.polish(a1,4)
image(mygrid1,mygrid2,mp1[[1]],zlim=range(x1),col=gray((64:20)/64), 
    xlab="",ylab="",main="residuals after \n median polishing")
plot(c(1,max(length(mp1[[2]]),length(mp1[[3]]))),range(c(mp1[[2]],mp1[[3]])), 
    xlab="row/column number",ylab="value",type="n") 
lines(mp1[[2]],lty=2) ##  row effect
lines(mp1[[3]],lty=1) ## column effect
legend(1,.5,legend=c("row effect","column effect"),lty=c(1,2))
mp1[[4]] ## overall mean