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

/*
Might create a map like this:
https://twitter.com/TimBroderick/status/971589994820374529

Scraping data from:
https://www.cdc.gov/drugoverdose/maps/rxrate-maps.html
*/

%let st=NC;
%let stname=North Carolina;

%macro getdata(year);

filename dataurl url "https://www.cdc.gov/drugoverdose/maps/rxcounty&year..html";

data opioid_data_&year (drop = whole_line county_line state_line fips_line rate_line);
infile dataurl pad;
year=.; year=&year;
length county_name $50;
length statecode $2;
length fips $5;
input whole_line $ 1-1000;
if index(whole_line,'<tr>')^=0 then do;
 input county_line $ 1-1000;
 county_name=scan(county_line,2,'<>');
 input state_line $ 1-1000;
 statecode=scan(state_line,2,'<>');
 input fips_line $ 1-1000;
 fips=trim(left(scan(fips_line,2,'<>')));
 state=.; state=substr(fips,1,2);
 county=.; county=substr(fips,3,3);
 input rate_line $ 1-1000;
 rate=.; 
 rate=scan(rate_line,2,'<>');
 if index(county_line,'County')=0 then output;
 end;
run;

%mend;

%getdata(2006);
%getdata(2007);
%getdata(2008);
%getdata(2009);
%getdata(2010);
%getdata(2011);
%getdata(2012);
%getdata(2013);
%getdata(2014);
%getdata(2015);
%getdata(2016);

data opioid_data; set 
 opioid_data_2006 
 opioid_data_2007 
 opioid_data_2008 
 opioid_data_2009 
 opioid_data_2010 
 opioid_data_2011 
 opioid_data_2012 
 opioid_data_2013 
 opioid_data_2014 
 opioid_data_2015 
 opioid_data_2016
 ;
run;

data all_counties; set mapsgfk.us_counties_attr;
do year=2006 to 2016;
 output;
 end;
run;

proc sql noprint;
create table opioid_data as
select unique 
 all_counties.year, 
 all_counties.state, 
 all_counties.county, 
 all_counties.statecode,
 all_counties.id1name as state_name,
 all_counties.idname as county_name,
 opioid_data.rate
from all_counties left join opioid_data
on 
all_counties.year=opioid_data.year and
all_counties.state=opioid_data.state and
all_counties.county=opioid_data.county;
quit; run;

proc sort data=opioid_data out=opioid_data;
by year state county;
run;

data opioid_data; set opioid_data;
if rate>125 then rate_bucket=6;
else if rate>100 then rate_bucket=5;
else if rate>75 then rate_bucket=4;
else if rate>50 then rate_bucket=3;
else if rate>1 then rate_bucket=2;
else rate_bucket=1;
run;

proc format;
 value ratefmt
 1 = 'no data'
 2 = '<=50'
 3 = '50-75'
 4 = '75-100'
 5 = '100-125'
 6 = '>125'
 ;
run;

data opioid_data; set opioid_data (where=(statecode="&st"));
length my_html $300;
my_html='title='||quote(trim(left(county_name)));
run;

data my_map; set mapsgfk.us_counties (
 where=(
 statecode="&st"
 and 
 density<=3
 ) 
 drop=resolution);
run;

proc gproject data=my_map out=my_map latlong eastlong degrees dupok
 parmout=work.projparm;
id state county;
run;

data anno_cities; set mapsgfk.uscity (where=(statecode="&st" and pop_type='over100k'));
run;

proc gproject data=anno_cities out=anno_cities latlong eastlong degrees dupok
 parmin=work.projparm parmentry=my_map;
id;
run;

data anno_cities; set anno_cities;
length function $8 color $12 style $35;
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; rotate=360; size=0.8; 
style='psolid'; color='gray88'; output;
style='pempty'; color='gray33'; output;
function='label'; position='2'; size=2.2; style=''; text=trim(left(city)); 
if city='Greensboro' then position='8';
output;
run;


options dev=sasprtc printerpath=gif animduration=1.0 animloop=0 
 animoverlay=no animate=start center nobyline;

goptions xpixels=900 ypixels=480;
goptions border;

ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Opioid Prescription Rate") 
 style=htmlblue;

goptions gunit=pct htitle=6.5 htext=3.2 ftitle="albany amt/bold" ftext="albany amt";
goptions ctext=gray33;

title1 ls=2.5 "Opioid Prescription Rate in &stname Counties";
title2 ls=0.8
 link="https://www.cdc.gov/drugoverdose/maps/rxrate-maps.html"
 c=gray "Data source: Centers for Disease Control";
footnote;
proc gslide / des='' name="&name";
run;

pattern6 v=s c=cxd7191c;
pattern5 v=s c=cxfdae61;
pattern4 v=s c=cxffffbf;
pattern3 v=s c=cxa6d96a;
pattern2 v=s c=cx1a9641;
pattern1 v=s c=graydd;

legend1 label=(position=top font='albany amt/bold' 
 j=c "Prescriptions per" 
 j=c "100 persons:") 
 value=(justify=center) position=(bottom left inside) mode=share
 across=1 order=descending shape=bar(.15in,.15in) offset=(5,1);

options nobyline;

proc gmap data=opioid_data map=my_map all anno=anno_cities;
by year;
note move=(35,10) font='albany amt/bold' height=13 "#byval(year)";
id state county;
format rate_bucket ratefmt.;
choro rate_bucket / midpoints=(1 2 3 4 5 6)
 coutline=gray88 cdefault=cyan 
 legend=legend1
 des='' name="&name";
run;

/* repeat the last frame a couple of extra times, for a 'delay' before repeating */
proc gmap data=opioid_data (where=(year=2016)) map=my_map all anno=anno_cities;
by year;
note move=(35,10) font='albany amt/bold' height=13 "#byval(year)";
id state county;
format rate_bucket ratefmt.;
choro rate_bucket / midpoints=(1 2 3 4 5 6)
 coutline=gray88 cdefault=cyan
 legend=legend1
 des='' name="&name";
run;
run;

/* Now, do 1 frame with the html mouse-over text for county names */
proc gmap data=opioid_data (where=(year=2016)) map=my_map all anno=anno_cities;
by year;
note move=(35,10) font='albany amt/bold' height=13 "#byval(year)";
id state county;
format rate_bucket ratefmt.;
choro rate_bucket / midpoints=(1 2 3 4 5 6)
 coutline=gray88 cdefault=cyan
 legend=legend1
 html=my_html
 des='' name="&name";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;