17.5 Exploring for Similar Data Structures

# Fig. 17.5: Robust and Classical DA for lithology 82, 83 based on PCs with CV library(StatDA) data(chorizon) data(kola.background) lit=chorizon[,"LITO"] Region=rep(1, length(lit)) Region[lit==82] <- 3 Region[lit==83] <- 3 X=chorizon[,"XCOO"] Y=chorizon[,"YCOO"] XY=cbind(X,Y) # cut region xcent=c(800000,7535000) sel=sqrt((X-xcent[1])^2+(Y-xcent[2])^2)<70000 # true representation of x and y axis of map for plot xwid=diff(range(X))/12e4 ywid=diff(range(Y))/12e4 pdf("fig-17-5.pdf",width=2*xwid,height=1*ywid) par(mfrow=c(1,2),mar=c(1.5,1.5,1.5,1.5)) el <- c("Al_XRF","Ca_XRF","Fe_XRF","K_XRF","Mg_XRF","Mn_XRF","Na_XRF","P_XRF","Si_XRF","Ti_XRF") x=chorizon[,el] # log-ratio x1=x/x[,"Ti_XRF"] x=log10(x1[,1:9]) x=scale(x) ####################################################################################### # LDA ####################################################################################### set.seed(100) x.lda=lda(x, Region, prior=c(0.5,0.5), CV=TRUE) cl=x.lda$class pt=x.lda$post plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n") plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T) points(XY[Region==3,],pch=17,cex=0.8) points(XY[Region==1,],pch=3,cex=0.8) legend("topright",pch=c(17,3),pt.cex=c(0.8,0.8), legend=c("Alkal. Intr.", "Others"), cex=0.8,title="Original Groups") # scalebar scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8) # North arrow Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6) ####################################################################################### # Results Group 1 # generate plot with background plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n") plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T) v.red=Region points(XY[v.red==3,],pch=17,cex=0.8,col=gray(1-pt[v.red==3,2])) points(XY[v.red==1,],pch=3,cex=0.8,col=gray(1-pt[v.red==1,2])) legend("topright",pch=c(17,3),pt.cex=c(0.8,0.8), legend=c("Alkal. Intr.","Others"), cex=0.8,title="Prediction Alkal. Intr.") # scalebar scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8) # North arrow Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6) dev.off()