goptions reset=all target=xl300ps ftext=swissb
rotate=landscape lfactor=2;


%let longmax=79;    /* western longitude boundary for map */
%let latmin=33.5;   /* southern latitude boundary for map */
%let latmax=36.75;  /* northern latitude boundary for map */
%let d2r=atan(1) / 45; /* algorithm for converting longitude
                          and latitude values from units of
                          degrees to Cartesian coordinate
                          values in units of radians */


data lights;
length name $17;
input name $ & range latdeg latmin longdeg longmin;
x=&d2r * (longdeg + longmin/60);
y=&d2r * (latdeg + latmin/60);
cards;
Currituck Beach  18 36 22.6 75 49.8
Bodie Island  18 35 49.2 75 33.8
Cape Hatteras  24 35 15.3 75 31.2
Ocracoke  14 35 06.5 75 59.2
Cape Lookout  25 34 37.3 76 31.5
Oak Island  24 33 53.6 78 02.1
;
run;


data ncscva;
set maps.counties;
if state in(37 45 51);
run;

data all;
set lights (in=lights) ncscva;
if lights then state=199;
run;

proc gproject data=all out=allproj
              longmax=&longmax
              latmin=&latmin
              latmax=&latmax
              ;
id state county;
run;


data mapproj annoproj;
set allproj;
if state=199 then output annoproj;
else output mapproj;
run;


data anno;
set annoproj;
length function color style $8 text $17;
maprange=(&latmax*&d2r)-(&latmin*&d2r);
mapdist=3949.99*maprange;
scale=maprange/mapdist;
piesize=((range*scale)/maprange) * 100;
retain hsys '1' lightx lighty;
lightx=x; lighty=y;

/* Place symbol on map at lighthouse location */
xsys='2'; ysys='2'; when='a'; color='red';
function='label'; style='special'; text='J'; size=3; position='+'; output;

/* Place lighthouse name and visibility range next to symbol */
ysys='9'; xsys='2';
if (latdeg gt 35) or (latdeg = 35 and latmin ge 30) then y=5;
else if name='Cape Hatteras' then y=-2;
else y=-3;
function='label'; style='swissb'; position='>'; size=2.5; color='black';
text=trim(left(name));
output;

ysys='9';  y=-1.5;
position='>'; text=trim(left(range))||' miles';  output;

/* Draw visibility ranges on the map */
function='pie'; hsys='1'; ysys='2'; xsys='2'; style='solid'; color='yellow';
angle=90; rotate=360; x=lightx; y=lighty; size=piesize;
when='b'; output;
when='a'; style='empty'; color='black'; output;
run;

pattern1 v=solid c=greyee;
pattern2 v=empty c=black;
pattern3 v=empty c=black;

title1 f=swissb h=3 'Selected Lighthouses on the North Carolina Coast:';
title2 h=1 ' ';
title3 f=swissb h=3 'How far can you see the light?';
footnote1 j=r h=3 pct ' ';


proc gmap map=mapproj data=mapproj anno=anno;
id state county;
choro state / discrete nolegend coutline=black;
run;
quit;