################################################################################### # # Fit theoretical models of variograms (ellipse fitting) # # Last review: 23/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### varcolXY <- c(1,2) varlistindex <- 6 fileNameOutMOD="C:/Users/Utilizador/Desktop/GTD/Dados/Solos_MP/Co.MOD" ndecimals <- 3 resolucao <- 200 unidadesdist="m" # # Models are inputed by proportion of variance (sum of Nugget effect + Sills = 1) # nugget_effect <- 0.05 nstructures <- 1 model_type <- c("Exp","Exp") sill <- c(0.95,1.0) range <- c(1800,500) anisotropy <- c(1.5,1.0) maindir <- c(60,90) maxX=2500 ################################################################################### par(mfrow=c(2,3), mar=c(3.5, 3.5, 2, 1), mgp=c(2.25, 0.8, 0)) # # Duplicate gtd -> gtdnew # gtdnew=gtd[,c(varcolXY,varlistindex)] names(gtdnew)=names(gtd)[c(varcolXY,varlistindex)] nlin=dim(gtdnew)[1] # # Initiate variables # anguloref = c(NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN) angulotol = c(NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN) legenda <- array(rep("",times=10),dim=c(10)) # # Read a variogram MOD file # con1=file(fileNameOutMOD,open="r") linha <- scan(file=con1, nlines=1, what=character(), quiet=TRUE) linha <- scan(file=con1, nlines=1, what=character(), quiet=TRUE) linha <- scan(file=con1, nlines=1, what=character(), quiet=TRUE) varname <- trimws(scan(file=con1, nlines=1, what=character(), quiet=TRUE)) nang <- scan(file=con1, nlines=1, what=integer(), quiet=TRUE) maxnpassos <- scan(file=con1, nlines=1, what=integer(), quiet=TRUE) for (i in 1:nang) { angulos <- scan(file=con1, nlines=1, sep = ";", what=character(), quiet=TRUE) anguloref[i]=strtoi(trimws(angulos[1])) angulotol[i]=strtoi(trimws(angulos[2])) } distanciah <- array(rep(0,times=nang*maxnpassos),dim=c(nang,maxnpassos)) variograma <- array(rep(0,times=nang*maxnpassos),dim=c(nang,maxnpassos)) npares <- array(rep(0,times=nang*maxnpassos),dim=c(nang,maxnpassos)) variancia <- scan(file=con1, nlines=1, what=single(), quiet=TRUE) for (i in 1:nang) { for (j in 1:maxnpassos) { dadosMOD <- scan(file=con1, nlines=1, what=character(), sep = " ", quiet=TRUE) distanciah[i,j]=as.numeric(dadosMOD[1]) variograma[i,j]=as.numeric(dadosMOD[2]) npares[i,j]=as.numeric(dadosMOD[3]) if (npares[i,j] == 0) { distanciah[i,j]=NaN variograma[i,j]=NaN } } } close(con1) # # Compute ellipses # rrr=range/anisotropy theta <- seq(0,359) xx <- rep(NaN, times=360) yy <- rep(NaN, times=360) xxx <- array(rep(0, 2*360), c(2,360)) yyy <- array(rep(0, 2*360), c(2,360)) rangecorr <- array(rep(0, 2*2*360), c(2,2,360)) if (nstructures > 0) { for(k in 1:nstructures){ aroda <- maindir[k]*pi/180 yy[1:90] <- sqrt(1/((tan(theta[1:90]*pi/180)^2/rrr[k]^2)+(1/range[k]^2))) yy[91:180] <- -sqrt(1/((tan(theta[91:180]*pi/180)^2/rrr[k]^2)+(1/range[k]^2))) yy[181:270] <- -sqrt(1/((tan(theta[181:270]*pi/180)^2/rrr[k]^2)+(1/range[k]^2))) yy[271:360] <- sqrt(1/((tan(theta[271:360]*pi/180)^2/rrr[k]^2)+(1/range[k]^2))) xx[1:180] <- sqrt((1-(yy[1:180]^2/range[k]^2))*rrr[k]^2) xx[181:360]<- -sqrt((1-(yy[181:360]^2/range[k]^2))*rrr[k]^2) xxx[k,] <- xx*cos(aroda)+yy*sin(aroda) yyy[k,] <- -xx*sin(aroda)+yy*cos(aroda) if (maindir[k] == 0) { rangecorr[k,1,] <- c(0:359) } else { rangecorr[k,1,] <- c(maindir[k]:359,0:(maindir[k]-1)) } rangecorr[k,2,] <- sqrt(xxx[k,]^2+yyy[k,]^2) rangecorr[k,,] <- rangecorr[k,,order(rangecorr[k,1,])] } } # # Compute models # maxdist=max(range) hhmodelo <- rep(0,times=resolucao) variogmodelo <- array(rep(0,times=10*resolucao),dim=c(10,resolucao)) for(j in 1:resolucao){ hhmodelo[j]=(j-1)*maxdist/resolucao } for (i in 1:nang) { azim = as.integer(ifelse(anguloref[i] < 0, anguloref[i]+360, anguloref[i]))+1 for(j in 1:resolucao){ variogmodelo[i,j]=variancia*nugget_effect if (nstructures > 0) { for(k in 1:nstructures){ if (angulotol[i] < 90) { if (model_type[k] == "Sph") { variogmodelo[i,j]= variogmodelo[i,j]+ ifelse (hhmodelo[j] <= rangecorr[k,2,azim], variancia*sill[k]*(((3*hhmodelo[j])/(2*rangecorr[k,2,azim]))-0.5*(hhmodelo[j]/rangecorr[k,2,azim])^3), variancia*sill[k]) } else if (model_type[k] == "Exp") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*(1-exp(-3*hhmodelo[j]/rangecorr[k,2,azim]))) } else if (model_type[k] == "Gauss") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*(1-exp(-3*hhmodelo[j]^2/rangecorr[k,2,azim]^2))) } else if (model_type[k] == "Power") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*hhmodelo[j]^rangecorr[k,2,azim]) } } else { if (model_type[k] == "Sph") { variogmodelo[i,j]= variogmodelo[i,j]+ ifelse (hhmodelo[j] <= mean(rangecorr[k,2,]), variancia*sill[k]*(((3*hhmodelo[j])/(2*mean(rangecorr[k,2,])))-0.5*(hhmodelo[j]/mean(rangecorr[k,2,]))^3), variancia*sill[k]) } else if (model_type[k] == "Exp") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*(1-exp(-3*hhmodelo[j]/mean(rangecorr[k,2,])))) } else if (model_type[k] == "Gauss") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*(1-exp(-3*hhmodelo[j]^2/mean(rangecorr[k,2,])^2))) } else if (model_type[k] == "Power") { variogmodelo[i,j] = variogmodelo[i,j] + variancia*(sill[k]*hhmodelo[j]^mean(rangecorr[k,2,])) } } } } } } # # Draw plots # cores <- c("red2","orange2","blue2","sienna2","green2","maroon2","black") maxdist=round(max(distanciah,na.rm=TRUE)) maxscale=round(max(variograma,na.rm=TRUE),3) plot(gtdnew[,1], gtdnew[,2], pch=20, lwd=0.5, asp=1, col="gray75", cex=0.75, main="Localizações amostras", xlab=names(gtdnew)[1], ylab=names(gtdnew)[2]) cgx=mean(gtdnew[,1]) cgy=mean(gtdnew[,2]) for (i in 1:nang) { azim = as.integer(ifelse(anguloref[i] < 0, anguloref[i]+360, anguloref[i]))+1 if (angulotol[i] < 90) { bx=cgx-sum(rangecorr[,2,azim])*cos((90-anguloref[i])*pi/180.0) by=cgy-sum(rangecorr[,2,azim])*sin((90-anguloref[i])*pi/180.0) cx=cgx+sum(rangecorr[,2,azim])*cos((90-anguloref[i])*pi/180.0) cy=cgy+sum(rangecorr[,2,azim])*sin((90-anguloref[i])*pi/180.0) lines(c(bx,cx),c(by,cy),col=cores[i],lwd=2) } } if (nstructures > 0) { for(k in 1:nstructures){ if (model_type[k] != "Power") { lines (xxx[k,]+cgx, yyy[k,]+cgy) lines (xxx[k,c(1,181)]+cgx,yyy[k,c(1,181)]+cgy,col="black",lty=1) lines (xxx[k,c(91,271)]+cgx,yyy[k,c(91,271)]+cgy,col="black",lty=1) } } } # # Text of the models # for (i in 1:nang) { azim = as.integer(ifelse(anguloref[i] < 0, anguloref[i]+360, anguloref[i]))+1 legenda[i]=paste("Co:",variancia*nugget_effect) if (nstructures > 0) { for(k in 1:nstructures){ if (model_type[k] == "Power") { legenda[i]=paste(legenda[i],"+",model_type[k],"(alfa=",sill[k],"; beta=",round(rangecorr[k,2,azim]),")") } else { legenda[i]=paste(legenda[i],"+",model_type[k],"(C=",variancia*sill[k],"; a=",round(rangecorr[k,2,azim]),")") } } } } # # Graphics # for (n in 1:nang) { plot(distanciah[n,],variograma[n,],type="b",xlab="",ylab="",xlim=c(0,maxX),ylim=c(0,maxscale),col=cores[n],axes=FALSE) axis(side=1, at=seq(0, maxdist, by=maxdist/10.0), pos=0) axis(side=2, at=seq(0, maxscale, by=maxscale), pos=0) lines (c(0,maxdist),c(variancia,variancia),col="red",lwd=2) lines (c(0,distanciah[n,1]),c(0,variograma[n,1]),col=cores[n],lwd=2,lty=3) lines (hhmodelo,variogmodelo[n,],col="gray50",lwd=2,lty=1) legend(28,2500, legenda[n], lty=1, col="gray50", bty='n', lwd=2, cex=1) title(main = paste("Variograma",names(gtdnew)[3]," max{Np}:",max(npares[n,],na.rm=TRUE)),xlab="h",ylab="gama(h)",cex.main=1.25) }