Created
August 15, 2012 12:51
-
-
Save jevy/3359889 to your computer and use it in GitHub Desktop.
R code to generate graph of Water Levels in Ottawa
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
## changeing the working directory. | |
diegoDir='/Users/jevin/Code/personal/r/waterLevel' | |
setwd(diegoDir) | |
require(XML) | |
## getting the data from the web | |
x = readLines("http://www.ottawariver.ca/arnprior.htm") | |
tableNodes = getNodeSet(htmlParse(paste('</br>', x)), "//table") | |
tb = readHTMLTable(tableNodes[[1]], colCasses = "character") | |
tb2 = tb[complete.cases(tb), ] | |
## giving the data the right format | |
class(tb2) | |
names(tb2) = c("Year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec", "Annual Means", "Maximum", "Minimum", "Year") | |
head(tb2) | |
tb2 = tb2[-1, ] | |
tb2 = tb2[-c(nrow(tb2)-1, nrow(tb2)), ] | |
for( i in 1:ncol(tb2)){ | |
is.na(tb2[,i]) = which(tb2[,i] == '---') | |
} | |
tb2 = tb2[, -ncol(tb2)] | |
## saving data as csv | |
write.table(tb2, 'waterLevels.csv', sep=';', row.names=FALSE) | |
#################### | |
rm(list = ls()) | |
#################### | |
wlevel = read.table('waterLevels.csv', sep = ';', header = TRUE) | |
pdf('DensityEstimationHistoricWaterLevels.pdf') | |
mnths = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec") | |
i = 1 | |
plot(density(wlevel[, mnths[i]], na.rm=TRUE), col = i, xlim = c(min(wlevel[, mnths], na.rm = TRUE), max(wlevel[, mnths], na.rm = TRUE)) ) | |
for(i in 2:length(mnths)){ | |
lines(density(wlevel[, mnths[i]], na.rm=TRUE), add = TRUE, col = i) | |
} | |
legend(x="topright", legend = mnths, col = 1:length(mnths), lwd = 1) | |
dev.off() | |
## current water Level | |
current = readLines("http://www.ottawariver.ca/emess.htm") | |
current = substr(current[grep('Chats Lake at Arnprior', current)], 91, 95) | |
actualMonth = months(Sys.Date(), abb=TRUE) | |
historic = wlevel[, actualMonth] | |
pdf('ActualAgainstHistoricWaterLevel.pdf') | |
plot(density(historic, na.rm = TRUE), main = paste('Density estimation of historical water level at', actualMonth)) | |
abline(v=current, col = 'red') | |
legend(x="topright", legend = "Current water level", col = 'red', lwd = 1) | |
dev.off() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment