%let name=la_ozone_calendar_chart; filename odsout "."; /* Data from: https://aqs.epa.gov/aqsweb/airdata/download_files.html#Raw AQI for Ozone, color ranges: http://www.nws.noaa.gov/ost/air_quality/FAQ_01172011.pdf */ libname here '.'; /* Just run this part once, and save the resulting dataset */ /* %macro readdata(year); PROC IMPORT OUT=ozone_&year DATAFILE="D:\Public\EPA\daily_44201_&year..csv" DBMS=CSV REPLACE; GETNAMES=YES; DATAROW=2; RUN; data ozone_&year; set ozone_&year (where=(state_name='Californi' and local_site_name='Pasadena')); day=.; day=input(date_local,yymmdd10.); monname=trim(left(put(day,monname.))); year=put(day,year.); ozone_ppm=arithmetic_mean; run; %mend; %readdata(2016); %readdata(2015); %readdata(2014); %readdata(2013); %readdata(2012); %readdata(2011); %readdata(2010); %readdata(2009); %readdata(2008); %readdata(2007); %readdata(2006); %readdata(2005); %readdata(2004); %readdata(2003); %readdata(2002); %readdata(2001); %readdata(2000); %readdata(1999); %readdata(1998); %readdata(1997); %readdata(1996); %readdata(1995); %readdata(1994); %readdata(1993); %readdata(1992); %readdata(1991); %readdata(1990); %readdata(1989); %readdata(1988); %readdata(1987); %readdata(1986); %readdata(1985); %readdata(1984); %readdata(1983); %readdata(1982); %readdata(1981); data la_ozone; set ozone_2016 ozone_2015 ozone_2014 ozone_2013 ozone_2012 ozone_2011 ozone_2010 ozone_2009 ozone_2008 ozone_2007 ozone_2006 ozone_2005 ozone_2004 ozone_2003 ozone_2002 ozone_2001 ozone_2000 ozone_1999 ozone_1998 ozone_1997 ozone_1996 ozone_1995 ozone_1994 ozone_1993 ozone_1992 ozone_1991 ozone_1990 ozone_1989 ozone_1988 ozone_1987 ozone_1986 ozone_1985 ozone_1984 ozone_1983 ozone_1982 ozone_1981 ; run; data here.la_ozone; set la_ozone; run; */ data la_ozone; set here.la_ozone; run; /* from http://www.nws.noaa.gov/ost/air_quality/FAQ_01172011.pdf */ proc format; value aqi_fmt 1='Good (0-59 ppb)' 2='Moderate (60-75 ppb)' 3='Unhealthy for Sensitive Groups (76-95 ppb)' 4='Unhealthy (96-115 ppb)' 5='Very Unhealthy (116-374 ppb)' ; run; data la_ozone; set la_ozone; format colorbin aqi_fmt.; ozone_ppb=round(_1st_Max_Value*1000); if ozone_ppb=. then colorbin=.; else if ozone_ppb<=59 then colorbin=1; else if ozone_ppb<=75 then colorbin=2; else if ozone_ppb<=95 then colorbin=3; else if ozone_ppb<=115 then colorbin=4; else colorbin=5; run; /* First, make sure the data is sorted, so you can use "by year" later */ proc sort data=la_ozone out=la_ozone; by year day; run; proc sql noprint; /* Macro variables containing minimum & maximum years */ select min(year) into :min_year from la_ozone; select max(year) into :max_year from la_ozone; quit; run; /* My algorithm assumes you have an obs for each day, so create a grid of all days */ data grid_days; format day date9.; do day="01jan.&min_year"d to "31dec.&max_year"d by 1; weekday=put(day,weekday.); downame=trim(left(put(day,downame.))); monname=trim(left(put(day,monname.))); year=put(day,year.); output; end; run; /* Join your data with the grid-of-days */ proc sql noprint; create table la_ozone as select * from grid_days left join la_ozone on grid_days.day eq la_ozone.day; quit; run; /* Add some html title= flyover data tip text (could also add href= drilldown here) */ /* The "0d"x is a carriage-return, which is honored by most web browsers */ data la_ozone; set la_ozone; length myhtmlvar $200; my_html= 'title='||quote( trim(left(put(day,downame.)))||'0D'x|| put(day,date9.)||'0D'x|| 'Ozone: '||put(ozone_ppb,comma5.4)||' ppm'); run; /* Create a "map" of the calendar days, suitable for use in gmap (ie, 4 X/Y obsns per each day) */ /* When using "by year", you can then use "first.year" */ /* You're starting with minimum date at top/left, max at bottom/right */ data datemap; set la_ozone; keep day x y; by year; if first.year then x_corner=1; else if trim(left(downame)) eq "Sunday" then x_corner+1; y_corner=((&min_year-year)*8.5)-weekday; /* output 4 X/Y coordinates per each day, forming a rectangle that GMAP can draw */ x=x_corner; y=y_corner; output; x=x+1; output; y=y-1; output; x=x-1; output; run; /* Create darker outline to annotate around each month */ /* (this is similar to annotating a state outline onto a county map) */ data anno_month_outline; set datemap; /* combination of year & month makes a unique id for these outlines */ length yr_mon $ 15; yr_mon=trim(left(put(day,year.)))||"_"||trim(left(put(day,month.))); order+1; run; /* Sort it, so you can use "by" in next step */ proc sort data=anno_month_outline out=anno_month_outline; by yr_mon order; run; /* Remove the internal borders, within each month */ proc gremove data=anno_month_outline out=anno_month_outline; by yr_mon; id day; run; /* Now, convert the gmap data set into annotate move/draw commands */ data anno_month_outline; length color function $8; xsys="2"; ysys="2"; size=1.75; when="A"; color="black"; set anno_month_outline; by yr_mon; if first.yr_mon then function='poly'; else function='polycont'; run; /* Annotate some text labels for year and day-of-week, along the left */ proc sql noprint; create table anno_year_and_weekday as select unique year from la_ozone; quit; run; data anno_year_and_weekday; set anno_year_and_weekday; length text $10; function="label"; position="4"; color='black'; xsys="2"; ysys="2"; hsys="d"; when="A"; x=-8; y=((&min_year-year)*8.5)-1.25; style=""; size=12; text=trim(left(year)); output; x=-.1; size=7; color=''; text="Sunday"; output; y=y-1; text="Monday"; output; y=y-1; text="Tuesday"; output; y=y-1; text="Wednesday"; output; y=y-1; text="Thursday"; output; y=y-1; text="Friday"; output; y=y-1; text="Saturday"; output; run; /* Annotate some labels for the 3-character month name, along the top */ data anno_month; length text $10; function="label"; position="5"; color='black'; xsys="2"; ysys="2"; hsys="d"; when="A"; size=10; y=1; spacing=4.5; x=(spacing/3)*-1; x=x+spacing; text="JAN"; output; x=x+spacing; text="FEB"; output; x=x+spacing; text="MAR"; output; x=x+spacing; text="APR"; output; x=x+spacing; text="MAY"; output; x=x+spacing; text="JUN"; output; x=x+spacing; text="JUL"; output; x=x+spacing; text="AUG"; output; x=x+spacing; text="SEP"; output; x=x+spacing; text="OCT"; output; x=x+spacing; text="NOV"; output; x=x+spacing; text="DEC"; output; run; /* Combine the weekday & month annotated text into 1 data set. */ data text_anno; set anno_year_and_weekday anno_month; run; /* Put a fake map area to the top/left of the map, to guarantee room for the annotated labels on the left & top */ /* You might need to adjust the "-10" value to add more space, depending on your text size, map proportions, etc. */ data fake; day=1; x=-10; y=1; output; x=x-.001; y=y+.001; output; x=x+.002; output; run; data datemap; set datemap fake; run; goptions xpixels=700 ypixels=3000; goptions device=png; goptions border; ODS LISTING CLOSE; ODS HTML path=odsout body="&name..htm" (title="Ozone Levels in Los Angeles") style=htmlblue; goptions gunit=pct ftitle="albany amt/bold" ftext="albany amt" htitle=16pt htext=11pt; goptions ctext=gray33; pattern1 v=s c=cx00e400; pattern2 v=s c=cxffff00; pattern3 v=s c=cxefa011; pattern4 v=s c=cxff0000; pattern5 v=s c=cx99004c; legend1 shape=bar(.15in,.15in) label=none value=(justify=left) across=1 order=descending; title1 ls=1.5 "Los Angeles, CA - Ozone Levels"; title2 "Daily Maximum Values"; proc gmap data=la_ozone map=datemap all anno=text_anno; id day; choro colorbin / discrete midpoints = 1 2 3 4 5 cdefault=graydd cempty=graycc coutline=graycc anno=anno_month_outline legend=legend1 html=my_html des='' name="&name"; run; quit; ODS HTML CLOSE; ODS LISTING;