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

/* Run this once, and after that use the copy you saved the 1st time.
   This reduces the resolution (and therefore the number of obsns) 
   in the map, and makes the program run a lot faster. */
libname here '.';
/*
proc greduce data=maps.asia out=fooasia;
id id; run;
proc greduce data=maps.spacific out=foospacific;
id id; run;
data here.typhmap; 
 set fooasia (where = (density <= 2))  
     foospacific (where = (density <= 2)); 
run;
*/

proc sql;
create table typhmap as
 select -1*long as x, lat as y, segment, id
from here.typhmap;


/* Add html title= chart tips with the country name */
proc sql;
create table my_mapdata as
select unique id
from typhmap;
quit; run;
options fmtsearch=(sashelp.mapfmts);
data my_mapdata; set my_mapdata;
length cname $ 50;
cname=put(id,glcnsm.);
if cname in ('Indonesia' 'Sri Lanka' 'India' 'Malaysia' 'Thailand') then red_one='y'; else red_one='n';
length htmlvar $200;
htmlvar=
'title='||quote( trim(left(cname))|| ' ')|| ' '||
'href='||quote( 'http://www.google.com/search?hl=en&q=earthquake+tsunami+'||trim(left(cname))||' ');
run;


data dotloc;
 lat=3.3000; /* deg north */
 long=95.7800; /* deg East */
run;

data dotloc; set dotloc;
  x=atan(1)/45 * long; 
  y=atan(1)/45 * lat;
run;

data circle_anno; 
   length function style color $ 8 position $ 1 text $ 20 html $1024;
   retain xsys ysys '2' hsys '3' when 'a';
  set dotloc; 
  anno_flag=1;

/* This is the html title= charttip and href= drilldown for the annotated circle */
html=
'title='||quote(
  'Earthquake: 8.9 magnitude, 7am 26dec2004 '||'0d'x||
  'Off coast of Indonesia Sumatra Island '||'0d'x||
  'Epicenter '||'0d'x||
  'Latitude: '||trim(left(lat))||'0d'x||
  'Longitude: '||trim(left(long))||
  ' ')||
' '||
'href="http://www.cnn.com/SPECIALS/2004/tsunami.disaster/"';

  position='5'; /* centered at this z/y coordinate */
  function='pie'; 
  line=0; /* only draw outside of circle arc - don't draw the radius */

 /* The width= is new for v9.2 -- in older versions it is ignored */
  color='white'; width=2; size=2; rotate=360; output;

  /* after the 1st inner-most circle, take away the html charttip/drilldown */
  html='';

  color='white'; width=2; size=size+4; rotate=360; output;

  color='white'; width=2; size=size+6; rotate=360; output;

  color='white'; width=2; size=size+9; angle=75; rotate=252; output;
  color='graybb'; width=1; angle=75+252; rotate=360-252; output;

  color='white'; width=2; size=size+11; angle=95; rotate=220; output;
  color='grayaa'; width=1; angle=95+220; rotate=360-220; output;

run;


/* Annotate the longitude/latitude grid lines */
%let min_long=70;
%let max_long=125;
%let min_lat=-15;
%let max_lat=20;
data grid_anno;
 length function style color $ 8;
 retain xsys ysys '2' hsys '3' when 'a';
 anno_flag=1;
 do long = &min_long to &max_long by 5;
   function='label'; color='white'; 
   longitude=long; latitude=&min_lat;
   degrees=abs(int(long));
   minutes=abs((long-int(long))*60);
   text=trim(left(degrees)); 
   text=trim(left(degrees))||'B0'x;  
   position='a'; size=2.5; angle=90; output;
   color='white'; size=.1; line=33; 
   function='move'; longitude=long; latitude=&min_lat; output;
   do lat = &min_lat to &max_lat by 5;
     function='draw'; longitude=long; latitude=lat; output;
   end;
 end;
 do lat = &min_lat to &max_lat by 5;
   function='label'; color='white'; 
   longitude=&max_long; latitude=lat;
   degrees=abs(int(lat));
   minutes=abs((lat-int(lat))*60);
   text=trim(left(degrees)); 
   text=trim(left(degrees))||'B0'x; 
   position='c'; size=2.5; angle=0; output;
   color='white'; size=.1; line=33; 
   function='move'; latitude=lat; longitude=&min_long; output;
   do long = &min_long to &max_long by 5;
     function='draw'; latitude=lat; longitude=long; output;
   end;
 end;
run;
data grid_anno; set grid_anno;
 x=atan(1)/45 * longitude;
 y=atan(1)/45 * latitude;
run;

data all_anno;
 set circle_anno grid_anno;
run;

/* project the map */
/* See examples of projections here ... 
   http://www.progonos.com/furuti/MapProj/Normal/ProjTbl/projTbl.html */
data combined; set typhmap all_anno; run;
proc gproject data=combined out=combined dupok
  eastlong
  project=robinson
  longmax=125
  longmin=65
  latmax=25
  latmin=-15
  ;
  id id;
run;
data my_map all_anno;
  set combined;
  if anno_flag=1 then output all_anno;
  else output my_map;
run;


GOPTIONS DEVICE=png;
ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm"
 (title="SAS/Graph gmap - Indonesia Earthquake (26dec2004)") style=minimal;

goptions 
 gunit=pct 
 htitle=5 htext=3.5
 ftitle="arial" ftext="arial"  
 ctitle='white' ctext='white'
;
goptions border cback=cx35b2e0;

* goptions xpixels=800 ypixels=700;
goptions xpixels=725 ypixels=600;

title "Indonesian Earthquake (26dec2004)";
title2 link='http://www.google.com/search?hl=en&q=indonesia+earthquake+tsunami+2004'
  "Click here to see web pages...";

/* Getting a little tricky with this title ;) */
/*
title3 h=4pct a=-90 r=180 "edutitaL  ";
*/
footnote h=8pct " ";

pattern1 v=s c=gray;  
pattern2 v=s c=red;
/*
pattern1 v=s c=cx10CA10;  
pattern2 v=s c=cx417231;
*/

proc gmap map=my_map data=my_mapdata anno=all_anno;
id id; 
choro red_one / 
 coutline=black  
 nolegend 
 discrete 
 html=htmlvar
 des="" name="&name";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;