/*  Program:  dem2hillsh.aml
/*
/*  Purpose:  Creates ARC elevation grid and hillshade grid from a DEM file.
/*
/*  Usage:    &r dem2hillsh <in_dem> <out_elevgrid> <out_hillshgrd> <hillsh_res>
/*                          <azimuth> <altitude> <vert_exagg> 
/*                          <smoothing_kernel_size> <angle> {clip_cover}
/*                          {datum} {projection} {horiz_units} {UTM_Zone}
/*     Ex:    &r dem2hillsh gol.dem elev hillshade 10 315 35 5 0 1.682 goln ~
/*                        nad27 UTM meters 13
/*            &r dem2hillsh mor.dem morz hillsh 30 310 40 4 10 -0.65 # NAD83 utm
/*
/*  Output: Outputs are 1) a raw elevation grid of the same resolution as the
/*          input DEM, and 2) a hillshaded elevation grid with cellsize 
/*          of <out_res>, clipped to the extent of <clip_cover>, rotated an 
/*          angle of <angle> (positive angle = clockwise, negative = counter),
/*          with shaded relief parameters: altitude of <alt>, azimuth of <az>,
/*          and vertical exaggeration of <ve>.
/*
/*  Notes:  Run this program without arguments to see the usage and examples.
/*
/*          The <smoothing_kernel_size> is the size, in cells, of the square
/*          filter which will be passed over the hillshade grid (use 0 for no
/*          smoothing, use 10 for significant smoothing, any integer okay).
/*          If an angle of 0 is entered, the hillshade grid will not be rotated.
/*          Use the normal "#" character to skip optional arguments.
/*          If no UTM zone is passed in, no zone information will exist
/*          in the projection file.
/*
/*  History:  JKosovich  03/05/97  Coding for abandoned mine lands, based on 
/*                                 AMLs image2grid, makeshade, drg2grid, etc.
/*            JKosovich  09/10/97  Renamed from old dem2grid.aml
/*            JKosovich  10/10/97  Cleans up leftover Grid trash (xyxy001...)
/*            JKosovich  01/13/98  Turned watch off, loop to clean up grid junk
/*            JKosovich  02/04/98  Added <smoothing_kernel_size> argument
/*
/*******************************************************************************
/*
&args indem outelevgrd outshgrid outshres az alt ve xy angle clipcov datum ~
  proj units zone
&severity &error &routine bailout
&s prog = [ENTRYNAME %aml$file% -NOEXT -FILE]     /* current AML being run
&watch %prog%.watch

/*
/***** Check for correct usage and correct argument types                  *****
/*
&if [NULL %angle%] &then
  &do
    &s error_case = 1
    &call error_cond
  &end

&if [LENGTH %outelevgrd%] > 12 | [LENGTH %outshgrid%] > 12 &then
  &do
    &s error_case = 2
    &call error_cond
  &end

&do xx &list outshres angle ve                      /* these can be int or real
  &if [TYPE [VALUE %xx%]] > 0 &then
    &do
      &s error_case = 3
      &call error_cond
    &end
&end

&do xx &list alt az xy zone                         /* these can only be integer
  &if ^ [NULL [VALUE %xx%]] &then
    &do
      &if [TYPE [VALUE %xx%]] ^= -1 &then
        &do
          &s error_case = 4
          &call error_cond
        &end
    &end
&end

/*
/***** Set toggle options                                                  *****
/*
&s smooth = .TRUE.                                    /* Don't smooth if 0 xy
&if %xy% = 0 &then &s smooth = .FALSE.

&s rotate = .TRUE.                                    /* Don't rotate if 0 ang
&if %angle% = 0 &then &s rotate = .FALSE.
&if %angle% = 0.0 &then &s rotate = .FALSE.

&do xx &list clipcov datum proj units zone            /* if #, set to null
  &if ^ [NULL [VALUE %xx%]] &then
    &do
      &if [TYPE [VALUE %xx%]] = 1 &then               /* if character
        &do
          &if [VALUE %xx%] = # &then
            &s %xx% =  
        &end
    &end
&end

&s projdef = .TRUE.
&if [NULL %proj%] | [NULL %datum%] &then
  &s projdef = .FALSE.

&s clip = .TRUE.
&if [NULL %clipcov%] &then
  &s clip = .FALSE.

/*
/***** Check for existence of input file and cover                         *****
/*
&if ^ [EXISTS %indem% -FILE] &then
  &do
    &s error_case = 5
    &call error_cond
  &end

&if %clip% &then
  &do
    &if ^ [EXISTS %clipcov% -COVER] &then
      &do
        &s error_case = 6
        &call error_cond
      &end
    &else &if ^ [EXISTS %clipcov% -POLYGON] &then
      &do
        &s error_case = 7
        &call error_cond
      &end
  &end

/*
/***** Set terminal                                                        *****
/*
&if [NULL [EXTRACT 1 [SHOW &terminal]]] &then
  &terminal 9999

/* 
/***** Convert DEM to Arc grid and convert to integer if necessary         *****
/*
&if [EXISTS %outelevgrd% -GRID] &then KILL %outelevgrd% ALL
&if [EXISTS xxgrd1 -GRID] &then KILL xxgrd1 ALL
DEMLATTICE %indem% %outelevgrd% USGS
&describe %outelevgrd%
&s indemres = %grd$dx%     /* resolution of incoming DEM = latt res
&s gridtype = %grd$type%   /* type of grid: 1 = int, 2 = floating pt.
&if %gridtype% = 2 &then
  &do
    DISPLAY 0
    GRID
    xxgrd1 = INT(%outelevgrd%)   /* Convert to integer, much smaller file size!
    QUIT
  &end

&if [EXISTS xxgrd1 -GRID] &then
  &do
    KILL %outelevgrd% ALL
    RENAME xxgrd1 %outelevgrd%
  &end

/*
/***** Define projection for grid                                          *****
/*
&if %projdef% &then
  &do
    PROJECTDEFINE GRID %outelevgrd%
      PROJECTION %proj%
      DATUM %datum%
      &if ^ [NULL %units%] &then
        UNITS %units%
      &if ^ [NULL %zone%] &then
        ZONE %zone%
      PARAMETERS
      &if [UPCASE [SHOW program]] = PRJDEF &then END
  &end

/* 
/***** Resample grid                                                       *****
/*
&if [EXISTS xxgrd2 -GRID] &then KILL xxgrd2 ALL
&if %indemres% ^= %outshres% &then
  &do
    DISPLAY 0
    GRID
    xxgrd2 = RESAMPLE(%outelevgrd%,%outshres%)
    QUIT
  &end

&if [EXISTS xxgrd2 -GRID] &then 
  PROJECTCOPY GRID %outelevgrd% GRID xxgrd2
&else
  COPY %outelevgrd% xxgrd2

/* 
/***** Clip grid to %clipcov% extent                                       *****
/*
&if [EXISTS xxgrd3 -GRID] &then KILL xxgrd3 ALL
&if %clip% &then
  LATTICECLIP xxgrd2 %clipcov% xxgrd3 MINIMUM 1
&if [EXISTS xxgrd3 -GRID] &then
  &do
    PROJECTCOPY GRID xxgrd2 GRID xxgrd3
    KILL xxgrd2 ALL
  &end
&else
  RENAME xxgrd2 xxgrd3

/*
/***** Rotate grid to due N-S orientation for plotting                     *****
/*
&if [EXISTS xxgrd4 -GRID] &then KILL xxgrd4 ALL
&if %rotate% &then
  GRIDROTATE xxgrd3 xxgrd4 %angle%
&if [EXISTS xxgrd4 -GRID] &then
  &do
    PROJECTCOPY GRID xxgrd3 GRID xxgrd4
    KILL xxgrd3 ALL
  &end
&else
  RENAME xxgrd3 xxgrd4

/*
/***** Apply HILLSHADE parameters                                          *****
/*
&if [EXISTS %outshgrid% -GRID] &then KILL %outshgrid% ALL
DISPLAY 0
GRID
&if %smooth% &then
  %outshgrid% = HILLSHADE(FOCALMEAN(xxgrd4,RECTANGLE,%xy%,%xy%),%az%,%alt%,SHADE,%ve%)
&else
  %outshgrid% = HILLSHADE(xxgrd4,%az%,%alt%,SHADE,%ve%)
QUIT

&if [EXISTS %outshgrid% -GRID] &then
  &do
    PROJECTCOPY GRID xxgrd4 GRID %outshgrid%
    KILL xxgrd4 ALL
    &type
    &type
    &type *** %prog% successful. Shaded relief grid is %outshgrid%
    &type
    &type
  &end
&else
  &do
    &type
    &type
    &type *** %prog% ERROR: Shaded relief grid %outshgrid% not created!!!
    &type *** Last temporary grid made was: xxgrd4.
    &type
    &type
  &end

/*
/***** Clean up Grid junk                                                  *****
/*
&messages &off &all
&do n = 1 &to 9
  &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL
&end
&messages &on
  
/*
/*
&watch &off
&return
/*
/*
/*------------------------------------------------------------------------------
/*------------------------------  ROUTINE LIST  --------------------------------
/*------------------------------------------------------------------------------
/*
/*
/*------------------
&routine ERROR_COND
/*------------------
/*
&severity &error &routine bailout
&if [UPCASE [SHOW PROGRAM]] = TABLES &then QUIT
&messages &on

&type
&type ***************  PROGRAM [TRANSLATE %prog%] HALTED  *****************
/*
&if %error_case% = 1 &then  /***** Incorrect usage
  &do
    &type USAGE: &r %prog% <in_dem> <out_elevgrid> <out_hillshgrid> <hillsh_res>
    &type           <azimuth> <altitude> <vert_exagg> <smoothing_kernel_size>
    &type           <rotation_angle> {clip_cover}
    &type           {datum} {projection} {horiz_units} {UTM_Zone}
    &type
    &type    Ex: &r %prog% gol.dem elev hillsh 10 315 35 5 0 1.682 goln 
    &type                      nad27 UTM meters 13
    &type        &r %prog% abc.dem abcelv hillsh 30 310 40 4 9 -0.65 # NAD83 utm
    &type
    &type OUTPUT ELEV and HILLSHADE GRID NAMES must be 12 characters or less.
    &type OUTPUT HILLSHADE RESOLUTION is in ground units (eg. meters).
    &type AZIMUTH is whole degrees clockwise from North, North being 0.
    &type ALTITUDE is whole degrees above the horizon, horizon being 0.
    &type VERTICAL EXAGGERATION can be integer or whole, typically 5 or less.
    &type SMOOTHING KERNEL: use 0 for no hillshade smoothing, any integer OK.
    &type Use positive ANGLE to rotate grid clockwise to grid north,
    &type Use negative ANGLE to rotate grid counterclockwise.
    &type Use ANGLE of 0 for no rotation.
    &type If used, CLIP COVER must have POLYGON topology.
    &type If used, DATUM is the input data's datum (upper/lower case is okay).
    &type If used, PROJECTION is the input data projection.
    &type If used, UNITS is the the input data horizontal units.
    &type If used, ZONE is the input data projection zone.
  &end
&if %error_case% = 2 &then  /***** If output grid name too long
  &type          Output grid name is too long - must be 12 characters or less!!!
&if %error_case% = 3 &then  /***** Argument must be int or real
  &type          Argument %xx% must be an integer or real value!!!
&if %error_case% = 4 &then  /***** Argument must be integer
  &type          Argument %xx% must be an integer value!!!
&if %error_case% = 5 &then  /***** If dem file not found
  &type          DEM file %indem% does not exist!!!
&if %error_case% = 6 &then  /***** If clip cover not found
  &type          Clip coverage %clipcov% does not exist!!!
&if %error_case% = 7 &then  /***** If clip cover not a POLY cover
  &type          Clip coverage %clipcov% does not have POLYGON topology!!!
/*
&type      Process halted [DATE -DOW] [DATE -CAL] [DATE -AMPM]
&type *********************************************************
&type
/*

/*
/***** Clean up Grid junk                                                  *****
/*
&messages &off &all
&do n = 1 &to 9
  &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL
&end
&messages &on

&watch &off
&return; &return &error
/*
/*
/*------------------
&routine BAILOUT
/*------------------
/*
&severity &error &ignore

/*
/***** Clean up Grid junk                                                  *****
/*
&messages &off &all
&do n = 1 &to 9
  &if [EXISTS xyxy00%n% -GRID] &then KILL xyxy00%n% ALL
&end
&messages &on

&type
&type *** Error occurred at line %aml$errorline% in %aml$errorfile%
&watch &off
&return; &return &error
/*
/*