R code to load Whitetop data from Neyman et al (1969) and recreate plots (Fig 2)

create-Neyman1969-plots.R

url <- 'https://gist.githubusercontent.com/richarddmorey/862ca2681afd3cd85b3b/raw/8de4af8d99cbff8745f7b52fd9009410ff450bd5/Neyman-etal-1969-whitetop.txt'
library(RCurl)
df <- getURL(url, ssl.verifypeer=FALSE)
dat.all <- read.table(textConnection(df))

## Break into three data frames
dat.freq.wet = dat.all[,c(1:2,3:6)]
dat.rainfall.wet = dat.all[,c(1:2,7:10)]
dat.rainfall.all = dat.all[,c(1:2,11:14)]

colnames(dat.freq.wet) = c("Region","# gages", "Seeded (%)", "Not Seeded (%)","Diff. (%)","p")
colnames(dat.rainfall.wet) = c("Region","# gages", "Seeded (inches)", "Not Seeded (inches)","Change (%)","p")
colnames(dat.rainfall.all) = c("Region","# gages", "Seeded (inches)", "Not Seeded (inches)","Change (%)","p")

### quick Monte Carlo function to get average distance in rings
avg.dist = function(N = 100000){
  v = runif(N,-180,180)
  u = runif(N, -180, 180)
  x = sqrt(u^2 + v^2)
  r = x[x<180]
  z = sapply(r, function(z) sum( z>(30*1:6) ))
  tapply(r,z,mean)
}

xx = avg.dist()

## Wet days

matplot(xx,dat.rainfall.wet[1:6,3:4],ylim=c(0,.3),ty='b',pch=19,
        col=c("red","black"),xlim=c(0,190), axes=FALSE,
        ylab="Average precipitation",lty=1,xlab="Average distance from center (miles)")
axis(3,at=xx,lab=LETTERS[1:6])
axis(1,at=1:6*30)
box()
#abline(h=dat.rainfall.wet[7,3],col="red",lty=2)
#abline(h=dat.rainfall.wet[7,4],col="black",lty=2)

## I believe there is a problem with the data in the table. The 
## numbers for the entire region on wet days do not make sense 
## given the individual region numbers. It should be a weighted
## average of the regions, and that works out in the "all" days
## table, but the numbers are way off for the "wet" days.

axis(2)
mtext("Wet days",3,-1.5,adj=.9,cex=1.3)
legend(0,.1,legend=c("Non-seeded","Seeded"),col=c("black","red"),lty=1,pch=19)

text(xx,dat.rainfall.wet[-7,3]-.03,paste0("p=",dat.rainfall.wet[-7,6]))
mtext(paste0("Overall p=",dat.rainfall.wet[7,6]),1,-1.5,adj=.9,cex=1.3)

## All days


matplot(xx,dat.rainfall.all[1:6,3:4],ylim=c(0,.3),ty='b',pch=19,
        col=c("red","black"),xlim=c(0,190), axes=FALSE,
        ylab="Average precipitation",lty=1,xlab="Average distance from center (miles)")
axis(3,at=xx,lab=LETTERS[1:6])
axis(1,at=1:6*30)
box()
#abline(h=dat.rainfall.all[7,3],col="red",lty=2)
#abline(h=dat.rainfall.all[7,4],col="black",lty=2)
axis(2)
mtext("All days",3,-1.5,adj=.9,cex=1.3)
legend(0,.1,legend=c("Non-seeded","Seeded"),col=c("black","red"),lty=1,pch=19)
text(xx,dat.rainfall.all[-7,4]+.03,paste0("p=",dat.rainfall.all[-7,6]))
mtext(paste0("Overall p=",dat.rainfall.all[7,6]),1,-1.5,adj=.9,cex=1.3)