%let name=coronavirus_2020;

/* 
Set your current-working-directory (to read/write files), if you need to ...
%let rc=%sysfunc(dlgcdir('c:\someplace\public_html')); 
*/
filename odsout '.';

/*
Imitatation/variation of a map from a NYT article:
https://www.nytimes.com/2020/01/21/health/cdc-coronavirus.html
*/

/* 
These Province names must exactly match with what's in the 
id1name variable in mapsgfk.china_attr.
You can see the names using code like this:
 proc sql;
 select unique id1name from mapsgfk.china_attr;
 quit; run;
*/
data china_data;
input cases deaths name $ 8-80;
datalines;
375  9 Hubei Sheng
2    0 Liaoning Sheng
10   0 Beijing Shi
9    0 Shanghai Shi
2    0 Tianjin Shi
5    0 Zhejiang Sheng
26   0 Guangdong Sheng
1    0 Shandong Sheng
1    0 Henan Sheng
2    0 Jiangxi Sheng
1    0 Hunan Sheng
2    0 Sichuang Sheng
5    0 Chongqing Shi
4    0 Hainan Sheng
1    0 Yunnan Sheng
;
run;
/* merge in the id */
proc sql noprint;
create table china_data as
select unique china_data.*, china_attr.id1 as my_id
from china_data left join mapsgfk.china_attr
on china_data.name=china_attr.id1name;
select quote(trim(left(my_id))) into :provinces separated by ' ' from china_data;
quit; run;

/* Had to look up the id1 numbers for the Province names in mapsgfk.china_attr */
data china_map; set mapsgfk.china (where=(id1 in (&provinces)) drop=resolution);
run;
proc gremove data=china_map out=china_map (rename=(id1=my_id));
by id1; id id;
run;

/* Let's add 'China' to the province names, so it will show in the mouse-over text */
data china_data; 
length name $100;
set china_data;
name=trim(left(name))||', China';
run;


/* Now for the country data for non-China countries ... */
data countries_data;
input cases deaths name $ 8-80;
my_id=name;
datalines;
1    0 Japan
1    0 Macau
3    0 Thailand
1    0 South Korea
1    0 Taiwan, Province of China
;
run;
/* get a list of all the countries */
proc sql noprint;
select quote(trim(left(name))) into :countries separated by ' ' from countries_data;
quit; run;

data countries_map (rename=(idname=my_id) drop=resolution);  
 set mapsgfk.world (where=(idname in (&countries)));
run;

/* Combine the maps */
data my_map; set china_map (where=(density<=2)) countries_map (where=(density<=3));
run;

/* 
Combine the attr datasets, and add some labels to the variables 
(these labels show up in the mouse-over text and legend) 
*/
data my_attr; 
label my_id='ID';
label name='Name';
label cases='Coronavirus cases';
label deaths='Coronavirus deaths';
set china_data countries_data;
run;

/* 
Note: there's not a country/province polygon to conveniently use for Macau,
and also it's so small it wouldn't be visible on a worls map, so I left it out.
*/




/*
The country names on the tile-map background generally get obscured
by the overlaid choro polygons, therefore I'm going to overlay some
text for the country names.
I hard-code the country locations, rather than using the %centroid macro,
so you can get them in the exact location you want.
*/
data country_names;
input lat long country_name $ 24-80;
datalines;
33.3684933 137.9375099 Japan
37.8507097 129.7470463 S. Korea
23.3540002 121.6184836 Taiwan
15.3247781  99.1174074 Thailand
35.2143691 105.3870128 China
;
run;



ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Coronavirus 2020") 
 style=htmlblue;

ods graphics /
 noscale /* if you don't use this option, the text will be resized */
 imagemap tipmax=2500 
 imagefmt=png imagename="&name"
 width=750px height=800px noborder; 

title1 "Preliminary map: Japan, Thailand, South Korea, and Taiwan";
proc sgmap mapdata=countries_map maprespdata=countries_data;
esrimap url="http://services.arcgisonline.com/arcgis/rest/services/Canvas/World_Light_Gray_Base";
choromap cases / mapid=my_id;
run;


title1 "Preliminary map: all areas";
proc sgmap mapdata=my_map maprespdata=my_attr;
esrimap url="http://services.arcgisonline.com/arcgis/rest/services/Canvas/World_Light_Gray_Base";
choromap cases / mapid=my_id colormodel=(cornsilk red);
run;

title1 c=gray33 h=18pt "Coronavirus Cases";
title2 c=gray33 h=10pt "(proof-of-concept map, using NYT data snapshot from Jan 22, 2020)";

proc sgmap mapdata=my_map maprespdata=my_attr plotdata=country_names;
esrimap url="http://services.arcgisonline.com/arcgis/rest/services/Canvas/World_Light_Gray_Base";
choromap cases / mapid=my_id colormodel=(cornsilk red)
 tip=(name cases deaths);
text x=long y=lat text=country_name / position=right tip=none
 textattrs=(color=dodgerblue size=12pt weight=bold);
run;

proc sort data=my_attr out=my_attr;
by descending cases;
run;
title;
proc print data=my_attr label noobs; 
var name cases;
sum cases;
run;

quit;
ODS HTML CLOSE;
ODS LISTING;