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

/*
Got this map from here:
http://water.weather.gov/precip/download.php

Shapefile - current data
Last 7 days
Observed
Copy & extract to D:\public\precip\

Using color gradient from:
http://2.bp.blogspot.com/-SI2qFlYpUEA/UR5mQAZhUiI/AAAAAAAAEAo/Syfa06ca2qI/s1600/Western+US+60+days.gif
*/

%let folder=d:\public\precip\;
libname robsdata "&folder";

/* 
globvalue = inches of rain 
There is 1 data point for each lat/lon, containing the sum of 7 days.
*/
proc mapimport datafile="&folder\nws_precip_last7days_observed_20151005.shp" out=precip_point_data; 
id id;
run; 

%let latmin=32.0;
%let latmax=39.4;
%let longmin=-84.8;
%let longmax=-74.0;

data anno_dots; set precip_point_data 
 (rename=(lon=long) drop=id
 where=(
 lat>&latmin and
 lat<&latmax and
 long>&longmin and
 long<&longmax
 ));
flag=1;
run;

data my_map; set mapsgfk.us_counties (where=(statecode not in ('AK' 'HI') and density<=2) drop=resolution);
run;

data combined; set my_map anno_dots;
run;
proc gproject data=combined out=combined latlong eastlong degrees dupok
 latmin=&latmin latmax=&latmax
 longmin=&longmin longmax=&longmax
 ;
id statecode county;
run;
data my_map anno_dots; set combined;
if flag=1 then output anno_dots;
else output my_map;
run;

/* the breakpoints between the legend colors */
%let leg1=1;
%let leg2=2;
%let leg3=3;
%let leg4=4;
%let leg5=5;
%let leg6=6;
%let leg7=7;
%let leg8=8;
%let leg9=9;
%let leg10=10;

/* the legend colors */
%let col1=cxbe0000;
%let col2=cxfa0000;
%let col3=cxf08229;
%let col4=cxe6b02e;
%let col5=cxe6e600;
%let col6=cx64c864;
%let col7=cx00aa00;
%let col8=cx00a1e6;
%let col9=cx6e00dc;
%let col10=cxa100c7;
%let col11=cxfa00fa;

proc sort data=anno_dots out=anno_dots;
by Globvalue;
run;

data anno_dots; set anno_dots;
length function $8 color $8 style $35;
xsys='2'; ysys='2'; hsys='3'; 
when='a';
/* you can use when='b' if you want the county outlines to show up */
function='pie'; rotate=360; size=0.25; style='psolid'; 
     if Globvalue<&leg1 then color="&col1"; 
else if Globvalue<&leg2 then color="&col2"; 
else if Globvalue<&leg3 then color="&col3"; 
else if Globvalue<&leg4 then color="&col4"; 
else if Globvalue<&leg5 then color="&col5"; 
else if Globvalue<&leg6 then color="&col6"; 
else if Globvalue<&leg7 then color="&col7"; 
else if Globvalue<&leg8 then color="&col8"; 
else if Globvalue<&leg9 then color="&col9"; 
else if Globvalue<&leg10 then color="&col10"; 
else color="&col11";
output;
run;

data anno_legend;
length inches $10;
number=1; inches=' '; output;
number=2; inches="&leg1"; output;
number=3; inches="&leg2"; output;
number=4; inches="&leg3"; output;
number=5; inches="&leg4"; output;
number=6; inches="&leg5"; output;
number=7; inches="&leg6"; output;
number=8; inches="&leg7"; output;
number=9; inches="&leg8"; output;
number=10; inches="&leg9"; output;
number=11; inches="&leg10"; output;
run;

data anno_legend; set anno_legend;
length function color $8;
xsys='3'; ysys='3'; hsys='3'; when='a';

/* text labels below the legend segments */
function='label'; position='e'; size=2.5;
x=17.5+number*5; y=90;
text=inches;
output;
size=.001;

/* colored segments */

/* bottom/left of each color chicklet */
x=17.5+number*5; y=90;
function='move'; 
output;

/* top/right of each color chicklet */
x=x+5; y=y+2.5;
function='bar'; style='solid'; line=0;
color=symget('col'||trim(left(number)));
output;

/* dark outline around colored segments */

x=17.5+number*5; y=90;
function='move'; 
output;

x=x+5; y=y+2.5;
function='bar'; style='empty'; line=0;
color="gray44";
output;

run;

data my_data; set mapsgfk.us_counties_attr;
length my_html $200;
my_html='title='||quote(trim(left(idname))||' county, '||trim(left(statecode)));
run;

/* annotate state outlines on the county map */
data my_map; set my_map;
original_order=_n_;
run;
proc sort data=my_map out=my_map;
by state original_order;
run;
proc gremove data=my_map out=anno_outline;
id state county;
by state;
run;
data anno_outline; set anno_outline;
by state segment notsorted;
xsys='2'; ysys='2'; hsys='3'; when='a';
length function $8;
color='gray44'; style='empty'; size=2.0;
if first.state or first.segment then function='poly';
else function='polycont';
run;


goptions device=png;
goptions xpixels=750 ypixels=825;
goptions border;
 
ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="1-Week Precipitation map") 
 style=htmlblue;

goptions gunit=pct htitle=4.0 ftitle="albany amt/bold" htext=2.25 ftext="albany amt";

title1 
 link="http://water.weather.gov/precip/download.php"
 ls=1.5 "Inches of Precipitation";

title2 h=3 ' ';

data anno_stuff;
 set anno_dots anno_legend anno_outline;
run;

pattern1 v=empty c=white repeat=5000;

proc gmap data=my_data map=my_map anno=anno_stuff;
note move=(60,4) "Sep 29, 2015 - Oct 5, 2015";
id statecode county;
choro statecode / nolegend
 coutline=gray99 
 html=my_html
 des='' name="&name";
run;

quit;
ODS HTML CLOSE;
ODS LISTING;