/* NAME 
/*
/* TOPOAREA.AML
/*
/*----------------------------------------------------------------------------
/*
/* AUTHOR Jacek Blaszczynski, Physical Scientist, BLM NARSC, 10/99, 
/* Original Coding
/*
/*----------------------------------------------------------------------------
/*
/* DESCRIPTION:
/*
/* This routine estimates the topographically correct area values for
/* particular classes or geographic zones within your project area
/* utilizing slope gradient map generated from a Digital Elevation 
/* Model (DEM) for the project area. The theory behind it is that the
/* slope gradient for each particular cell can be understood as directly
/* proportional to the elongation of the cell in one direction.  This 
/* elongation can be can calculate using trigonometric equation:
/*
/*   h (elongated side) = x (current cell size) / cosine x
/*
/* and the new area of the cell which accounts for it's slope can be
/* calculated as:
/* 
/*   True area estimate per cell (At) = x * h
/*
/* while the existing area estimate per cell (no topo correction) = x * x
/*
/* For any particular zone we can use the ARC GRID ZONALSUM function
/* to add the areas for each cell of a particular class or zone. 
/* 
/* USAGE: &R TOPOAREA1 <SOURCE DEM GRID> <CLASS GRID> <VALUE FOR CELL SIZE>
/*
/*------------------------------------------------------------------------------
/*
/* Default map names:
/*
/* slopedeg (slope gradient map in degrees), cosine%cell_size% (cosine x map),
/* topoelong%cell_size% (map with values for length of elongated cell side),
/* topoarea%cell_size% (map with values for topo corrected cell area)
/* areameter%cell_size% (map with values for uncorrected cell areas x*x)
/* arearatio%cell_size% (map of ratios of corrected to uncorrected areas/cell)
/* classtopo%cell_size% (map of topographically corrected areas per class)
/* classacre%cell_size% (map of uncorrected areas per class using ZONALSUM)
/* classratio%cell_size% (map of ratios of corrected to uncorrected areas/zone)
/* areadiff%cell_size%(map of differences in corrected - uncorrected areas/zone)
/*
/*-----------------------------------------------------------------------------------
/*
/* WARNING:
/*
/* THIS PROGRAM ASSUMES THAT YOU ARE USING ONLY ONE DIGITAL ELEVATION MAP
/* PER PROJECT AREA WITHIN A SINGLE WORKSPACE, AND ANY SUBSETS OF THE PROJECT
/* AREA ARE PLACED IN SEPARATE WORKSPACES.  BASED ON THIS ASSUMPTION THE AML
/* CREATES A NUMBER OF DEFAULT GRIDS WHICH WILL EITHER BE REUSED BY ANOTHER RUN
/* OF THE AML IN THE FUTURE, OR RECREATED AGAIN.  IT IS ESSENTIAL THAT THE EXTENT
/* OF THE DIGITAL ELEVATION MODEL (DEM) IS THE LARGEST AREA EXTENT THAT IS BEING
/* ANALYZED WITHIN A WORKSPACE.
/*
/*-----------------------------------------------------------------------------------

&severity &error &routine bailout
&args source_dem class_grid cell_size
/* Check arguments
&call chkargs

/* Do the work
&call main

&call exit
&return /* End of topoarea1.aml

/*--------------
&routine chkargs
/*--------------
&if [SHOW PROGRAM] <> GRID &then &do
  &type This aml must be run from GRID.
  &call bailout
&end
&do arg &list source_dem class_grid cell_size
  &if [NULL [VALUE %arg%]] &then &do
    &call usage
    &call bailout
  &end
&end
&if not [EXISTS %source_dem% -GRID] &then &do
  &type Input grid: %source_dem%, does not exist.
  &call bailout
&end
&return /* End of routine chkarg


/*------------------------------------------------------------------------------
&routine main
/*------------------------------------------------------------------------------
/*
/* Set the environment
/*
/*---------------------------------------------------------------------
&terminal 9999
display 9999 3
setwindow %source_dem% %source_dem%
setcell %source_dem%
mape %source_dem%
gridnodatasymbol transparent
gridpaint %source_dem% value linear wrap gray

/*---------------------------------------------------------------------
/* Clean the class topographically corrected area calculations
/* from previous runs for the same cell size and prepare default
/* data preparation grids if not yet made and in workspace
/* CAUTION: There should only be one set of default data prep grids
/* per project area and particular cell size - these grids remain
/* in the workspace until manually deleted by the user.
/*---------------------------------------------------------------------
&if not [EXISTS slopedeg -GRID] &then &do 
  SLOPEDEG = slope ( %source_dem%, degree )
&end
&if [EXISTS cosine%cell_size% -GRID] &then &do 
  kill cosine%cell_size% all
&end
&if [EXISTS areameter%cell_size% -GRID] &then &do
  kill areameter%cell_size% all
&end
&if [EXISTS topoelong%cell_size% -GRID] &then &do
  kill topoelong%cell_size% all
&end
&if [EXISTS topoarea%cell_size% -GRID] &then &do
  kill topoarea%cell_size% all
&end
&if [EXISTS arearatio%cell_size% -GRID] &then &do
  kill arearatio%cell_size% all
&end
&if [EXISTS classtopo%cell_size% -GRID] &then &do
  kill classtopo%cell_size% all
&end
&if [EXISTS classacre%cell_size% -GRID] &then &do
  kill classacre%cell_size% all
&end
&if [EXISTS classratio%cell_size% -GRID] &then &do
  kill classratio%cell_size% all
&end
&if [EXISTS topoarea%cell_size% -GRID] &then &do
  kill topoarea%cell_size% all
&end
&if [EXISTS areadiff%cell_size% -GRID] &then &do
  kill areadiff%cell_size% all
&end

/*-------------------------------------------------------------------
/* Develop cosine map
/*-------------------------------------------------------------------

COSINE%cell_size% = cos (slopedeg div deg)
gridpaint cosine%cell_size% value linear nowrap gray

/*-------------------------------------------------------------------
/* Develop uncorrected area per cell map
/*-------------------------------------------------------------------

AREAMETER%cell_size% = %cell_size% * %cell_size%
gridpaint areameter%cell_size% value linear nowrap gray

/*-------------------------------------------------------------------
/* Develop a map of elongated cell lengths
/*------------------------------------------------------------------- 

TOPOELONG%cell_size% = ( %cell_size% / cosine%cell_size% )

/*-------------------------------------------------------------------
/*  Develop a map of topographically corrected cell areas
/*-------------------------------------------------------------------

TOPOAREA%cell_size% = ( %cell_size% * topoelong%cell_size% )
gridpaint topoarea%cell_size% value linear nowrap gray
kill topoelong%cell_size% all

/*-------------------------------------------------------------------
/*  Develop a map of ratios of corrected to uncorrected area per cell
/*-------------------------------------------------------------------

AREARATIO%cell_size% = topoarea%cell_size% / areameter%cell_size%
gridpaint arearatio%cell_size% value linear nowrap gray

/*--------------------------------------------------------------------------
/*  Develop a map of topographically corrected areas per class or zone
/*--------------------------------------------------------------------------

CLASSTOPO%cell_size% = ((zonalsum (%class_grid%, topoarea%cell_size%, data)) / 4046.852)

/*-----------------------------------------------------------------------------
/*  Develop a map of uncorrected areas per class or zone using ZONALSUM
/*-----------------------------------------------------------------------------

CLASSACRE%cell_size% = (( zonalsum (%class_grid%, areameter%cell_size%, data)) / 4046.852)

/*------------------------------------------------------------------------------
/*  Develop a map of ratios between topographically corrected and
/*  uncorrected area values of each class or zone
/*------------------------------------------------------------------------------

CLASSRATIO%cell_size% = classtopo%cell_size% / classacre%cell_size%
gridpaint classratio%cell_size% value linear nowrap gray

/*------------------------------------------------------------------------------
/*  Develop a map of differences between areas for each class or zone
/*------------------------------------------------------------------------------

AREADIFF%cell_size% = classtopo%cell_size% - classacre%cell_size%
gridpaint areadiff%cell_size% value linear nowrap gray
&return /* end of routine main



/*------------
&routine usage
/*------------
&type Usage: &r topoarea1 <source_dem> <class_grid> <cell_size>>
&return /* End of routine usage

/*--------------
&routine exit
/*---------------
&return /* End of routine exit

/*--------------
&routine bailout
/*--------------
&severity &error &ignore
&call exit
&return &error Bailing out of topoarea1.aml...