%let name=co2;
filename odsout '.';

/*
http://en.wikipedia.org/wiki/Image:Mauna_Loa_Carbon_Dioxide.png
http://www.cmdl.noaa.gov/projects/src/web/trends/co2_mm_mlo.dat
http://www.esrl.noaa.gov/gmd/ccgg/trends/
*/

%let urlname=ftp://ftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt;
*%let urlname=ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt;

%include '../data_vault/proxy.sas';
/*
I %include a file that sets the my_proxy macro variable.
You won't be able to use my proxy, therefore set your own
macro variable, which will be something like the following:
%let my_proxy=http://yourproxy.com:80;
*/

filename dataurl url "&urlname"
 proxy="&my_proxy";

/* or, if you need to use a snapshot copy, saved onto disk */
/*
filename dataurl "co2_mm_mlo.txt" termstr=crlf;
*/

data co2data;
label co2='Carbon dioxide concentration (ppmv)';
format month z2.;
infile dataurl dlm=' ' firstobs=73 lrecl=256 pad;
input year month decimal_date average_co2 interpolated_co2 season_corrected_co2;
co2=average_co2;
if (co2 eq -99.99) then co2=.;
/* Always the same 'MLO' in this data (Mauna Loa, Hawaii) */
site='MLO';
x_variable=year+((month-1)/12);

datestring=trim(left(put(month,z2.)))||'/15/'||substr(trim(left(year)),3,2);
format date date9.;
date=input(datestring,mmddyy8.);

/* Get the month number for the cycle graph */
monthnum=0; monthnum=put(date,year4.);

run;

proc sql;
select year into :year separated by ' ' from co2data having decimal_date=max(decimal_date);
select month into :month separated by ' ' from co2data having decimal_date=max(decimal_date);
quit; run;

 
/* First, create the average cycle plot, to be used as the inset graph */
proc sql;
create table cycle_data as 
select unique month, avg(co2) as co2
from co2data 
group by month;
quit; run;
data cycle_data; set cycle_data;
output;
if (month eq 1) then do;
 month=13;
 output;
 end;
run;
proc sort data=cycle_data out=cycle_data;
by month;
run;

/* Create user-defined format to print numeric month as month abbreviation */
proc format;
value monfmt
1='Jan'
2='Feb'
3='Mar'
4='Apr'
5='May'
6='Jun'
7='Jul'
8='Aug'
9='Sep'
10='Oct'
11='Nov'
12='Dec'
13='Jan'
;
run;





/* Create the small inset graph to show the annual cycle */

goptions device=png;
goptions xpixels=250 ypixels=190;
goptions cback=cxFFFFCC;

ODS LISTING CLOSE;
ODS HTML path=odsout body="cycle.htm" (title="Small inset 'cycle' graph") style=htmlblue;

goptions gunit=pct htitle=8 htext=8 ftitle="albany amt/bold" ftext="albany amt/bold";

symbol1 i=join width=1.5 value=dot height=5.5 mode=include;

axis1 label=none major=none minor=none value=none offset=(2,2);
axis2 label=none major=(height=2pct) minor=none order=(1 to 13 by 3) offset=(0,0);

title1 ls=1.5 "Averaged Annual Cycle";
title2 a=90 h=5 " ";
title3 a=-90 h=5 " ";
symbol1 value=dot interpol=join color=dodgerblue mode=include;
proc gplot data=cycle_data;
format month monfmt.;
plot co2*month / 
 vaxis=axis1 haxis=axis2
 cframe=white 
 des='' name="cycle";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;






data my_anno;
length function $8 style $12 text $100;
xsys='1'; ysys='1'; hsys='3'; when='a';

/* Annotate the 'annual cycle' image in bottom/right */
function='move'; 
x=67.5; y=3; output;
function='image'; imgpath='cycle.png'; style='fit';
x=x+31; y=y+41; output;

/* Annotate a yellow box, with title text, in upper/left */
html='href='||quote("&urlname");
function='move'; 
x=1; y=75; output;
function='bar'; color='cxFFFFCC'; style='solid'; line=0;
x=76; y=99; output;
html='';
function='move'; 
x=1; y=75; output;
function='bar'; color='grayaa'; style='empty'; line=0;
x=76; y=99; output;

function='label'; position='5'; color='';
style=''; 
x=38.5; y=94; size=8; text="Atmospheric Carbon Dioxide"; output;
x=38.5; y=85; size=6; text="Measured at Mauna Loa, Hawaii"; output;
run;


/* 12-month moving average, with 6 months trimmed off the end */
/* You'll need SAS/ETS for Proc Expand */
proc expand DATA=co2data OUT=co2data;
convert co2 = co2_12 / METHOD = none TRANSFORMOUT = (cmovave 12 trim 6);
run;


/* Now, create the actual/big graph */

goptions device=png;
goptions xpixels=900 ypixels=550;
goptions noborder;
goptions cback=white;  /* over-ride the yellow cback I declared earliers for the inset graph */

ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Atmospheric Carbon Dioxide (at Mauna Loa, Hawaii)") style=htmlblue;

goptions gunit=pct htitle=5.5 htext=3.25 ftitle="albany amt" ftext="albany amt";
goptions ctext=gray33;

axis1 label=(angle=90) order=(300 to 440 by 20) minor=none offset=(0,0);

axis2 label=none order=('01jan1960'd to '01jan2030'd by year10) minor=(number=1) offset=(0,0);

symbol1 interpol=spline color=red width=.1 value=none mode=include;
symbol2 interpol=sm12 color=blue width=.1 value=none mode=include;

/* This is the 'invisible' line for the plot1 statement.
   You must have a plot1, in order to have a plot2 (we want
   plot2 for the axis on the right-hand side.
*/
symbol3 i=none value=none mode=include;

/* Here's the 'invisible' axis on the left-hand side 
   I will use with plot1 */
axis3 label=none value=none major=none minor=none;


footnote1 
 link="&urlname"
 j=r move=(-7,+0) font="albany amt" c=gray "As of year=&year month=&month";

/* Plot the red data line, and then overlay the blue smoothed line */
title h=2 " ";
proc gplot data=co2data anno=my_anno;
format date year.;
plot co2*date=3 /    /* invisible/fake plot */
 vaxis=axis3 haxis=axis2
 autohref chref=gray77 lhref=33
 des='' name="&name";
plot2 co2*date=1 co2_12*date=2 / overlay
 vaxis=axis1 
 autovref cvref=gray77 lvref=33
 ;
run;

title1 ls=1.5 "Same data, plotted to a zero-scale...";

axis1 order=(0 to 450 by 50) label=(angle=90) minor=none offset=(0,0);
axis2 label=none order=('01jan1960'd to '01jan2030'd by year10) minor=(number=1) offset=(0,0);
axis3 label=none value=none major=none minor=none;

proc gplot data=co2data;
format date year.;
plot co2*date=3 / 
 vaxis=axis1 haxis=axis2 
 autohref chref=graydd lhref=1
 autovref cvref=graydd lvref=1
 des='' name="&name";
plot2 co2*date=1 / 
 vaxis=axis1;
run;

quit;
ODS HTML CLOSE;
ODS LISTING;