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


/* I started with an example in the book
   "The Encore Series, SAS/Graph Software", p. 141.
   It was originally from Observations, First Quarter 1994.
   
   I got the code from ...
   ftp ftp.sas.com
     anonymous
     email address
   cd observations
   cd picturet
   get picturet.sas

   This version of the sample didn't have all the data,
   so I went looking for it under the author's unix directory.
   I found it here ...

   Here's the original data:  /u/sasjzr/data/topo.dat
   Here's the original sas program:  /u/sasjzr/saspgms/picthis_issue10.sas

   Then, I modified it to use gmap to draw the map, and annotate to 
   draw the elevation colors 'behind' the map (whereas the original
   used annotate to draw the map outline on top of a gcontour of
   the elevation data).

   In this example, I'm using the entire US data, rather than just
   the subset used in the 'encore series' example.  Here's how I
   got the data ...

   http://www.ngdc.noaa.gov/mgg/gdas/gd_designagrid.html
   Your Grid Id: bobsas1
   Grid Database: ETOPO2 - 2-minute global relief
   x - Windows
   Lat/Lon (left defaults - lat(up/down) = 50N-25N , lon (across)=125W-67W
   Cell size - 10-minute  (the 2-minute is too big to handle on a pc)
   Cell params - default (2-byte integers, whole meters precision)
   Format: XYZ (lon,dat,depth), No Header, Space delimited, omit empty grid cells
 
   Clicked 'Design-A-Grid' button.
   Clicked 'Next-->' button
   Clicked 'Compress and Retrieve' button
   Clicked 'Retrieve' button
   Clicked 'save' 
   Saved bobsas1.zip to my U:
   Then, after saving it, I unzipped it, which created a directory structure...
   There's a lot of 'junk', but the file I wanted was ...

      /u/realliso/bobsas1_data/bobsas1/bobsas1.xyz
      (I renamed that to topo_10min.dat)
  
*/

/**********************************************************************/

filename elev 'topo_10min.dat';
data rawelev;
/* the CRLF allows dos sas to easily read the pc-format topo_10min.dat file */
infile elev dlm=',' termstr=CRLF;
   input long lat elev;
   x=-1*(long/57.295); 
   y=lat/57.295; 
   drop long lat;
run;

data annoelev; set rawelev;
   length function style color $ 8 position $ 1 text $ 30;
   retain xsys ysys '2' hsys '3' when 'b';
  anno_flag=1;
   function='pie';
   style='psolid';
   position='5';
   rotate=360;
   /* Make the size of the annotated colored pies big enough so that
   they overlap, and there is no white-space visible between them.
   If your elevation data is less dense, then you might want to annotate
   a rectangular areas rather than circles, so that they will meet-up
   better on the edges.  My elevation data is very dense, so we're almost
   down to the pixel-level, so I'm just using simple dot/pie's instead
   of rectangles.  (if I make the pie size smaller than .3, then I do
   get some white-space showing through... */
   size=.3;
if elev < -250 then color="cx0000ff";       /* deep ocean */
else if elev <= 0 then color="cyan";        /* below sea level */
else if elev <= 250 then color="cx005a32";  /* dark green - grassland */
else if elev <= 750 then color="cx238b45";
else if elev <= 1250 then color="cx41ab5d";
else if elev <= 1750 then color="cxfdae6b";
else if elev <= 2250 then color="cxfe9929";
else if elev <= 2750 then color="cxec7014";
else if elev <= 3250 then color="cxcc4c02";
else if elev > 3250 then color="cx8c2d04";  /* highest mountain - brown */
   output;
run;

/* Annotate a legend */
data color_legend;
input mytext $ 1-12 color $ 14-21 x y;
cards;
> 3,250      cx8c2d04 88 55
2,750-3,250  cxcc4c02 88 50
2,250-2,750  cxec7014 88 45
1,750-2,250  cxfe9929 88 40
1,250-1,750  cxfdae6b 88 35
7,50-1,250   cx41ab5d 88 30
250-750      cx238b45 88 25
0-250        cx005a32 88 20
-250 - 0     cyan     88 15
< -250       cx0000ff 88 10
;
run;

data color_legend; set color_legend;
   length function color $ 8 style $25 position $ 1 text $ 20;
   retain xsys ysys '3' hsys '3' when 'a';
   position='B';

   function='label'; size=4; text='U';
   /* color=color; */ style='marker';  output;
   color="black"; style='markere';  output;

   function='label'; color="black"; 
   text=mytext;
   x=x+7; y=y+.5; 
   style='"arial"'; size=2.75; output;
run;

data legend_title;
   length function color $ 8 style $25 position $ 1 text $ 20;
   retain xsys ysys '3' hsys '3' when 'a';
   function='label'; size=4; text='U';
   position='B';
   style='"arial"'; size=3.5; x=88+7-2;
   text='Elevation';   y=64.5; output;
   text='(in meters)'; y=60.5; output;
run;

data color_legend; set color_legend legend_title;
run;

data legend_box; 
   length function style color $ 8;
   retain xsys ysys '3' hsys '3' when 'b';
   color='white';
   style='msolid';
   function='poly'; 
   x=85; y=0; output;
   function='polycont'; 
   x=100; y=0; output;
   x=100; y=88; output;
   x=85; y=88; output;
run;
data color_legend; set legend_box color_legend;
run;


data mymap; set maps.states 
 (
  where=
    (
      (fipstate(state) not in ('AK' 'HI' 'PR'))
      and
      (density < 3)
    ) 
 );
run;

data combined; set mymap annoelev; run;
proc gproject data=combined out=combined;
id state;
run;
data mymap annoelev;
 set combined;
  if anno_flag=1 then output annoelev;
  else output mymap;
run;

proc sql;
create table mydata as
select unique state 
from mymap;
quit; run;
data mydata; set mydata;
length htmlvar $200.;
length altvar $100.;
length drillvar $100.;
/*Create the drilldown url */
 drillvar='HREF="http://www.ngdc.noaa.gov/cgi-bin/seg/topo/state2.pl?region='||
          trim(left(lowcase(fipstate(state))))||'.jpg"';
/*Create tooltips that display when you hover over a state. */
/*These tooltips will only appear when the output is viewed in */
/*Internet Explorer, not Netscape.  This is expected behavior. */
/*The 'OD'x hex value creates a line break in the text string */
 altvar = 'title="State FIPs code: ' || trim(left(state)) || '0D'x ||
           'Abbreviation: '||trim(left(fipstate(state)))||'0D'x||
           'State name: ' || trim(left(fipnamel(state))) || '"';
/*Combine the tooltip and drill down into one variable */
 htmlvar = drillvar || ' ' || altvar;
run;

 

/**********************************************************************/

 GOPTIONS DEVICE=gif;

 ODS LISTING CLOSE;
 ODS HTML path=odsout body="&name..htm"
 style=minimal
 ;
goptions border;

goptions xpixels=800 ypixels=515;

goptions gunit=pct ftitle="arial/bo" htitle=6.1 ftext="arial" htext=3.5;
title1 'US Elevation Map';
title2 'Click on states to see detailed NGDC maps';
/* blank title along the right-hand side, to make room for annotated legend/key */
title3 a=-90 h=24 ' ';

/* Make the map areas 'empty/see-through' so the annotated elevation 
   colors that are behind them will show through. */
pattern v=e c=black repeat=100;

proc gmap data=mydata map=mymap anno=annoelev;
 id state;
 choro state /
 nolegend
 anno=color_legend
 html=htmlvar
 des="" name="&name";
run;
quit;

 ODS HTML CLOSE;
 ODS LISTING;