%let name=albany_ny_loess; filename odsout '.'; /* This is a slight variation of the SAS analysis/graph Rick Wicklin wrote, pertaining to the following "SAS and R" blog: http://sas-and-r.blogspot.com/2012/04/example-925-its-been-mighty-warm-winter.html */ %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 webfile url "http://academic.udayton.edu/kissock/http/Weather/gsod95-current/NYALBANY.txt" proxy="&my_proxy"; data TempData; infile webfile; input month day year temperature; format date date9.; date=MDY(month, day, year); dayofyear=0; dayofyear=put(date,julday3.); label Proportion = "Day of Year (Normalized)"; Proportion = dayofyear / put(mdy(12, 31, year), julday3.); /* incorporate leap years */ if Temperature^=-99 then output; run; /* extend data by translating left and right */ data Periodic; set TempData(in=before) TempData TempData(in=after); if before then proportion = proportion - 1; /* (-1,0] */ if after then proportion = proportion + 1; /* (1,2] */ run; data Score; do proportion = 0 to 1 by 1/365; output; end; /* 3 times the data, so use 1/3 the smoothing parameter! */ proc loess data=Periodic plots(only maxpoints=none)=(FitPlot CriterionPlot); model Temperature = Proportion/ smooth=0.0557 interp=cubic; score data=Score; ods output ScoreResults=Fit; run; /* Estimate the date for each p_temperature */ data Fit; set Fit; dayofyear=proportion*365; run; quit; run; /* Combine back with the main data, so you can plot them together */ data combined; set TempData Fit (keep=proportion dayofyear p_temperature); run; /* Create an extra variable for the current winter's temperature */ data combined; set combined; length my_html $200; if (date >='01dec2011'd & date<='15mar2012'd) then do; CurrentWinterTemperature=temperature; my_html= 'title='||quote(put(date,date9.)||': '||trim(left(temperature)))|| ' href='||quote('http://academic.udayton.edu/kissock/http/Weather/gsod95-current/NYALBANY.txt');; end; run; goptions device=png; ODS LISTING CLOSE; ODS HTML path=odsout body="&name..htm" (title="Albany NY Temperature") style=htmlblue; goptions htitle=13pt htext=9pt ftitle="albany amt/bold" ftext="albany amt"; goptions ctext=gray33; axis1 label=(angle=90 font="albnay amt/bold" 'Temperature') minor=none offset=(0,0); axis2 label=none order=('01jan1960'd to '01jan1961'd by month) minor=none value=(t=13 ''); symbol1 value=circle height=.7 interpol=none color=A2c56ab33; symbol2 value=dot height=1.0 interpol=none color=gray55; symbol3 value=none interpol=join width=2 color=cx002a94; legend1 label=none position=(bottom center) cborder=gray55; proc sql noprint; select unique min(year) into :minyear separated by ' ' from TempData; select unique max(year) into :maxyear separated by ' ' from TempData; quit; run; title ls=1.5 link="http://academic.udayton.edu/kissock/http/Weather/gsod95-current/NYALBANY.txt" "Temperature in Albany, NY (&minyear-&maxyear)"; proc gplot data=combined; format dayofyear monname3.; label temperature='Temperature'; label CurrentWinterTemperature='Winter 2011-2012'; label p_temperature='Periodic smoother'; plot temperature*dayofyear=1 CurrentWinterTemperature*dayofyear=2 p_temperature*dayofyear=3 / overlay legend=legend1 autovref cvref=graydd autohref chref=graydd vaxis=axis1 haxis=axis2 html=my_html des='' name="&name"; run; quit; ODS HTML CLOSE; ODS LISTING;