################################################################################### # # H scatterplots # # Last review: 23/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### varcolXY <- c(1,2) varlistindex <- 11 ndecimals <- 3 lag <- 75 n.hplots <- 8 # acresce o gráfico do correlograma sizePOINTS=0.75 titleGRAPH="H-scatterplots and correlogram" ################################################################################### par(mfrow=c(3,3), mar=c(3.5, 3.5, 2, 2), oma=c(0, 0, 2, 0), mgp=c(2.0, 0.8, 0)) # # Duplicate gtd -> gtdnew # gtdnew <- data.frame(gtd[,c(varcolXY,varlistindex)]) names(gtdnew)=names(gtd)[c(varcolXY,varlistindex)] nlin=dim(gtdnew)[1] # # Initiate variables # nXYvar <- rep(0, n.hplots) corrPearson <- rep(0, n.hplots) hs <- seq(from = lag/2., to = (lag/2.)+(n.hplots-1)*lag, by = lag) Xvar <- array(rep(NaN, n.hplots*0.5*(nlin^2)), c((nlin^2)/2.,n.hplots)) Yvar <- array(rep(NaN, n.hplots*0.5*(nlin^2)), c((nlin^2)/2.,n.hplots)) # # Compute H-data # for (j in 1:(nlin-1)) { for (i in (j+1):nlin) { dist=sqrt((gtdnew[i,1]-gtdnew[j,1])^2+(gtdnew[i,2]-gtdnew[j,2])^2) class=round((dist+0.5*lag)/lag) if (class <= n.hplots) { nXYvar[class]=nXYvar[class]+1 Xvar[nXYvar[class],class]=gtdnew[i,3] Yvar[nXYvar[class],class]=gtdnew[j,3] } } } # # Plot results # minvar=min(gtdnew[,3]) maxvar=max(gtdnew[,3]) for (i in 1:n.hplots) { corrPearson[i] <- round(cor(Xvar[1:nXYvar[i],i], Yvar[1:nXYvar[i],i], method="pearson"),digits=ndecimals) plot (Xvar[1:nXYvar[i],i], Yvar[1:nXYvar[i],i], xlim=c(minvar,maxvar), ylim=c(minvar,maxvar), col=rgb(0,0,1,0.5),pch=16,cex=sizePOINTS, xlab=names(gtdnew)[3], ylab=names(gtdnew)[3], main=paste("lag=]",lag*(i-1),";",lag*i,"] # ",nXYvar[i],"Pearson:",corrPearson[i])) } plot (hs, corrPearson, type="b", col=rgb(0,0,1,0.5), pch=16, cex=sizePOINTS, xlim=c(0,n.hplots*lag), main="Correlogram", xlab="h", ylab="Pearson") mtext(paste(titleGRAPH,names(gtdnew)[3]),side=3,cex=1.5,outer=TRUE)