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

/*
Creating a map similar to the New York Times map David Mintz
had on p. 6 of his paper: 
"Five Crazy Good Visualizationa and How to Plot Them"
https://www.lexjansen.com/sesug/2018/SESUG2018_Paper-214_Final_PDF.pdf

Using data from:
https://www.ers.usda.gov/data-products/county-level-data-sets/
https://www.ers.usda.gov/data-products/county-level-data-sets/download-data/
*/

/* Poverty */

proc import out=poverty datafile="PovertyEstimates.xls" dbms=xls replace;
range='Poverty Data 2017$A4:K0';
getnames=yes;
run;

data poverty (rename=(FIPStxt=FIPS_Code PCTPOVALL_2017=Poverty Area_Name=county_name)); 
 set poverty (keep = FIPStxt State Area_Name PCTPOVALL_2017 rename=(state=statecode));
county=.;
county=substr(FIPStxt,3,3);
format PCTPOVALL_2017 percent7.1;
PCTPOVALL_2017=PCTPOVALL_2017/100;
if county^=0 and PCTPOVALL_2017^=. then output;
run;

proc sort data=poverty out=poverty;
by Poverty;
run;
data poverty; set poverty;
poverty_rank=_n_;
run;


/* Education */
/* Percent of adults with a bachelor's degree or higher, 2013-17 */

proc import out=education datafile="Education.xls" dbms=xls replace;
range='Education 1970 to 2017$A5:AU0';
getnames=yes;
run;

data education (rename=(Percent_of_adults_with_a_bachel2=Education Area_Name=county_name) drop=FIPS_Code); 
 set education (keep = FIPS_Code State Area_Name Percent_of_adults_with_a_bachel2 rename=(state=statecode));
county=.; 
county=substr(FIPS_Code,3,3);
format Percent_of_adults_with_a_bachel2 percent7.1;
Percent_of_adults_with_a_bachel2=Percent_of_adults_with_a_bachel2/100;
if county^=0 and Percent_of_adults_with_a_bachel2^=. then output;
run;

proc sort data=education out=education;
by descending education;
run;
data education; set education;
education_rank=_n_;
run;


/* Unemployment */

proc import out=unemployment datafile="Unemployment.xls" dbms=xls replace;
range='Unemployment Med HH Inc$A8:BB0';
getnames=yes;
run;

data unemployment (rename=(Unemployment_rate_2018=Unemployment Area_Name=county_name) drop=FIPS); 
 set unemployment (keep = FIPS State Area_Name Unemployment_rate_2018 rename=(state=statecode));
county=.; 
county=substr(trim(left(put(FIPS,z5.))),3,3);
format Unemployment_rate_2018 percent7.1;
Unemployment_rate_2018=Unemployment_rate_2018/100;
if county^=0 and Unemployment_rate_2018^=. then output;
run;

proc sort data=unemployment out=unemployment;
by unemployment;
run;
data unemployment; set unemployment;
unemployment_rank=_n_;
run;


/* Combine */

proc sql noprint;
create table my_data as
select unique poverty.*, education.education, education.education_rank
from poverty left join education
on poverty.statecode=education.statecode and poverty.county=education.county;

create table my_data as
select unique my_data.*, unemployment.unemployment, unemployment.unemployment_rank
from my_data left join unemployment
on my_data.statecode=unemployment.statecode and my_data.county=unemployment.county;

quit; run;

data my_data; set my_data;
if poverty_rank^=. and education_rank^=. and unemployment_rank^=.
 then combined=(poverty_rank+education_rank+unemployment_rank);
run;

proc sort data=my_data out=my_data;
by combined;
run;
data my_data; set my_data;
overall_rank=_n_;
length my_html $300;
my_html='title='||quote(
 trim(left(county_name))||', '||trim(left(statecode))||' '||'(overall rank #'||trim(left(put(overall_rank,comma8.0)))||')'||'0d'x||
 '-------------------'||'0d'x||
 'Poverty: '||put(poverty,percent7.1)||' (#'||trim(left(put(poverty_rank,comma8.0)))||')'||'0d'x||
 'College degree: '||put(education,percent7.1)||' (#'||trim(left(put(education_rank,comma8.0)))||')'||'0d'x||
 'Unemployment: '||put(unemployment,percent7.1)||' (#'||trim(left(put(unemployment_rank,comma8.0)))||')'
 );
run;


libname robsmaps '../democd97';
data my_map; set robsmaps.uscounty (where=(density<=1));
run;

proc gremove data=my_map out=anno_outline;
by state; 
id county;
run;
data anno_outline; set anno_outline; 
by state segment notsorted;
length function $8 color $8;
color='white'; style='mempty'; when='a'; xsys='2'; ysys='2';
if first.segment then function='poly';
else function='polycont';
run;


goptions device=png border;
goptions xpixels=1100 ypixels=800;

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

goptions gunit=pct htitle=28pt htext=12pt ftitle="albany amt/bold" ftext="albany amt";
goptions ctext=gray33;

pattern1 v=s c=cx397b7e;
pattern2 v=s c=cx77a5a8;
pattern3 v=s c=cxa7c5c6;
pattern4 v=s c=cxebe3d8;
pattern5 v=s c=cxf8c6a1;
pattern6 v=s c=cxf3a267;
pattern7 v=s c=cxf07f32;

legend1 mode=share across=1 position=(bottom right)
 shape=bar(.15in,.15in) offset=(-4,12)
 label=none
 value=(
 t=1 'Best'
 t=2 ''
 t=3 ''
 t=4 ''
 t=5 ''
 t=6 ''
 t=7 'Worst'
 );
 

title1 ls=1.5 "Ranking the US Counties";

/* use a 'note' instead, so it can share what is technically the map's space */
*title2 ls=0.8 h=16pt "Based on: Poverty, Education, and Unemployment";

footnote 
 link='https://www.ers.usda.gov/data-products/county-level-data-sets/download-data/'
 c=gray "Data source: USDA Economic Research Service (2019 snapshot)";

proc gmap data=my_data map=my_map all anno=anno_outline;
id statecode county;
note move=(29,89.5) h=16pt "Based on: Poverty, Education, and Unemployment";
choro overall_rank / levels=7 legend=legend1 
 cdefault=gray99 coutline=graycc
 html=my_html
 des='' name="&name";
run;

title;
footnote;

title c=gray33 h=18pt "Best 20 Overall";
proc print data=my_data (obs=20) label noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt};
label overall_rank='Overall Rank' county_name='County' statecode='State'
 poverty='Poverty' poverty_rank='Poverty Rank' 
 education='College Degree' education_rank='Education Rank' 
 unemployment='Unemployment' unemployment_rank='Unemployment Rank';
var overall_rank county_name statecode poverty poverty_rank education education_rank unemployment unemployment_rank;
run;

title c=gray33 h=18pt "Best 20 Poverty";
proc print data=poverty (obs=20) label noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt}; 
label county_name='County' statecode='State'
 poverty='Poverty' poverty_rank='Poverty Rank';
var poverty_rank county_name statecode poverty;
run;

title c=gray33 h=18pt "Best 20 Education";
proc print data=education (obs=20) label noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt}; 
label county_name='County' statecode='State'
 education='College Degree' education_rank='Education Rank';
var education_rank county_name statecode education;
run;

title c=gray33 h=18pt "Best 20 Unemployment";
proc print data=unemployment (obs=20) label noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt}; 
label county_name='County' statecode='State'
 unemployment='Unemployment' unemployment_rank='Unemployment Rank';
var unemployment_rank county_name statecode unemployment;
run;

quit;
ODS HTML CLOSE;
ODS LISTING;