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

%let minyear=1900;

%let landcolor=cxf1e9e0;
%let outcolor=gray99;

/* For color-binning of magnitude */
%let upper=7.5;
%let lower=4.0;

/* transparent version */
%let color1=A78c92199;   /* weakest/smallest magnitude  */
%let color2=A007fff44;   /*                             */
%let color3=Aff000066;   /* strongest/biggest magnitude */

/* solid version */
%let color1s=cx78c921;   /* weakest/smallest magnitude  (green) */
%let color2s=cx007fff;   /*                             (blue) */
%let color3s=cxff0000;   /* strongest/biggest magnitude (red) */

/* Reduce density, so it's easier to work with */
data whole_map; set maps.world (where=(density<=1 and cont^=97) drop=x y);
 x=long;
 y=lat;
run;

/* Use gproject to clip-and-close polygons for left & right side */
proc gproject data=whole_map out=right_map dupok project=none
  longmin=0 longmax=4;
  id id;
run;
proc gproject data=whole_map out=left_map dupok project=none
  longmin=-4 longmax=0;
  id id;
run;

/*
Move the right half of the map to the left, and make the segment
variables unique, since some segments got split into two.
*/
data right_map; set right_map;
segment=segment+5000;
x=x-(3.14159*2);
run;
data left_map; set left_map;
run;

/* Put map back together. */
data my_map; set right_map left_map;
run;


proc sql;
create table my_mapdata as
select unique id, '1' as colorval
from my_map;
quit; run;


/* 
Go to this page:
http://www.ngdc.noaa.gov/nndc/struts/form?t=101650&s=1&d=1
Scroll to bottom and select
"Return All Selected Events in tab-delimited format for import into Excel"
Click 'Search database'
Then Save-as
worldquakes.txt
*/
PROC IMPORT OUT=quake_data DATAFILE="worldquakes.txt" DBMS=TAB REPLACE;
 GETNAMES=YES;
 DATAROW=2; 
RUN;

data quake_data; set quake_data (where=(year>=&minyear));
datestring=put(month,z2.)||'/'||put(day,z2.)||'/'||trim(left(year));
format date date9.;
date=input(datestring,mmddyy10.);
datestring2=lowcase(put(date,date9.));
/* 
Over the years, magnitudes have been measured in different ways.
Here's the priority list of the order in which I use them.
*/
magnitude=.;
if eq_mag_ms^=. then magnitude=eq_mag_ms;
else
if eq_mag_mw^=. then magnitude=eq_mag_mw;
else
if eq_mag_mb^=. then magnitude=eq_mag_mb;
else
if eq_mag_ml^=. then magnitude=eq_mag_ml;
else
if eq_mag_mfa^=. then magnitude=eq_mag_mfa;

/* Convert eastlong degrees to westlong radians */
x=atan(1)/45 * -1*longitude;
y=atan(1)/45 * latitude;

/* 
move data that was to left of prime meridian (n & s america, etc)
to the right of the map 
*/
if x>0 then x=x-(3.14159*2);
if x^=. and y^=. then output;
run;

/*
proc sort data=quake_data out=quake_data;
 by descending magnitude;
run;
*/

proc sql;
select min(year) into :minyear from quake_data;
select max(year) into :maxyear from quake_data;
select unique datestring2 into :maxdate from quake_data having date=max(date);
quit; run;
%let minyear=%trim(&minyear);
%let maxyear=%trim(&maxyear);


/* Create an annotate dataset, containing pie/circle for each location */
data circle_anno; set quake_data; 
anno_flag=1;
run;

/* combine, project, and separate */
data combined; set my_map circle_anno; run;
proc gproject data=combined out=combined dupok project=cylindri;
  id id;
run;
data my_map circle_anno;
  set combined;
  if anno_flag=1 then output circle_anno;
  else output my_map;
run;


%let max_area=15;
proc sql noprint;
select max(magnitude) into :max_val from circle_anno;
quit; run;

data circle_anno; set circle_anno; 
length function $8 style color $12 text $20 html $1024;
xsys='2'; ysys='2'; hsys='3'; when='a';
if magnitude lt &lower then color="&color1";
else if magnitude lt &upper then color="&color2";
else color="&color3";
function='pie'; style='pempty'; rotate=360; 
size=sqrt((magnitude/&max_val)*&max_area/3.14);


run;


/* Now, add a "fake legend" using annotated stuff ... */
data legend; 
 length function $8 text $20 style $35;
 xsys='3'; ysys='3'; hsys='3'; when='A';
 function='PIE'; line=0; angle=0.0; rotate=360.0; size=1.5; style='SOLID';
 /* add circles ... */
 size=sqrt(((&upper+1)/&max_val)*&max_area/3.14); 
 color="&color3s"; x=3; y=89; style='pempty'; output; 
 size=sqrt((((&upper+&lower)/2)/&max_val)*&max_area/3.14); 
 color="&color2s";  x=3; y=84-.5; style='pempty'; output; 
 size=sqrt(((&lower-1.5)/&max_val)*&max_area/3.14); 
 color="&color1s";  x=3; y=79-.5; style='pempty'; output; 
 /* add text ... */
 function='LABEL'; position='6'; color='dag';
 style='albany amt/bold'; size=2.8;
 text='Magnitude:';  x=2;  y=95; output;
 style='albany amt';  size=2.5;
 text=">= &upper";  x=5; y=89; output;
 text=" < &upper";  x=5; y=84; output;
 text=" < &lower";  x=5; y=79; output;
run;

/* Move the annotated legend to a better location */
data legend; set legend;
 y=y-75;
 x=x+18;
run;


goptions device=png;
goptions xpixels=1200 ypixels=600;
goptions border cback=white;
 
ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Significant Earthquakes") options(pagebreak='no') 
 style=htmlblue;

goptions gunit=pct htitle=4 htext=2.5 ftitle="albany amt" ftext="albany amt" ctitle='black' ctext='black';

pattern1 color=&landcolor v=s;

footnote c=gray j=right link="http://www.ngdc.noaa.gov/nndc/struts/form?t=101650&s=1&d=1"
 "Data from www.ngdc.noaa.gov (as of &maxdate)";

title1 ls=1.5 "Significant Earthquakes";
title2 "Date Range:  year &minyear to &maxdate";
proc gmap map=my_map data=my_mapdata anno=circle_anno;
id id; 
choro colorval / levels=1 coutline=&outcolor  nolegend discrete 
 anno=legend des='' name="&name";
run;

proc sql noprint;
create table summary as
select unique year, count(*) as count
from quake_data
group by year;
quit; run;
data summary; set summary;
length my_html $300;
my_html='title='||quote(trim(left(year))||': '||trim(left(count)));
run;
title1 ls=1.5 "Count of Significant Earthquakes per year";
title2 "Date Range:  year &minyear to &maxdate";
symbol1 v=circle h=1.5 i=needle c=blue;
axis1 minor=none offset=(0,0);
axis2 minor=none order=(&minyear to 2025 by 25) offset=(0,0);
proc gplot data=summary;
plot count*year /
 vaxis=axis1 haxis=axis2 noframe vzero
 autovref cvref=graydd
 html=my_html
 des='' name="&name._graph";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;