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

/* Written by Robert Allison (Robert.Allison@sas.com) */

/*
Reproduction of Minard's famous map of Napoleon's march on Moscow,
similar to the following...

Links:
http://www.math.yorku.ca/SCS/Gallery/re-minard.html  (links to many versions)

http://www.spss.com/research/wilkinson/TheGrammarOfGraphics/minard.txt    <-- data
http://www.spss.com/research/wilkinson/TheGrammarOfGraphics/GOG.html

http://www.math.yorku.ca/SCS/Gallery/minard/
http://www.math.yorku.ca/SCS/Gallery/minard/Minard-IML.jpg (lat/long!)
http://www.math.yorku.ca/SCS/Gallery/minard/NapoleonsMarch.iml (iml code w/ lat/long)
http://www.math.yorku.ca/SCS/Gallery/minard/march-animated.gif
http://www.math.yorku.ca/SCS/Gallery/minard/minard-nvizn.gif
http://www.math.yorku.ca/SCS/Gallery/minard/minard-odt.jpg

http://www.edwardtufte.com/tufte/posters
http://www.edwardtufte.com/tufte/minard
http://www.csiss.org/classics/content/58
*/

goptions reset=global;

options fmtsearch=(sashelp.mapfmts);
data my_map;
 set maps.europe maps.asia;
 /* we want the unprojected values, since we're projecting them ourselves */
 /* Also convert longitude to a proper 'eastlong' value, so it matches 
    the annotate long/lat values */
 x=long*-1;  
 y=lat;
 country=id;  /* This just makes the code easier to follow */
   length  countryname $20;
   countryname=put(country,glcnsm.);
run;
proc sql;
create table countries as 
select unique country, countryname
from my_map;
quit; run;
data countries; set countries;
length myhtml $ 1024;
myhtml= 'title='|| quote( trim(left(countryname))||' ');
run;


/*
Got data from...
http://www.spss.com/research/wilkinson/TheGrammarOfGraphics/minard.txt

I added one extra obsn for group=2 during the river-crossing, so it
was more evident that the 2 groups combined...
28.3  54.4   28000  R  2
*/
data armyline;
length direc $1;
length direction $10;
input lon   lat surviv direc group ;
  n+1;
  if direc eq 'A' then direction='advance';
  else if direc eq 'R' then direction='retreat';
  /* Convert degrees to radians */
  x=atan(1)/45 * lon;
  y=atan(1)/45 * lat;
datalines;
24.0  54.9  340000  A  1
24.5  55.0  340000  A  1
25.5  54.5  340000  A  1
26.0  54.7  320000  A  1
27.0  54.8  300000  A  1
28.0  54.9  280000  A  1
28.5  55.0  240000  A  1
29.0  55.1  210000  A  1
30.0  55.2  180000  A  1
30.3  55.3  175000  A  1
32.0  54.8  145000  A  1
33.2  54.9  140000  A  1
34.4  55.5  127100  A  1
35.5  55.4  100000  A  1
36.0  55.5  100000  A  1
37.6  55.8  100000  A  1
37.5  55.7   98000  R  1
37.0  55.0   97000  R  1
36.8  55.0   96000  R  1
35.4  55.3   87000  R  1
34.3  55.2   55000  R  1
33.3  54.8   37000  R  1
32.0  54.6   24000  R  1
30.4  54.4   20000  R  1
29.2  54.4   20000  R  1
28.5  54.3   20000  R  1
28.3  54.4   20000  R  1
24.0  55.1   60000  A  2
24.5  55.2   60000  A  2
25.5  54.7   60000  A  2
26.6  55.7   40000  A  2
27.4  55.6   33000  A  2
28.7  55.5   30000  A  2
29.2  54.3   30000  R  2
28.5  54.2   30000  R  2
28.3  54.3   28000  R  2
28.3  54.4   28000  R  2
27.5  54.5   20000  R  2
26.8  54.3   12000  R  2
26.4  54.4   14000  R  2
24.6  54.5    8000  R  2
24.4  54.4    4000  R  2
24.2  54.4    4000  R  2
24.1  54.3    4000  R  2
24.0  55.2   22000  A  3
24.5  55.3   22000  A  3
24.6  55.8    6000  A  3
24.2  54.4    6000  R  3
24.1  54.3    6000  R  3
;
run;

proc sort data=armyline out=armyline;
 by group n;
run;

%let maxwidth=3;   /* Maximum width of line */
%let maxdot=1.5;   /* Maximum size (radius) of vertex fill-in dot */
data armyline;
   length function style color $ 12 position $ 1 text $ 20 html $1024;
   retain xsys ysys '2' hsys '3' when 'a';
  set armyline;
  by group;
  anno_flag=1;
  if first.group then do;
    function='move';
    output;
    end;
  else do;
    if direc eq 'A' then color='cxaddd8e';
    else if direc eq 'R' then color='cxfe0000';
    function='draw';
    size=(surviv/340000) * &maxwidth;
    output;
    /* Vertex fill-in dot -- this gives smoother transition between line segments, like a kneecap/elbow, 
       and also gives a place for html charttips (since line segments don't support charttips) */
    size=(surviv/340000) * &maxdot;
    function='pie'; style='psolid'; rotate=360; 
    length html $ 100;
    html= 'title='|| quote( trim(left(put(surviv,comma7.0))||' survivors during '||trim(left(direction)) )||' ');
    output;
    end;
run;

data drillbox;
   length function style color $ 12 position $ 1 text $ 20 html $1024;
   retain xsys ysys '2' hsys '3' when 'a';
  anno_flag=6;
  size=.1;
  line=2;
  color='black';
  html= 'title='||quote('zoom in')||' '||'href="zoom.htm"';
  function='poly';
  x=atan(1)/45 * 23.4;
  y=atan(1)/45 * 52.0;
  output;
  function='polycont';
  x=atan(1)/45 * 39.0;
  y=atan(1)/45 * 52.0;
  output;
  x=atan(1)/45 * 39.0;
  y=atan(1)/45 * 56.5;
  output;
  x=atan(1)/45 * 23.4;
  y=atan(1)/45 * 56.5;
  output;
run;


data cities;
input lon lat city $ 11-41;
  /* Convert degrees to radians */
  x=atan(1)/45 * lon;
  y=atan(1)/45 * lat;
/*
30.4 53.9 Mohilow               :
27.6 53.9 Minsk                 :
*/
datalines;
24.0 55.0 Kowno          
25.3 54.7 Wilna         
26.4 54.4 Smorgoni     
26.8 54.3 Molodexno   
27.7 55.2 Gloubokoe      
28.5 54.3 Studienska   
28.7 55.5 Polotzk     
29.2 54.4 Bobr       
30.2 55.3 Witebsk                  
30.4 54.5 Orscha                  
32.0 54.8 Smolensk              
33.2 54.9 Dorogobouge          
34.3 55.2 Wixma               
34.4 55.5 Chjat              
36.0 55.5 Mojaisk           
37.6 55.8 Moscou           
36.6 55.3 Tarantino       
36.5 55.0 Malo-jarosewli 
;
run;
/* Annotate a black dot/pie at each city, with mouseover html charttip
   showing the city name.  (This variable *must* be called 'html'.) */
data cities;
   length function style color $ 12 position $ 1 text $ 20 html $1024;
   retain xsys ysys '2' hsys '3' when 'a';
  set cities;
length html $ 100;
html= 'title='|| quote( trim(left(city))||' ');
   function='pie';
   color='black';
   style='psolid';
   position='5';
   rotate=360;
   size=.3;
   anno_flag=2;
   output;
run;

/* Blue rectangle behind map area (ie, water) */
/* Along certain edges, I had to use these .01 offsets,
   because gproject clipping was trimming off the edge otherwise */
data blue_anno;
 length function style color $ 8;
 retain xsys ysys '2' hsys '3' when 'b';
 anno_flag=3;
 color='cx35b2e0';
 style='msolid';
 function='poly';
  long=5; lat=39+.01; output;
 function='polycont';
 lat=39+.01;
 do long=5 to 40 by 5;
   output;
 end;
 long=39.99;
 do lat=39+5 to 60 by 5;
   output;
 end;
 lat=60;
 do long=40 to 5 by -5;
   output;
 end;
 long=5;
 do lat=60 to 5+5 by -5;
   output;
 end;
 long=5; lat=39+.01; output;
run;
data blue_anno; set blue_anno;
 x=atan(1)/45 * long;
 y=atan(1)/45 * lat;
run;


data names;
input lon lat countryname $ 11-31;
  /* Convert degrees to radians */
  x=atan(1)/45 * lon;
  y=atan(1)/45 * lat;
datalines;
32.0 57.5 Russia                
 6.4 48.8 France                
10.0 51.5 Germany               
19.0 52.5 Poland                
27.5 53.2 Belarus               
24.0 56.0 Lithuania             
;
run;
data names;
   length function style color $ 12 position $ 1 text $ 20 html $1024;
   retain xsys ysys '2' hsys '3' when 'a';
  set names;
   function='label';
   text=trim(left(countryname));
   color='black';
   style='"arial"';
   position='5';
   size=3;
   anno_flag=4;
run;





data combined;
 set 
 my_map
 armyline
 cities
 blue_anno
 names
 drillbox
 ;
run;

proc gproject data=combined out=combined dupok
  eastlong
  project=robinson
  latmax=60
  latmin=39
  longmin=5
  longmax=40
  ;
  id country;
run;

data my_map armyline cities blue_anno names drillbox;
  set combined;
  if anno_flag=1 then output armyline;
  else if anno_flag=2 then output cities;
  else if anno_flag=3 then output blue_anno;
  else if anno_flag=4 then output names;
  else if anno_flag=6 then output drillbox;
  else output my_map;
run;

data img;
 length function text $ 8; length color $ 8;
 xsys='3'; ysys='3'; hsys='3';
 when='A';
 html='title="Tell me more about SAS/Graph..." href="http://www.sas.com/technologies/bi/query_reporting/graph/factsheet.pdf" ';
 function='move'; x=82.5; y=5.5; output;
 function='image'; x=x+11; y=y+6; imgpath='powered_stat.gif'; style='fit'; output;
run;


data anno1; 
 set blue_anno drillbox armyline cities names img;
run;


 GOPTIONS DEVICE=gif;
 
 ODS LISTING CLOSE;
 ODS HTML path=odsout body="&name..htm"
 (title="Napoleon's March 50,000-mile view (SAS Version of Minard's Map)")
 style=minimal
 gtitle nogfootnote
 ;

 goptions noborder;
 goptions gunit=pct htitle=5 htext=3 ftitle="arial/bo" ftext="arial";

 pattern1 v=s c=tan r=500;
 title "Napoleon's Russian Campaign, 1812";
 title2 "Plotted on modern map";
 footnote h=10pt f="arial" "Mouse over dots and campaign path to see more info";
 proc gmap map=my_map data=countries anno=anno1;
  id country;
  choro countryname / 
  levels=1
  coutline=gray
  nolegend
  html=myhtml
  des="" name="&name";
 run;

 quit;
 ODS HTML CLOSE;
 ODS LISTING;