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

/*
This example uses data/scenario from a contest on kaggle.com
https://www.kaggle.com/c/predict-west-nile-virus/data?west_nile.zip
*/

PROC IMPORT OUT=train (rename=(longitude=lon_deg latitude=lat_deg))
 DATAFILE="train.csv"
 DBMS=CSV REPLACE;
 GETNAMES=YES;
 DATAROW=2;
RUN;

%let datef=01aug2007;

data train; set train;
date=round(date);
if date="&datef"d then output;
run;


/* get the map tiles ... */

%let zoomlvl=11;
%let min_x=524;
%let max_x=526;
%let min_y=760; /* top of map */
%let max_y=762; /* bottom of map */

/*
%let zoomlvl=12;
%let min_x=1048;
%let max_x=1051;
%let min_y=1520; 
%let max_y=1525; 
*/

/*
http://maps.stamen.com
Attribution:
For Toner: Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL. 
For everything else: Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA. 
%let tileserv=a.tile.stamen.com/toner;
*/
%let tileserv=a.tile.stamen.com/terrain-lines;

%include '../data_vault/proxy.sas';
/*
I %include a file that sets the my_proxy macro variable.
You won't be able to use my proxy, therefore set your own
macro variable, which will be something like the following:
%let my_proxy=http://yourproxy.com:80;
*/

/* 
Programmatically read in those map images, and save them
to my hard drive so I can later annotate them (it would be
nice if annotate could just refer to them without actually
downloading them ... if anyone figures out how to do that, 
let me know! :)
(This bit of code was adapted from some code Jim Metcalf
was using to download a spreadsheet from the web.  Note that
your 'proxy' url will be different, for your own Intranet.))
*/
%macro download_map(zoom, x, y);
filename osm_pic url "http://&tileserv./&zoom./&x./&y..png" 
 proxy="&my_proxy";
data _null_;
 n=-1;
 infile osm_pic recfm=s nbyte=n length=len _infile_=tmp;
 input;
 file "opnsta_&zoom._&x._&y..png" recfm=n;
 put tmp $varying32767. len;
run;   
%mend;


/* loop through all the desired map tile x/y numbers, and call the macro for each */
/* 
If you're re-running the same map several times, you probably want to comment
this out, after the first time, since the png files are already downloaded.
*/
data _null_;
do loop_x=&min_x to &max_x;
 do loop_y=&min_y to &max_y;
   call execute('%download_map('|| &zoomlvl ||','|| loop_x ||','|| loop_y ||');');
   call execute('run;');
  end;
 end;
run;
/*
*/


/*
Create a blank rectangular map area of the same proportions as your images.
Use the OpenStreetMap tile numbers.
*/
data blank_map;
x=&min_x;   y=&min_y;   output;
x=&max_x+1; y=&min_y;   output;
x=&max_x+1; y=&max_y+1; output;
x=&min_x;   y=&max_y+1; output;
run;


/*
Create an annotate data set to 'stitch' the images together.
*/
data anno_maps;
length function $8 color $20 style $35;
xsys='2'; ysys='2'; when='b'; style='fit';
do loop_x=&min_x to &max_x;
 do loop_y=&min_y to &max_y;
  imgpath="opnsta_&zoomlvl._"||trim(left(loop_x))||"_"||trim(left(loop_y))||".png";
  /* the y-coordinates might seem "reversed", but you will be "flipping" them later */
  function='move';  x=loop_x;   y=loop_y+1; output;
  function='image'; x=loop_x+1; y=loop_y;   output;
  end;
 end;
run;

/*
Convert the long & lat degrees to x & y in the OSM tile-based coordinate system,
using equations from:
http://wiki.openstreetmap.org/wiki/Tilenames 
SAS doesn't have a sec() function, so I used 1/cos(value) 
*/
data train; set train;
lon_rad=(atan(1)/45)*lon_deg;
lat_rad=(atan(1)/45)*lat_deg;
n=2**&zoomlvl;
x = ((lon_deg + 180) / 360) * n;
y = (1 - (log(tan(lat_rad) + (1/cos(lat_rad))) / constant('pi') )) / 2 * n;
run;

/*
proc sort data=train out=train;
by WnvPresent;
run;
*/

proc sql;
create table train_summary as
select unique date, trap, AddressNumberAndStreet, lat_deg, lon_deg, 
 sum(NumMosquitos) as NumMosquitos, sum(WnvPresent) as WnvPresent, x, y
from train
group by date, trap
order by NumMosquitos descending;
quit; run;

%let green=A9AFF9Acc;
%let olive=ABCEE68cc;
%let pink=cxFF0000;

/* Annotate red 'dot' at each location */
data anno_mosquitos; set train_summary;
length function $8 color $20 style $35;
xsys='2'; ysys='2'; hsys='3'; when='a'; 
function='pie'; rotate=360; style='psolid'; 
size=.15+sqrt(NumMosquitos/3.14)*.3;
color="&green"; 
if WnvPresent>0 then color="&olive";
output;
/* Rather than html for each bubble, we just do a html for each trap (because of gif animation */
/* hopefully the last frame has an obsn for each trap! */
length html $300;
html='title='||quote(
  'Trap: '||trim(left(trap))||'0d'x||
/*
  'Mosquitos: '||trim(left(NumMosquitos))||'0d'x||
  'WnvPresent: '||trim(left(WnvPresent))||'0d'x||
*/
  trim(left(AddressNumberAndStreet))
  );
style='pempty'; color='gray33'; output;
run;

data anno_wnv; set anno_mosquitos (where=(WnvPresent>0));
size=.15+sqrt(WnvPresent/3.14)*.3;
if style='psolid' then color="&pink";
html='';
run;

data anno_mosquitos; set anno_mosquitos anno_wnv;
run;


/* Add html drilldown */
data anno_mosquitos; set anno_mosquitos;
/*
length html $500;
if style='pempty' then 
html=
 'title='||quote(trim(left(area))||': '||trim(left(description)))||
 ' href='||quote('https://maps.google.com/maps?'||
  'hl=en&ll='||trim(left(lat_deg))||','||trim(left(lon_deg))||
  '&spn=0.000939,0.000916&sspn=0.084982,0.117245&t=h&z=20');
*/
run;


/* 
Flip the y-value, since SAS/GRAPH's coordinate system's origin is at bottom/left,
and OpenStreetMaps' coordinate system's origin is at the top/left
*/
data blank_map; set blank_map;
y=-y;
run;
data anno_maps; set anno_maps;
y=-y;
run;
data anno_mosquitos; set anno_mosquitos;
y=-y;
run;

/*
Some of the dots are outside of the plot axes area I'm showing.
This generates an 'error' in annotate, and once a certain limit
is reached, annotate stops drawing.  I bump that limit way up,
to make sure all my annotated markers and arrows are drawn :)
*/
data ignore_error;
function='seterror'; size=100000; output;
run;


data foo;
zoomlvl=&zoomlvl;
min_x=&min_x;
min_y=&min_y;
max_x=&max_x;
max_y=&max_y;
run;
proc sql noprint;
/* gplot axis length, in inches (each tile is 256x256 pixels, and output is 96ppi) */
select ((max_x-min_x+1)*256)/96 into :x_length from foo;
select ((max_y-min_y+1)*256)/96 into :y_length from foo;
/* size for entire page ... needs to be bigger, to hold titles & axis text, etc */
select (((max_x-min_x+1)*256))+30 into :x_pixels from foo;
select (((max_y-min_y+1)*256))+60 into :y_pixels from foo;
quit; run;



goptions xpixels=&x_pixels ypixels=&y_pixels;
goptions noborder; 

goptions device=png;

ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Mosquitos in Chicago") style=htmlblue;

goptions gunit=pct ftitle='albany amt/bold' ftext='albany amt/bold' htitle=4.0 htext=10pt;

/* the 4 corner points will be basically 'invisible' */
symbol1 value=point interpol=none color=black;

title;
footnote;

/*
Draw the plot & map, and suppress all the axis & refline stuff
that makes it look like a gplot (those things help you see what's going on
while you're first setting the map up, but for the final map you'll probably 
want to remove them (as shown below), so it just looks like a map & not a plot.
(style=0 would also suppress the border/axis, if you don't want that)
*/

data anno_text;
length function $8 color $20 text $50 style $35;
xsys='3'; ysys='3'; hsys='3'; when='a';

function='label'; color='black'; style='albany amt/bold'; position='5'; 
x=75; y=96;  
size=3.0;
y=y-3.5; text="Chicago, IL"; output;
size=3.0;
y=y-4.0; text="Mosquito Monitoring"; output;
size=12.0;
y=y-6; text=substr("&datef",6,9); output;
size=3.5;
y=y-9.5; text=substr("&datef",1,5); output;

function='move'; x=69.0; y=42.0; output;
function='image'; x=x+25; y=y-19.6; imgpath='moggie2.png'; style='fit'; output;

run;

/* Draw a white transparent curtain in front of the map, so it will be 'lighter' */
data anno_curtain;
length color $12 style $35;
xsys='1'; ysys='1'; when='a';
function='move'; x=0; y=0; output;
function='bar'; x=100; y=100; style='solid'; line=0; color="Affffffcc"; output;
run;

/* Now, add a "fake legend" using annotate ... */
data anno_legend;
length function color $12 text $50 style $35;
xsys='3'; ysys='3'; hsys='3'; when='A';
function='pie'; rotate=360.0; 

x=64; y=61; style='psolid';
y=y-3.8; size=1.5; color="&green"; output;
y=y-3.8; size=0.6; color="&pink"; output;

x=64; y=61; style='pempty';
y=y-3.8; size=1.5; color="gray33"; output;
y=y-3.8; size=0.6; color="gray33"; output;

function='label'; position='>'; rotate=.; size=2.5; color='';
x=64+3; y=61; style='albany amt'; 
y=y-3.8; text='Total mosquitos'; output;
y=y-3.8; text='West Nile mosquitos'; output;
run;

/* combine the markers and text ... into anno_stuff */
data anno_stuff; set ignore_error anno_curtain anno_mosquitos (where=(date="&datef"d)) anno_legend anno_text;
run;

axis1 color=gray label=none offset=(0,0) length=&y_length in major=none minor=none value=none /*style=0*/;

axis2 color=gray label=none offset=(0,0) length=&x_length in major=none minor=none value=none /*style=0*/;

footnote link="http://maps.stamen.com/"
  c=gray ls=1.0 font='albany amt'
  "Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA.";

title;
proc gplot data=blank_map anno=anno_maps;
plot y*x / anno=anno_stuff 
 vaxis=axis1 haxis=axis2
 des='' name="&name";
run;





quit;
ODS HTML CLOSE;
ODS LISTING;