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

/*
Imitation/simplification of graph & maps from:
https://www.washingtonpost.com/graphics/national/power-plants/
Using data from:
https://www.eia.gov/electricity/data/state/
(Net Generation by State by Type of Producer by Energy)
*/

proc import
 file="D:\Public\EIA\Power\annual_generation_state.xls"
 out=my_data dbms=xls replace;
getnames=yes;
range='$A2:E0';
run;

proc sql noprint;
select max(year) into :maxyear separated by ' ' from my_data;
create table my_data as
select unique state as statecode, energy_source, generation__megawatthours_
from my_data
having type_of_producer='Total Electric Power Industry' and year=max(year)
 and statecode not in ('DC' 'US-Total');
quit; run;

proc transpose data=my_data out=my_data (drop = _name_ _label_);
by statecode;
id energy_source;
run;

data my_data; set my_data (rename=(
 hydroelectric_conventional=hydro 
 solar_thermal_and_photovoltaic=solar))
;
if coal=. then coal=0;
if natural_gas=. then natural_gas=0;
if nuclear=. then nuclear=0;
if hydro=. then hydro=0;
if wind=. then wind=0;
if solar=. then solar=0;
if petroleum=. then petroleum=0;
coal_pct=coal/total;
natural_gas_pct=natural_gas/total;
nuclear_pct=nuclear/total;
hydro_pct=hydro/total;
wind_pct=wind/total;
solar_pct=solar/total;
petroleum_pct=petroleum/total;
other=total-(coal+natural_gas+nuclear+hydro+wind+solar+petroleum);
other_pct=other/total;
run;


data my_map; set mapsgfk.us;
original_order=_n_;
run;

/* Get the state centroids to position the annotated pies */
%annomac;
proc sort data=my_map out=my_map;
by statecode original_order;
run;
/* get the x/y center of each state */
%centroid(my_map, pie_centers, statecode);





options mprint;

%macro do_plots(varname,source,rgb);
ods html anchor="&source";

proc sort data=my_data out=bar_temp;
by descending &varname statecode;
run;

data bar_temp; set bar_temp;
bar_order=_n_;
run;

proc sql noprint;
create table control as 
select unique bar_order as start, statecode as label 
from bar_temp;
quit; run;
data control; set control;
fmtname='barfmt';
type='N';
run;
proc format lib=work cntlin=control;
run;

data bar_temp; set bar_temp;
length category $50;
format total_value comma20.0;
format pct_value percentn7.2;
label category='Source';
label total_value='MWh';
label pct_value='Percent';
category='Coal';        total_value=coal;        pct_value=coal_pct; output;
category='Natural gas'; total_value=natural_gas; pct_value=natural_gas_pct; output;
category='Nuclear';     total_value=nuclear;     pct_value=nuclear_pct; output;
category='Hydro';       total_value=hydro;       pct_value=hydro_pct; output;
category='Wind';        total_value=wind;        pct_value=wind_pct; output;
category='Solar';       total_value=solar;       pct_value=solar_pct; output;
category='Petroleum';   total_value=petroleum;   pct_value=petroleum_pct; output;
category='Other';       total_value=other;       pct_value=other_pct; output;
run;

data bar_temp; set bar_temp;
length my_html my_html_legend $300;
my_html=
 'title='||quote(
   trim(left(fipnamel(stfips(statecode))))||': '||'0d'x||
   trim(left(put(total_value,comma10.0)))||' MWh of electric power comes from '||trim(left(category))||'.'||'0d'x||
   "That's "||trim(left(put(pct_value,percent7.0)))||" of this state's electric power."
  )||
 ' href='||quote('#'||translate(trim(left(category)),'_',' '));
my_html_legend=
 'title='||quote(trim(left(category)))||
 ' href='||quote('#'||translate(trim(left(category)),'_',' '));

/* make the non-selected ones show up below the axis (negative values) */
if category^=translate("&source",' ','_') then pct_value=-1*pct_value;
run;

%let outclr=graybb;

/* 
These colors are assigned by the default alphabetical legend order
(and then I manually change the order of the legend in the legend1 statement)
*/
pattern1 v=s c=cx555555; /* Coal */
pattern2 v=s c=cx0081c5; /* Hydro */
pattern3 v=s c=cxf68b28; /* Natural gas */
pattern4 v=s c=cxcf4a9a; /* Nuclear */
pattern5 v=s c=cxa6761d; /* Other */
pattern6 v=s c=cxed1c24; /* Petroleum */
pattern7 v=s c=cxd7c944; /* Solar */
pattern8 v=s c=cx0db14b; /* Wind */

axis1 label=none order=(-1.00 to 1.00 by .20) 
 major=none minor=none style=0 offset=(0,0)
 value=(
 t=11 "100%" 
 t=10 "80" 
 t=9 "60" 
 t=8 "40" 
 t=7 "20" 
 t=6 "0" 
 t=5 " " 
 t=4 " " 
 t=3 " " 
 t=2 " " 
 t=1 " "
 );

axis2 label=none style=0 /*value=(height=.1 color=white)*/ value=none offset=(1.5,1.5)pct;

legend1 label=none position=(top center inside) shape=bar(.18in,.18in) 
 order=('Coal' 'Natural gas' 'Nuclear' 'Hydro' 'Wind' 'Solar' 'Petroleum' 'Other')
 value=(font='albany amt/bold' height=2.9) offset=(10,-5);
/* foofoo */

/*
Gchart doesn't think it can fit the state abbreviations without 
automatically 'stacking' the characters one over the other.
Therefore I have to manually annotate them.
*/
data anno_statecodes; set bar_temp (where=(category=translate("&source",' ','_')));
if pct_value^=0 then do;
 xsys='2'; ysys='2'; hsys='3'; when='a';
 x=bar_order; subgroup=category; /*y=0;*/
 function='label'; position='2'; size=1.7; text=trim(left(statecode));
 end;
else do;
 xsys='2'; ysys='2'; hsys='3'; when='a';
 x=bar_order; y=0;
 function='label'; position='2'; size=1.7; text=trim(left(statecode));
 end;
run;

data anno_axis;
do y=0 to -1 by -.2;
 output;
 end;
run;
data anno_axis; set anno_axis;
xsys='1'; ysys='2'; hsys='3'; when='a';
x=100.5;
function='label'; position='>';
if y=-1 then text=trim(left(put(abs(y),percentn7.0)));
else text=trim(left(put(abs(y)*100,comma8.0)));
run;



goptions xpixels=1100;
title1 ls=1.5 "US Electricity Generation by Power Source";
title2 a=-90 h=4pct ' ';
footnote 
 link='https://www.eia.gov/electricity/data/state/'
 j=l move=(+6,+0) c=gray33 h=2.2
 "Data source: eia.gov/electricity/data/state/  (&maxyear)";
proc gchart data=bar_temp anno=anno_statecodes;
format pct_value percent7.0;
format bar_order barfmt.;
vbar bar_order / discrete 
 width=2.0 space=0 coutline=&outclr
 type=sum sumvar=pct_value
 subgroup=category
 raxis=axis1 maxis=axis2 noframe
 anno=anno_axis /* annotate the values along right axis */
 autoref lref=33 cref=gray55 clipref
 legend=legend1
 html=my_html
 html_legend=my_html_legend
 des='' name="&name._bar_&source";
run;

footnote;
title2;

/* maximum circle area */
%let max_area=26; 
/* maximum possible data value, to scale all circles to */
%let max_val=1.00;

proc sql noprint;
create table anno_pies as
select unique statecode, category, pct_value, total_value
from bar_temp
where category=translate("&source",' ','_');
/*
where category="&source";
*/
create table anno_pies as
select unique pie_centers.*, anno_pies.pct_value, anno_pies.total_value
from pie_centers left join anno_pies
on pie_centers.statecode=anno_pies.statecode;
quit; run;

/* hmm ... since the pies are all the same size, this doesn't reall matter/help */
proc sort data=anno_pies out=anno_temp;
by descending pct_value;
run;

data anno_temp; set anno_temp (where=(pct_value^=.));
length function $8 color $15;
xsys='2'; ysys='2'; hsys='3'; when='a';

function='pie'; 
size=.5 + sqrt((&max_val/&max_val)*&max_area/3.14);
rotate=-1*pct_value*360; style='psolid'; color="cx&rgb"; angle=90; 
if pct_value^=0 then output;
 
length html $300;
html=
 'title='||quote(
   trim(left(fipnamel(stfips(statecode))))||': '||'0d'x||
   trim(left(put(total_value,comma10.0)))||' MWh of electric power comes from '||trim(left("&source"))||'.'||'0d'x||
   "That's "||trim(left(put(pct_value,percent7.0)))||" of this state's electric power."
  )||
 ' href='||quote('https://www.eia.gov/state/analysis.php?sid='||trim(left(upcase(statecode))));
/*
 ' href='||quote('us_power_sources_info.htm');
*/

rotate=360; style='pempty'; color="&outclr"; output;
run;

/*
Create the fake data this way, using a user-defined-format for the text,
so that the colors will always be assigned in the same/predictable order,
rather than alphabetical by the text.
*/
data fake_data;
length legend_text $50;
statecode='AA'; choro_value=1; legend_text=translate("&source",' ','_');; output;
statecode='BB'; choro_value=2; legend_text="Other than "||trim(left(translate("&source",' ','_'))); output;
run;
proc sql noprint;
create table control as 
select unique choro_value as start, legend_text as label 
from fake_data;
quit; run;
data control; set control;
fmtname='myfmt';
type='N';
run;
proc format lib=work cntlin=control;
run;

pattern1 v=s c=cx&rgb;
pattern2 v=s c=white;
legend2 label=none position=(top) shape=bar(.18in,.18in) value=(height=3.0)
 mode=share across=2;

goptions xpixels=1000;
proc gmap data=fake_data map=my_map anno=anno_temp all;
/* Use note rather than footnote, to share the map's space */
note 
 link='https://www.eia.gov/electricity/data/state/'
 move=(49.5,2) c=gray33 h=2.2
 "Data source: eia.gov/electricity/data/state/  (&maxyear)";
id statecode;
format choro_value myfmt.;
choro choro_value / 
 cdefault=white 
 coutline=&outclr
 cempty=gray99 
 legend=legend2
 des='' name="&name._map_&source";
run;

/*
proc print data=anno_axis; run;
proc print data=bar_temp; run;
proc print data=anno_temp; run;
*/

%mend;


ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="US Power Sources") 
 style=htmlblue;

goptions device=png noborder;
goptions ypixels=600;

goptions gunit=pct ftitle="albany amt/bold" ftext="albany amt";
goptions htitle=4.8 htext=2.0 ctext=gray22 ctitle=gray22;

/*
*/
%do_plots(coal_pct,Coal,555555);
%do_plots(hydro_pct,Hydro,0081c5);
%do_plots(natural_gas_pct,Natural_gas,f68b28);
%do_plots(nuclear_pct,Nuclear,cf4a9a);
%do_plots(petroleum_pct,Petroleum,ed1c24);
%do_plots(other_pct,Other,a6761d);
%do_plots(solar_pct,Solar,d7c944);
%do_plots(wind_pct,Wind,0db14b);

quit;
ODS HTML CLOSE;
ODS LISTING;