#This is the source file for the graph
#http://home.casema.nl/errenwijlens/co2/co2lt_2007.gif
#on
#http://home.casema.nl/errenwijlens/co2/sink.htm

#Source data is from 
#http://vortex.nsstc.uah.edu/data/msu/t2lt/tltglhmam_5.2
#and
#http://www.esrl.noaa.gov/gmd/ccgg/trends/co2_mm_mlo.dat
#Nov 2007: CO2 data source changed to
#ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.dat
#Dec 2007: CO2 data source changed to
#ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt

#Script and graphs can be used freely for non profit purposes with citation of the following source:
#"Graphs/script made available by Hans Erren from http://home.casema.nl/errenwijlens/co2/co2_lt_noaa.R"

#© Hans Erren, July 2007, updated October 2007, updated November 2007, updated December2007, updated June 2008, updated February 2009




#Read Lower Tropospheric temperature from University of Alabama in Huntsville
lt2t <- read.table("http://vortex.nsstc.uah.edu/data/msu/t2lt/tltglhmam_5.2", header=TRUE,  skip=4, nrows=362, na.strings = "-99.990")

#calculate decimal time
lt2t$time <- lt2t$YEAR+(lt2t$MON-0.5)/12

# show last part of data
tail(lt2t)

#plot data
x<-lt2t$time
y<-lt2t$GLOBAL
plot (x, y, ann=FALSE,  type = "n")
points (x, y, bg = "limegreen", pch = 21)
y<-lt2t$GLOBAL.1
lines(x,y, col = "blue", lwd=2)
 title(main = "Lower Troposphere Temperature by UAH", xlab = "Year",ylab = "Temperature (K)",  
    col.main = "blue", col.lab = "black", cex.main = 1.2, cex.lab = 1, 
    font.main = 4, font.lab = 3)
readline("press return to see next plot")


#Read Monthly CO2 data Mauna Loa from NOAA server
#Old html file structure
#mloco2 <- read.table("http://www.esrl.noaa.gov/gmd/webdata/ccgg/trends/co2_mm_mlo.dat", header=FALSE,  skip=43, nrows=583, na.strings = "-99.990")
#read data post 1977
#mloco2 <- read.table("http://www.esrl.noaa.gov/gmd/webdata/ccgg/trends/co2_mm_mlo.dat", header=FALSE,  skip=263, na.strings = "-99.990")
#colnames(mloco2)<- c("OBS", "YEAR", "MON", "PPMCO2")
#mloco2$time <- mloco2$YEAR+(mloco2$MON-0.5)/12
#New ftp file structure
#mloco2 <- read.table("ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt", header=FALSE,  skip=290,nrows=377, na.strings = "-99.99")
#colnames(mloco2)<- c("YEAR", "MON","time", "PPMCO2","PPMCO2_int","PPMCO2_trend")
#New ftp file structure
mloco2 <- read.table("ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt", header=FALSE,  skip=307,nrows=374, na.strings = "-99.99")
colnames(mloco2)<- c("YEAR", "MON","time", "PPMCO2","PPMCO2_int","PPMCO2_trend", "days")



#test plot
x<-mloco2$time
y<-mloco2$PPMCO2
plot (x, y, ann=FALSE,  type = "n")
lines(x,y, col = "blue", lwd=2)
 title(main = "Mauna Loa CO2 measurements", xlab = "Year",ylab = "CO2 (ppm)",  
    col.main = "blue", col.lab = "black", cex.main = 1.2, cex.lab = 1, 
    font.main = 4, font.lab = 3)
readline("press return to see next plot")

# calculating annual CO2 difference
mloco2$delCO2<-(mloco2$PPMCO2[mloco2$time > 1978.875])-mloco2$PPMCO2 

#set 12 months data tail to void (clumsy method)

mloco2$delCO2[363]<-NA
mloco2$delCO2[364]<-NA
mloco2$delCO2[365]<-NA
mloco2$delCO2[366]<-NA
mloco2$delCO2[367]<-NA
mloco2$delCO2[368]<-NA
mloco2$delCO2[369]<-NA
mloco2$delCO2[370]<-NA
mloco2$delCO2[371]<-NA
mloco2$delCO2[372]<-NA
mloco2$delCO2[373]<-NA
mloco2$delCO2[374]<-NA


#show data head
head(mloco2, n=24)

#show data tail
tail(mloco2, n=24)



#testplot
   x<-c(1978,2009)
   y<-c(0,4)
   plot (x, y, ann=FALSE,  type = "n")
   x<-mloco2$time
   y<-mloco2$delCO2
   lines(x,y, col = "blue", lwd=1)
   #points (x, y, bg = "limegreen", pch = 21, cex=0.3)
    title(main = "Mauna Loa CO2 measurements annual increment", xlab = "Year",ylab = "CO2 (ppm)",  
       col.main = "blue", col.lab = "black", cex.main = 1.2, cex.lab = 1, 
       font.main = 4, font.lab = 3)
readline("press return to see next plot") 

#plot co2 increment on last month
mloco2$time2<- mloco2$time +1
       
#Final plot       
plot(lt2t$time,((lt2t$GLOBAL+0.3)*4), type="n", ann=FALSE)
lines(lt2t$time,((lt2t$GLOBAL+0.3)*4),type="l", col="blue")
lines(mloco2$time2,(mloco2$delCO2), col="black")

title(main = "Comparison of annual CO2 growth and lower troposphere temperature \n Black: Moving annual CO2 growth (Mauna Loa), \nBlue: Lower troposphere temperature (MSU UAH)", xlab = "year",ylab = "ppm CO2/year , Temperature x 4",  
    col.main = "black", col.lab = "black", cex.main = 0.7, cex.lab = 0.8, 
    font.main = 2, font.lab = 4)
 
legend (1989,-0.7, "R-script at http://home.casema.nl/errenwijlens/co2/co2_lt_noaa.R", cex=0.7, bty="n") 

readline("press return to see alternatives")
 
##testing alternatives:
##plot co2 increment midyear
#    mloco2$time2<- mloco2$time +0.5
#    lines(mloco2$time2,(mloco2$delCO2), col="red")     
##plot on co2 increment on first month 
#    mloco2$time2<- mloco2$time 
#    lines(mloco2$time2,(mloco2$delCO2), col="green")