################################################################################### # # Draw elippses for variogram modelling purposes # # Last review: 23/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### nstructures <- 2 range <- c(2.0, 5.0) anisotropy <- c(1, 2) maindir <- c(0, 120) maindir <- ifelse (maindir < 0, maindir+180, maindir) cor_elipse <- c("blue","red") ################################################################################### par(mfrow=c(1,1), mar=c(3.5, 3.5, 2, 1), mgp=c(2.25, 0.8, 0)) # # 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,])] } } # # Graphics # rangegraf=1.1*max(max(xxx)-min(xxx),max(yyy)-min(yyy)) plot (0, 0, xlim=c(-rangegraf/2,rangegraf/2), ylim=c(-rangegraf/2,rangegraf/2), xlab="X", ylab="Y", col=rgb(0,0,1,0.5), pch=19, cex=0.1, asp=1) abline(h=0, lwd=1, col = "gray50", lty = 3) abline(v=0, lwd=1, col = "gray50", lty = 3) if (nstructures > 0) { for(k in 1:nstructures){ lines (xxx[k,], yyy[k,], col=cor_elipse[k], lty=1, lwd=2) lines (c(xxx[k,1],xxx[k,359]), c(yyy[k,1],yyy[k,359]), col=cor_elipse[k], lty=1, lwd=2) lines (c(xxx[k,1],xxx[k,181]), c(yyy[k,1],yyy[k,181]), col=cor_elipse[k], lty=1, lwd=1) lines (c(xxx[k,91],xxx[k,271]), c(yyy[k,91],yyy[k,271]), col=cor_elipse[k], lty=1, lwd=1) } } title(main = "Elipse de um modelo de variograma", xlab="X", ylab="Y", cex.main = 1)