Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
R code to load Whitetop data from Neyman et al (1969) and recreate plots (Fig 2)
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.