File:       a13_satelliteImage.txt
By:         Dan Carr
Date:       copyright 2005, 2005, 2006, Modified for R April 2007, 2008
              

Purpose:    Illustrate different paired comparisons
            1) pairing values by common location
            2) pairing quantiles by common cumulative probabilities	         
            
            Working with image data
            A little string processing
 
Files:      Obtain files from the web site
               ndviSept82.txt
               ndviSept92.txt 

Due:        Run script as is.
            Find the lines with ### and ? marks
            Make reasonable changes and rerun
            Comment out (#) the previous lines for the same purpose
            Your pdf plot from 1,
            2 plots from 2,
            Your plot from 3.  

Notes:      Study the script enough to understand the
	       general steps and know where to look if
            you need similar code later on.		

Warning     Some plots created below are very similar to
            plots in a copyrighted encyclopedia article.  


Warning:    If you have a pdf plot open with Acrobat and 
            try to produce a new pdf file with the same name
            using R, it won't work.  Close Acrobat first
            you want to rewrite the file 


1. Images plots and hexagon binned mean and difference of paired quantiles

NDVI stands for normalized difference vegetation index.
NDVI  = (NIR - RED)/(NIR+RED)
NIR   = near infrared spectral reflectance measurements 
RED   = infrared spectral reflectance measures

Living green plants cannot make good use of solar radiation in the 
near infrared part of the spectrum. The leaf surfaces have evolved
to scatter radiation in this part of the spectrum and avoid plant
overheating. 

Snow and clouds scatter radiation in both the NIR and RED parts
of the spectrum. Large differences, NIR - Red, are indicative
of living green plants. 

Normalization by the sum produces a unitless index.         


1.1 Read NDVI datafor processing  NDVI data, fix data and labels, and re-arranges

Read files
Replace ocean values of -9.9990 with -0.06
Edit the longitude labels
Rearrange matrix for use in the image() function
    Transpose and put y values in the reverse order

## Run

# Read file Sept 1982___________________________________

ndvi82= read.table("ndviSept82.txt", header= TRUE,row.names=1)
ndvi82.mat = as.matrix(ndvi82)

# Fix column labels

lab82 = dimnames(ndvi82.mat)[[2]]    # column labels
#   [1] "X.120.0" "X.119.0" "X.118.0" "X.117.0" "X.116.0" "X.115.0"
# [261] "X140.0"  "X141.0"  "X142.0"  "X143.0"  "X144.0"  "X145.0"

# Remove X
# Remove trailing .0
# R has substituted . for -, put it back  

nch = nchar(lab82)
lab82 = substring(lab82,2,nch-2)    #remove leading X
                                    #and trailing ".0"
nch =nchar(lab82)
tst = substring(lab82,1,1)
pre = ifelse(tst=='.','-','')  # find -
base = ifelse(tst=='.',substring(lab82,2,nch),substring(lab82,1,nch))
nlab82 = paste(pre,base,sep='') # add prefix to base
dimnames(ndvi82.mat) = list(dimnames(ndvi82.mat)[[1]],nlab82)

# Transpose the matrix and reverse the y axis for 
# the R image function
s82 = t(ndvi82.mat)      
s82 = s82[,ncol(s82):1]
mean82 = mean(s82[s82 > -0.06])  

# Sept 1992___________________

ndvi92 = read.table("ndviSept92.txt", header= TRUE,row.names=1)
ndvi92.mat = as.matrix(ndvi92)

lab92 = dimnames(ndvi92.mat)[[2]]
nch = nchar(lab92)
lab92 = substring(lab92,2,nch-2)               # remove leading X
                                               # remove trailing .0
nch =nchar(lab92)
tst = substring(lab92,1,1)
pre = ifelse(tst=='.','-','')  # find -
base = ifelse(tst=='.',substring(lab82,2,nch),substring(lab82,1,nch))
nlab92 = paste(pre,base,sep='') # add prefix to base
dimnames(ndvi92.mat) = list(dimnames(ndvi92.mat)[[1]],nlab92)

# Transpose the matrix and reverse the y axis for 
# the R image() function
s92 = t(ndvi92.mat)
s92 = s92[,ncol(s92):1]
mean92 = mean(s92[s92 > -0.06 ])

## End

1.2 Define the break points and color classes+++++++++++++++++++++++++
  
In image(), break points are used with cut()
to define a subscript for each image pixel
that indexes an element in a vector of colors.

##Run
breaks = c(-.061,-.059,-.000001,.1,.2,.3,.4,.5,.6,.7,.8)
# ndvi = -.06 is ocean and maps into subscript 1
# -.059 < ndvi <= -.000001 maps into subscript 2
# -.000001 < ndvi <= .1    maps into subscript 3 
# ...
#  .7 < ndvi < .8          maps into subscript 10               
#  The large value in the files is < .8
## End

Define your own colors corresponding to the subscripts

## Run to define 10 colors to match with the subscripts

cmat = matrix(c(
 .00, .00,.87,   # blue for ocean
 .00, .75,.75,
1.00, .95,.70,
 .84, .74,.41,
 .69, .48,.21,
 .53, .22,.00,
 .83,1.00,.51,
 .58, .83,.36,
 .31, .67,.20,
 .04, .40,.05),  # dark green
 ncol=3,byrow=T)
mycolors = rgb(cmat[,1],cmat[,2],cmat[,3],max=1)

### cmat = ?  Made a different color table

## End

1.3 Device call, page layout, graphics parameters  

pdf(file='a13_imagemd.pdf',width=7.5,height=9.2)

pan = panelLayout(nrow=3,ncol=2,
       rowSep=c(0,.15,.5,.07),
       colSize= c(20,1),
       topMar=.5,leftMar=.5,bottomMar=.5)

pan2 =panelLayout(nrow=3,ncol=3,
       rowSep=c(0,.15,.5,.07),
       colSize= c(3,14,4),    #note sum is 21 as above
       topMar=.5,leftMar=.5,bottomMar=.5)

# Graphics parameters 

ygrid = c(-50,-25,0,25,50)
labs = dimnames(s82)
x = as.numeric(labs[[1]])
y = as.numeric(labs[[2]])
rx = c(-120.5,150.5)
ry = c(-50.5,50.5)

#_________________________Panel 1______________________

panelSelect(pan,1,1)
panelScale(rx=rx,ry=ry)

# Maps often need to have a specific aspect ratio
# Here I wanted to represent degrees equally 
# 	both horizontally and vertically
# One method is to compare the panel
#       range(longitude) and range(latitude) ratio
#       to the ratio of x and y inches used in panel 
# In Splus this is a par()$uin vector that shows
#       x axis user units per inch and
#       y axis user units per inch
#       

tmp = par()$pin
tmp[1]/tmp[2]
diff(range(x))/diff(range(y))

image(x,y,z=s92,breaks=breaks,col=mycolors,
      xlab="",ylab="",xaxs='i',yaxs='i',axes=F,add=T)
axis(side=2,at=ygrid,labels=as.character(ygrid),mgp=c(2,.5,0),
	cex.axis=.8,las=2,tck=-.03)
mtext('Sept. 1992',side=1,line=-1.2,cex=1,col="white")
text(rx[1]-.1*diff(rx),mean(ry),'Sept. 1992',cex=.8)
panelOutline()

#__________________________Legend____________________

panelSelect(pan,1,2)
panelScale(rx=c(0,1),ry=c(-.1,.8))
px = rep(c(.4,.4,1,1,NA),9)
yval = seq(-.1,.8,length=10)
py = yval[ c(1,2,2,1,NA)+ rep(seq(0,8),rep(5,9))]
polygon(px,py,col=mycolors[2:10])
polygon(px,py,density=0)
par(xpd=T)
text(1.15,yval,format(yval),cex=.8,adj=0)
mtext(' NDVI',side=3,line=.3,cex=.8)

#_________________________Panel 2_____________________

panelSelect(pan,2,1)
panelScale(rx=rx,ry=ry)
image(x,y,z=s92,breaks=breaks,col=mycolors,
      xlab="",ylab="",xaxs='i', yaxs='i',add=T)
axis(side=2,at=ygrid,labels=F,mgp=c(2,.5,0),
	cex.axis=.8,las=2,tck=-.03)
axis(side=2,at=ygrid[-5],labels=,as.character(ygrid[-5]),mgp=c(2,.5,0),
	cex.axis=.8,las=2,tck=-.03)

axis(side=1,mgp=c(2,.3,0),tck=-.025,cex.axis=.8)
mtext('Sept. 1982',side=1,line=-1.2,cex=1,col="white")

panelOutline()

#________________________Panel 3________________________________

panelSelect(pan2,3,2)
v82 = as.vector(s82)
v92 = as.vector(s92)
nx = (v92+v82)/2
ny = v92-v82
# Select land values
land = nx > -.05 
nx = nx[land]
ny = ny[land]

# Determine panel shape so hexagons have the right shape
tmp = par()$pin
aspect = tmp[1]/tmp[2]

# Hexagon binning
ans = hbin(nx,ny,xbins=40,aspectRatio=aspect)
###ans = hbin(nx,ny,xbins=?,aspectRatio=aspect) # Try a difference # of bins

# Smoothing with sqrt count weights
sm = loess(ans$yMean~ans$xMean,weights=sqrt(ans$count),span=.3)
###sm = loess(ans$yMean~ans$xMean,weights=sqrt(ans$count),span=?) # different span

# Calculate limits for hexagon binned plot
tmp = attributes(ans)    #  access additional object values  
nrx = mean(tmp$xlim)+ 1.06*diff(tmp$xlim)*c(-.5,.5)
nry = mean(tmp$ylim)+ 1.06*diff(tmp$ylim)*c(-.5,.5)

panelScale(rx=nrx,ry=nry)
temp.ans = ans
hbin.draw(ans,style='lattice',minarea=1,maxarea=1,
          col.one=mycolors[2]) #  bin background 
hbin.draw(ans,style='lattice',col.one="black")
panelOutline()
hgrid = c(-.4,-.2,0,.2,.4)

axis(side=1,mgp=c(2,.3,0),cex.axis=.8,tck=-.03)
axis(side=2,at=hgrid,labels=format(c(-.4,-.2,0.0,.2,.4)),
     mgp=c(2,.5,0),las=2,cex.axis=.8,tck=-.025)

# abline(y=0,lwd=2)
tx = unique(sort(ans$xMean))
ty = predict(sm,tx)
lines(range(tx),c(0,0),col='black',lwd=1)
lines(tx,ty,col="red",lwd=3)

mtext('Difference: 1992 - 1982 NDVI',side=2,line=2,cex=.8)
mtext('Grid Cell Average of 1992 and 1982 NDVI',side=1,line=1.4,cex=.8)

# Add label at the top
panelSelect(pan,1,1)
panelScale()
par(xpd=T)
text(.5,1.2,'Binned Grid Cell Based Comparison with Smooth',adj=.5)
text(.5,1.08,'For Two NDVI Images',adj=.5)

dev.off()

2. Construct a qqplot comparing distributions==========================

Find pixels that are not water in both images.
Sort to obtain quantiles from each images
Use points to obtain the common probability points.
Put in labels in the plot command

2.1  QQplot

## Run
v82 = as.vector(s82)
v92 = as.vector(s92)
good = v82 > -.06 & v92 > -.06

q82 = sort(v82[good])
q92 = sort(v92[good])

p82 = ppoints(q82)
p92 = ppoints(q92)

pcommon = ppoints(500)  # restrict qqplot to 500 points
###pcommon = ppoints(?)

q82match = approx(p82,q82,pcommon)$y
q92match = approx(p92,q92,pcommon)$y

# pdf(file='a9_imageqq.pdf',width=7.5,height=8.7)
windows()
plot(q82match,q92match,type='l',pch=19,col="#00C000",
     xlab='1982 NDVI quantiles',ylab='1992 NDVI Quantiles',
     main='QQplot For Distibution Comparison',
     sub='Red Line = Equality, Gray Line Robust Fit',lwd=2)
abline(0,1,lwd=3,col="red")

# fit line to 1st and 3rd quartiles
q82ref = approx(p82,q82,c(.25,.75))$y
q92ref = approx(p92,q92,c(.25,.75))$y

ans = lm(q92ref ~ q82ref)
abline(coef=ans$coef,col="#808080",lwd=2)

# ans = rreg(q82match,q92match) not in R
# abline(ans$coef)

## End

2.2  Paired location Plot

## Run

plot(v82[good],v92[good],main="Paired by Pixel Location",
xlab="1982 NDVI",ylab="1992 NDVI",sub="Red line = Equality")
abline(0,1,col='red',lwd=2)

## End

3.  Box Plots from scratch

A close variation on this was copyrighted. 

Select a different part of the map for comparison
and change the label. 

## Run
# Select South America
sam82 = s82[39:86,1:63]
sam92 = s92[39:86,1:63]

### Select some other place
### look at dimnames(s92) to see the
### correspondence between degrees and subscripts
###sam82 = s82[?:?,?:?]
###sam92 = s92[?:?,?:?]

samlabs = dimnames(sam82)
samx = as.numeric(samlabs[[1]])
samy = as.numeric(samlabs[[2]])

# check image
image(samx,samy,sam82,breaks=breaks,col=mycolors,
      xlab="Logitude",ylab="Latitude",main="1982 NDVI")

# remove ocean
sam82 = sam82[sam82>-.059]
sam92 = sam92[sam92>-.059]

# calculate box plot statistics
boxsum = boxplot(sam82,sam92,notch=T,plot=F)
stats = boxsum$stats
conf = boxsum$conf

windows()

panels = panelLayout(nrow=1,ncol=2,
	      colSize=c(1,4.8),
              leftMar=0,rightMar=0,topMar=.5,bottomMar=.5)

#___________________________plot labels_______________________

panelSelect(panels,1,1)
panelScale()
cex = 1
tcex = 1
a = .25
y = (2:1-a)/(3-2*a)
text(c(.1,.1),y,c('Sept. 1992','Sept. 1982'),cex=cex,adj=0)

#____________________Plot NDVI boxplot _________________________

rx = range(stats)
rx = c(-.1,.8)
panelSelect(panels,1,2)
panelScale(rx=rx,ry=c(0,1))
panelFill(col='#A0A0A0')
panelGrid(x=c(.2,.5),col="white")

# construct rectangles for boxplots
a = .25
seps = diff((1:2-a)/(3-2*a))/2
medline = .85*c(-seps,seps)
seps = c(-seps,-seps,seps,seps)
thick = .50*seps
thin  = .25*seps
ext = .125*seps
   
y = (2:1-a)/(3-2*a)
for (i in 1:2){
       vals = stats[,i]
       cv = conf[,i]
       polygon(vals[c(1,5,5,1)],y[i]+thin,col="blue",border=F)
       polygon(vals[c(1,5,5,1)],y[i]+thin,col="black",density=0,border=T)
       polygon(vals[c(2,4,4,2)],y[i]+thick,col="#00B0F0",border=F)
       polygon(vals[c(2,4,4,2)],y[i]+thick,col=1,density=0,border=T)
       polygon(c(cv,rev(cv)),y[i]+ext,col="white",border=F)
       lines(vals[c(3,3)],y[i]+medline,col=1,lwd=1)
   }
axis(side=1,at=c(-.1,.2,.5,.8),labels=c('-.1','.2','.5','.8'),
     mgp=c(2,.3,0),cex=tcex)

mtext('NDVI from South America',side=1,line=1.4,cex=cex)
### mtext('NDVI from ??',side=1,line=1.4,cex=cex)

panelOutline()

##End