Is it possible to add lines/text/minor tick marks/personalized axes to a variogram plot?

Question

I'm trying to add text, lines, minor tick marks, and/or personalized axes to a variogram using the plot.variogram function. When I try to add any of these using the traditional function calls (i.e. text("Text Here")) it returns the error of plot.new has not been called yet even though the plotting window for the variogram is open.

Here is my code:

#v is sample variogram, vmf is fitted model
plot(v, model=vmf, xlim=c(0, 65), ylim=c(0,25), xlab="Distance between Point Pairs (km)",
ylab="Semivariance ((C/km) )", cex.xlab=6.5, cex.ylab=6.5, cex.xaxis=2.5, cex.main=5.5)

#Add a 2 to the y label that is in 10 pt. font so it looks like it is (C/km)^2
par(ps=10, cex=1, cex.main=1)
text(-2, 16, labels=2, srt=90)

#Add lines showing the desired point pair distance and semivariance for the problem
par(new=TRUE, lines(c(53,53),c(0,15),col='red'))
par(new=TRUE, lines(c(0, 53),c(15,15),col='red'))

#Add axis minor tick marks in increments of 5
axis(side=1, at=c(0, 5, 15, 25, 35, 45, 55, 65), labels = NA, tck=-0.01, pos=0) 
axis(side=2, at=c(0, 2.5, 7.5, 12.5, 17.5, 22.5, 25),labels = NA, tck=-0.01, pos=0)

I have tried to "trick" R by calling:

plot(c(0,65), c(0,25))

and then running the code above. This allows for the traditional functions to work, but they are unfortunately not in the appropriate locations (i.e. x=5 is not located at 5 on the x axis).

Any recommendations for better ways to "trick" R to plotting correctly? Any functions that add text, axes, etc. automatically to variogram plots?

Please let me know if there's anything else you would like to know.

Thanks!

Answer 1

I'm not sure how you got your variogram, but using some of the code and info from this link, I was able to get something that may help you.

Using the geoR package, function variog you can manipulate the plot as usual.

> sampleV <- 
    read.table(header = TRUE, text = "Station   Av8top      Lat       Lon
  1       60 7.225806 34.13583 -117.9236
  2       69 5.899194 34.17611 -118.3153
  3       72 4.052885 33.82361 -118.1875
  4       74 7.181452 34.19944 -118.5347
  5       75 6.076613 34.06694 -117.7514
  6       84 3.157258 33.92917 -118.2097
  7       85 5.201613 34.01500 -118.0597
  8       87 4.717742 34.06722 -118.2264
  9       88 6.532258 34.08333 -118.1069
  10      89 7.540323 34.38750 -118.5347", row.names = 1)

> library(geoR)
> sampleVMF <- variog(coords = sampleV[,3:4], data = sampleV[,2], 
                      breaks = seq(0, 1.5, length = 11))
> plot(sampleVMF, axes = FALSE,  
       xlab="Distance between Point Pairs (km)",
       ylab="Semivariance ((C/km) )")
> axis(1, at = sampleVMF$u)
> axis(2, at = sampleVMF$v)
> box()
> text(median(sampleVMF$u), median(sampleVMF$v), "Hello world!")
> lines(sampleVMF$u, sampleVMF$v)

Answer 2

Richard Scriven's answer works well when using the geoR package.

If using the gstat package, the plots will all need to be modified using trellis graphics commands (lattice package and latticeExtra package). The reason for the error of plot.new has not been called yet is because of base graphics use on trellis graphics.

Example plot modifications:

Create the plot using list to modify all input parameters. The scales modifies the x and y axis, and you can add in tick marks (at) without labels (see y-axis list). The key adds the legend.

p1 = plot(v.GF, model=vmf.GF, lwd=2, col.line="black", pch=16, cex=1.5, ylim=c(0,25), xlim=c(0,150), 
      main=list("Gaussian Semivariogram Model for Geothermal Gradient", cex=2),
      xlab=list("Distance between Point Pairs (km)", cex=2), 
      ylab=list(expression("Geothermal Gradient Semivariance (°C/km)"^2), cex=2), 
      scales=list(x=list(at=c(0,25,50,75,100,125,150), labels=c(0,25,50,75,100,125,150)), y=list(at=c(0,2.5,5,7.5,10,12.5,15,17.5,20,22.5,25), labels=c(0,"",5,"",10,"",15,"",20,"",25)), cex=2), 
      key=list(text=list(lab=c("Gaussian Model","Sill","Maximum Interpolation Distance")), space="top", lines=list(col=c("black","black","red"), lwd=2, lty=c(1,2,1)), columns=3, cex=1.5))  

Then plot and modify the plotting area using trellis.focus, and llines to add lines and ltext to add text:

trellis.focus("panel",1,1)
plot(p1)
trellis.focus("panel",1,1)
llines(x=c(50,50), y=c(0,14.5), col="red", lwd=2, lty=1)
llines(x=c(0,50), y=c(14.5,14.5), col="red", lwd=2, lty=1)
llines(x=c(0,150), y=c(vmf.GF$psill[2]+vmf.GF$psill[1],vmf.GF$psill[2]+vmf.GF$psill[1]), col="black", lty=2, lwd=2)
ltext(x=12,y=5,"Nugget ~6.0", cex=1.5)
trellis.unfocus()

Here is the resulting plot: