# Oil Production Graphs
# Rev 20 December 2015
# Michael Minn

data = read.csv("oil-production.csv", as.is=T, skip=8)
years = data$YEAR
data = data * 365 / 1000 # convert to MM barrels/year
data$YEAR = years

# Fill out data unavailable from EIA with data from USCB and Alaska Tax
data$LOWER48 = ifelse(!is.na(data$LOWER48), data$LOWER48, data$USCB_PETROLEUM)
data$NGPL = ifelse(!is.na(data$NGPL), data$NGPL, data$USCB_NGPL)
data$ALASKA = ifelse(!is.na(data$ALASKA), data$ALASKA, data$ALASKA_TAX)


# Conventional is Lower 48 minus tight oil and natural gas plant liquids
data$TIGHT[is.na(data$TIGHT)] = 0
data$CONVENTIONAL = data$LOWER48 - data$TIGHT
data$CONVENTIONAL[data$YEAR >= 2014] = NA

png(filename="oil-production-conventional.png", width=600, height=400, pointsize=12)

plot(data$YEAR, data$CONVENTIONAL, type="l", lwd=3, col="navy", fg="white", 
	ylim=c(0,4000), xlim=c(1860, 2060), las=1, 
	xlab="Year", ylab="MM Barrels / Year", main="US Mainland Conventional Petroleum Production",
	panel.first=grid(nx=NA, ny=NULL, lwd=1, col='gray', lty='solid'))

lines(c(1850, 2100), c(0, 0), lwd=3, col="black") # X-axis

# Fit a logistic distribution curve to conventional production
# Peak year 1972, Ultimate production 208.2 Billion Barrels
# model = nls(CONVENTIONAL ~ SSlogis(YEAR, peak, half, scale), data=data, start=list(peak = 3400, half=1970, scale=70))

model = nls(CONVENTIONAL ~ ultimate * dlogis(YEAR, half, scale), data=data, start=list(ultimate = 100000, half=1970, scale=7))
fityears = 1850:2100
fitted = predict(model, data.frame(YEAR = fityears))
lines(fityears, fitted, lwd="3", col="darkred", lty="dotted")


# Hubbert's curves
fit150 = 150000 * dlogis(fityears, 1965, 14)
lines(fityears, fit150, lwd="3", col="darkgreen", lty="dashed")

fit200 = 200000 * dlogis(fityears, 1970, 16.7)
lines(fityears, fit200, lwd="3", col="orange", lty="dotdash")

legend(x="topleft", legend=c("Actual (EIA 2015, USCB 1976)", 
	"Fitted Logistic Model", "Hubbert (1956) 150 BBL", "Hubbert (1956) 200 BBL"),
	lwd=3, col=c("navy", "darkred", "darkgreen", "orange"),
	lty=c("solid", "dotted", "dashed", "dotdash"), bg="white")

dev.off()

# Total US Production using current EIA accounting

png(filename="oil-production-total.png", width=600, height=400, pointsize=12)

data$US = data$LOWER48 + ifelse(is.na(data$NGPL), 0, data$NGPL) + ifelse(is.na(data$ALASKA), 0, data$ALASKA)
plot(data$YEAR, data$CONVENTIONAL, type="l", lwd=3, col="navy", fg="white", 
	ylim=c(0,6000), xlim=c(1860, 2060), las=1, 
	xlab="Year", ylab="MM Barrels / Year", main="Historic and Projected US Liquid Fuels Production",
	panel.first=grid(nx=NA, ny=NULL, lwd=1, col='gray', lty='solid'))

lines(data$YEAR, data$US, lwd=3, col="darkred")
# lines(data$YEAR, data$US_BP, lwd=3, col="darkgreen")
lines(fityears, fit200, lwd="3", col="orange", lty="dotdash")

# Reference case projections
ref = data[!is.na(data$REF_US_TOTAL), c("YEAR", "REF_US_TOTAL")]
lines(ref$YEAR, ref$REF_US_TOTAL, lwd=3, lty="dotted", col="darkgreen")

# High growth projections
data$PROJ_HI = data$PROJ_HI_LOWER48 + data$PROJ_HI_ALASKA + data$PROJ_HI_NGPL
hi = data[!is.na(data$PROJ_HI), c("YEAR", "PROJ_HI")]
lines(hi$YEAR, hi$PROJ_HI, lwd=3, lty="dashed", col="red")

legend(x="topleft", legend=c("Conventional (EIA 2015, USCB 1976)", 
	"US Total (EIA 2015, USCB 1976)",
	"Hubbert (1956) 200 BBL Projection", "EIA (2015) Reference Case Projection", 
	"EIA (2015) High Oil Price Projection"), 
	lwd=3, col=c("navy", "darkred", "orange", "darkgreen", "red"),
	lty=c("solid", "solid", "dotdash", "dotted", "dashed"), bg="white")

dev.off()


# Global Oil Production

png(filename="oil-production-world.png", width=600, height=400, pointsize=12)
# par(mar=c(5, 6, 4, 2))

plot(data$YEAR, data$WORLD_BP, type="l", lwd=3, col="navy", fg="white", 
	ylim=c(0,40000), xlim=c(1900, 2100), las=1, 
	xlab="Year", ylab="MM Barrels / Year", main="Global Oil Production",
	panel.first=grid(nx=NA, ny=NULL, lwd=1, col='gray', lty='solid'))

lines(c(1850, 2100), c(0, 0), lwd=3, col="black") # X-axis

# Hubbert (1969) curves
fit1350 = 1350000* dlogis(fityears, 1990, 13)
lines(fityears, fit1350, lwd="3", col="darkgreen", lty="dashed")

fit2100 = 2100000* dlogis(fityears, 2000, 14)
lines(fityears, fit2100, lwd="3", col="orange", lty="dotdash")

legend(x="topleft", bg="white", legend=c("Total World (BP 1951; 2015)",
	"Hubbert (1969) 1,350 BBL Proj.", 
	"Hubbert (1969) 2,100 BBL Proj."),
	lwd=3, col=c("navy", "darkgreen", "orange"),
	lty=c("solid", "dashed", "dotdash"))

dev.off()