/* demsr.aml
/* macro to make a color shaded-relief image
/* with light cyan water and even lighter glaciers  
/*    Ralph Haugerud, U.S. Geological Survey
/*    April 1998
/*    haugerud@usgs.gov
/* to be run from ARC
/*

&args DEM WaterCov GlacierCov

/* DEM is input DEM, defines extent and resolution of output.
/* GlacierCov and WaterCov are coverages with closed polys around glaciers
/* and open water bodies, respectively, and NO BORDER. Not-glacier and 
/* not-water should be parts of the universe polygons.

&s azimuth = 315
&s inclination = 45
&describe %DEM%
&set cellsize = %GRD$DX%

&set OldDisplay = [SHOW DISPLAY]
DISPLAY 0
&call MakeVGrid
&call MakeHGrid
&call MakeSGrid
&if [exists %GlacierCov% -cover] &then &call glaciers
&if [exists %WaterCov% -cover] &then &call water

DISPLAY %OldDisplay%
&return
/**************************************
&routine MakeVGrid  
  GRID
  v1 = hillshade(%DEM%, 280, 45, shade)
  &describe v1
  v2 =  (v1 - [VALUE GRD$MEAN]) / [value GRD$STDV] * 15 + 85
  %DEM%v = con(v2 < 60, 60, v2 > 99, 99, int(v2))
  QUIT
  kill v1 all
  kill v2 all
  &return
  
&routine MakeHGrid
  GRID
  %DEM%h = con(%DEM% < 4000, 120 - %DEM% * 3 / 100) 
  /* for DEM in meters, this maps 1 to 4000 meters onto 
  /* Green (hue 120) to red (hue 0)
   QUIT
   &return
   
&routine MakeSGrid
  GRID
  %DEM%s = con(%DEM%h > 80, 40, INT(%DEM%h / 2) )
  QUIT  
  &return
  
&routine Glaciers
  polygrid %GlacierCov% GlacGrid
  %CellSize%    /* set cell size to match DEM
  Y             /* convert entire coverage
  GRID
  GlacierSat = con(not isnull(GlacGrid), 12)
  GlacierHue = con(not isnull(GlacGrid), 180)
  %DEM%s2 = merge(GlacierSat, %DEM%s)
  kill %DEM%s all
  rename %DEM%s2 %DEM%s
  %DEM%h2 = merge(GlacierHue, %DEM%h)
  kill %DEM%h all
  rename %DEM%h2 %DEM%h
  QUIT
  kill GlacGrid all
  &return
  
&routine Water
  polygrid %WaterCov% WatGrid
  %CellSize%    /* set cell size to match DEM
  Y             /* convert entire coverage
  GRID
  WaterSat = con(not isnull(WatGrid), 30)
  WaterHue = con(not isnull(WatGrid), 180)
  %DEM%s2 = merge(WaterSat, %DEM%s)
  kill %DEM%s all
  rename %DEM%s2 %DEM%s
  %DEM%h2 = merge(WaterHue, %DEM%h)
  kill %DEM%h all
  rename %DEM%h2 %DEM%h
  QUIT
  kill WatGrid all
  &return
  
&return