; $Id: aspect.pro,v 5.4 2001/04/4
;
; Copyright (c) 2001 Javier G Corripio.
;
;+
; NAME:
;   aspect
;
; PURPOSE:
;   Returns a m columns x n rows array with aspects values for every grid cell
;
; CALLING SEQUENCE:
;   Result = aspect(DEM [,SLOPE=variable])
;
; INPUTS:
;   DEM : input variable array(m,n).  Geotiff format
;
; KEYWORD PARAMETERS:
;   SLOPE: Returns a named variable containing the slopes of every cell
;
;
; OUTPUTS:
;   result: array with the values of
;    every grid cell aspect in degrees.
;
;   slope:  array with the slope value
;    in degrees of every grid cell in the DEM
;
;
; EXAMPLE
;   asp=ASPECT(gradientz, SLOPES=slp)
;
; MODIFICATION HISTORY:
;   Written by:  Javier Corripio, April 2001.
;
;-

@cell_gradient

FUNCTION aspect, DEM, dl, SLOPE=slope

on_error,2  ;Return to caller if error
nparams=N_PARAMS()
IF (nparams EQ 0) THEN BEGIN
MESSAGE, 'USAGE: Result = ASPECT(DEM, Dl, [,SLOPE=variable])',/INFO
    return,!VALUES.F_NaN
ENDIF

sz=size(DEM,/dimensions)
cols=sz[0]  ;** number of columns
rows=sz[1]  ;** number of rows

aspecta = findgen(cols,rows)
aspecto = findgen(cols,rows)
xcoordinate = findgen(cols,rows)
ycoordinate = findgen(cols,rows)

;*** The next lines prevent arithmetic errors when x and y are zero
;*** in a=TAN(y/x).  It makes x=1/10^30 so resulting aspect negligible
;*** from original, but avoid singularities in function TAN

Grad = cell_gradient(dem, dl)
xcoordinate(*,*) = Grad(0,*,*)
ycoordinate(*,*) = Grad(1,*,*)
xzero = WHERE(xcoordinate EQ 0, count)
If count NE 0 THEN xcoordinate(xzero) = 1e-20

aspecta(*,*) = (ATAN(ycoordinate, xcoordinate)*!RADEG) + 90.0

aspecto(*,*) = aspecta(*,*)
negs = WHERE(aspecta LT 0, COUNT)
IF COUNT NE 0 THEN aspecto(negs) = (360.0 + aspecta(negs))

IF ARG_PRESENT(slope) THEN BEGIN
  zcoordinate = findgen(cols,rows)*0.0
  slope  = findgen(cols,rows)*0.0
  zcoordinate(*,*) = grad(2,*,*)
  slope = (ACOS(zcoordinate))*!RADEG
ENDIF

RETURN, aspecto

END