Last active
November 12, 2015 15:44
-
-
Save richarddmorey/5f22fc742535d078c25f to your computer and use it in GitHub Desktop.
R code to load Whitetop data from Neyman et al (1969) and recreate plots (Fig 2)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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