################################################################################### # # Compute variograms for 1D grids # # Last review: 23/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### fileNameGrid="C:/temp/DemoGrid1d.prn" fileNameMOD="C:/temp/Demogrid1D.MOD" xmincenter=0 ncellsx=12 ncellsy=1 cellsizex=2 titleX="Coordinate X" legendaeixoX="Distance (m)" titleGRAPH="Water depth" ################################################################################### maxcells <- max(ncellsx) xgrid <- array(rep(NaN, ncellsx), c(ncellsx)) ygrid <- array(rep(0, ncellsx), c(ncellsx)) vgrid <- array(rep(NaN, ncellsx), c(ncellsx)) variograma <- array(rep(0, maxcells), c(maxcells)) npares <- array(rep(0, maxcells), c(maxcells)) hh <- array(rep(0, maxcells), c(maxcells)) # # Read 1D grid file # con1=file(fileNameGrid,open="r") for (iy in 1:ncellsy) { for (ix in 1:ncellsx) { cc <- scan (con1,nlines=1,quiet = TRUE) xgrid[ix] <- xmincenter+(ix-1)*cellsizex vgrid[ix] <- cc[1] } } close(con1) # # Variograms in OX direction # for (iy in 1:ncellsy) { for (ix in 1:ncellsx) { if (is.finite(vgrid[ix])) { for (ih in 1:maxcells) { if (ix+ih <= ncellsx) { if (is.finite(vgrid[ix+ih])) { variograma[ih]=variograma[ih]+(vgrid[ix]-vgrid[ix+ih])^2 npares[ih]=npares[ih]+1 } } } } } } for (j in 1:maxcells) { hh[j]=cellsizex+cellsizex*(j-1) variograma[j] <- ifelse(npares[j] > 0, variograma[j]/(2*npares[j]), NaN) } # # Gráficos # par(mfrow=c(1,2), mar=c(3.5, 3.5, 2, 1), mgp=c(2.25, 0.8, 0)) variancia=sd(vgrid,na.rm=TRUE)^2 maxgama=max(variograma,variancia,na.rm = TRUE) minx=min(xgrid) maxx=max(xgrid) dx=maxx-minx minx=minx-0.01*dx maxx=maxx+0.01*dx minv=min(vgrid) maxv=max(vgrid) dv=maxv-minv minv=minv-0.01*dv maxv=maxv+0.01*dv plot(xgrid, ygrid, type="p", xlim=c(minx,maxx), ylim=c(0,0), pch=19, col="red", main=titleGRAPH, xlab=titleX, ylab="", cex.main=1.0, cex.sub=0.5, cex.lab=1.0, cex.axis=1.0) text(xgrid, ygrid, labels=vgrid, pos=3, col="black", cex=0.5) plot(hh,variograma,type="p",xlim=c(0,dx*1.1),ylim=c(0,maxgama*1.1),pch=20, main = "Variograma experimental 1D", xlab=legendaeixoX, ylab="Variograma", cex.main=1.0, cex.sub=0.5, cex.lab=1.0, cex.axis=1.0) abline(h=variancia, col="red", lwd=2) text(hh, variograma, labels=npares, pos=3, col="black", cex=0.8) # # Output a variogram MOD file # con1=file(fileNameMOD,open="w") writeLines("VARG", con=con1, sep = "\n") writeLines("geoVAG", con=con1, sep = "\n") writeLines("1", con=con1, sep = "\n") writeLines(varNAME, con=con1, sep = "\n") writeLines("1", con=con1, sep = "\n") writeLines(paste(maxcells), con=con1, sep = "\n") writeLines(paste("90",";","0"), con=con1, sep = "\n") writeLines(paste(variancia), con=con1, sep = "\n") for (j in 1:maxcells) { writeLines(paste(hh[j],variograma[j],npares[j]), con=con1, sep = "\n") } close(con1)