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

/*
The data for this one has moved around over time...

http://www.globalwarmingart.com/wiki/Image:Sunspot_Numbers_png
ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/MONTHLY.PLT
ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/GROUP_SUNSPOT_NUMBERS/monthrg.dat
http://www.ngdc.noaa.gov/stp/SOLAR/solarda3.html
filename dataurl url "ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/MONTHLY.PLT"
filename dataurl url "ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/INTERNATIONAL/monthly/MONTHLY.PLT"
http://www.sidc.be/silso/DATA/SN_m_tot_V2.0.txt
*/

%let datafile=https://wwwbis.sidc.be/silso/DATA/SN_m_tot_V2.0.txt;
filename dataurl url "&datafile";

data newer_data;
infile dataurl dlm=' ';
input year month yymm sunspot foo1 foo2;
date=input('15-'||trim(left(put(month,z2.)))||'-'||trim(left(year)),ddmmyy10.);
run;


/*
Justin Mabie <justin.mabie@noaa.gov>
says they no longer plan to keep this old-old data on the web,
so I'm saving a snapshot...
filename data2url url "ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/GROUP_SUNSPOT_NUMBERS/monthrg.dat"
 proxy="&my_proxy";
*/
filename data2url "monthrg.dat";
data older_data;
infile data2url dlm=' ';
input year month foo1 sunspot foo2;
date=input(trim(left(put(month,z2.)))||'/15/'||trim(left(year)),mmddyy10.);
if (sunspot eq -99.0) then sunspot=.;
if (year<1749) then output;
run;


data sundata; set older_data newer_data;
label sunspot='Number of Sunspots';
if date<='01aug1748'd then colorvar=2;
else colorvar=1;
run;

proc sql noprint;
select unique max(date) format=monyy8. into :maxdate separated by ' ' from sundata;
select unique max(date) into :max_x_num separated by ' ' from sundata;
quit; run;

 

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

function='label'; 
angle=0;
x='01jan1680'd; y=95; 
position='2'; text='Maunder'; output;
position='5'; text='Minimum'; output;

x='01jan1971'd; y=345; 
position='2'; text='Modern'; output;
position='5'; text='Maximum'; output;

angle=90;
x='01jan1815'd; y=250; 
position='2'; text='Dalton'; output;
position='5'; text='Minimum'; output;

xsys='2'; ysys='3';
x=&max_x_num; y=3;
position='4'; angle=.; rotate=.; style=''; color="gray77"; size=3.5; style="albany amt";
text="Data through: &maxdate"; 
html='href='||quote("&datafile");
output;

run;



goptions device=png;
goptions xpixels=1000 ypixels=400;
goptions noborder;
goptions cback=cxe0e0e0;

ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm"
 (title="400 Years of Sunspot Observations") style=htmlblue;

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

axis1 label=none value=none 
 minor=none offset=(0,0) order=(0 to 400 by 100);

axis2 label=(angle=90 'Sunspots per Month')
 minor=none offset=(0,0) order=(0 to 400 by 100) value=(t=5 ' ');

axis3 
 order=('01jan1600'd, '01jan1650'd, '01jan1700'd, '01jan1750'd, '01jan1800'd, 
        '01jan1850'd, '01jan1900'd, '01jan1950'd, '01jan2000'd, '01jan2050'd) 
 label=none /* minor=none */ minor=(number=4) offset=(0,0);

/* There's probably some extra leap-year stuff in the 1600s and 1700s that's 
   throwing off my days-per-year, when I try to do it this way.
   365.25 x 50 = 18262.5 days per 50-year
 order=('01jan1600'd to '10jan2050'd by 18262.5) 
*/

symbol1 interpol=join mode=include color=cx436EEE width=.01 value=none;
symbol2 interpol=none color=red height=2.0 font="albany amt" value='X';
symbol3 interpol=sm12 color=black width=2 value=none;

/* Compute the moving average (black line) using a 
   240-month (20-year) window, and trim off 120-month (10 years) 
   at each end (otherwise the ends of the black line will be unduly
   influenced by the data at the ends of the graph
 */
/* You'll need SAS/ETS for Proc Expand */
proc expand DATA=sundata OUT=sundata;
convert sunspot = sunspot_mov_avg / METHOD = none TRANSFORMOUT = (cmovave 240 trim 120);
run;

footnote " ";

title1 link="&datafile" ls=1.5 "400 Years of Sunspot Observations";

proc gplot data=sundata anno=my_anno;
format date year4.;
plot sunspot*date=colorvar / nolegend
 vaxis=axis1 haxis=axis3
 autovref cvref=gray77 lvref=2 
 cframe=white
 des='' name="&name";
plot2 sunspot_mov_avg*date=3 / 
 vaxis=axis2;
run;



symbol4 interpol=join mode=include color=cx436EEE width=.01 value=circle height=2.0;

axis4 
 order=(
 '01jan2009'd, 
 '01jan2010'd, 
 '01jan2011'd, 
 '01jan2012'd, 
 '01jan2013'd, 
 '01jan2014'd, 
 '01jan2015'd, 
 '01jan2016'd, 
 '01jan2017'd, 
 '01jan2018'd, 
 '01jan2019'd, 
 '01jan2020'd
 '01jan2021'd
 '01jan2022'd
/*
*/
 ) 
 label=none minor=none offset=(0,0);

title1 link="&datafile" ls=1.5 "Current Sunspot Cycle";

/* don't need to subset here - subsetting will be done by axis statement */
data plot2; set sundata /*(where=(date>='01jan2010'd))*/;
length my_html $500;
my_html=
 'title='||quote(trim(left(put(sunspot,comma8.1)))||' = mean sunspots during month of '||trim(left(put(date,monyy7.))))||
 ' href='||quote('sunspot_info.htm');
run;

axis1 label=none 
 minor=none offset=(0,0) order=(0 to 400 by 100);

proc gplot data=plot2;
format date year4.;
plot sunspot*date=4 / nolegend
 vaxis=axis1 haxis=axis4
 autovref cvref=gray77 lvref=2 
 autohref chref=graydd
 cframe=white
 html=my_html
 des='' name="&name";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;