Aspect map..

[William A Perkins]

>>>>> "Tarek" == Tarek Sadek <sadek at aqua.civag.unimelb.EDU.AU> writes:

    Tarek> Hi Grassers,
    Tarek> I am looking for a code which can determine the aspect map
    Tarek> from a DEM in degrees azimuth (measured every 1
    Tarek> degree). All programs avialble in the current version
    Tarek> calculate the aspect of the cells every 45 degrees
    Tarek> (i.e. 0,45,90,.....,360 degrees) such as r.watershed,
    Tarek> r.slope.aspect,..etc.  Also I can't find r.water.aspect
    Tarek> which I think it might do the job.

    Tarek> tny help would be so appreciated.

    Tarek> Cheers,

    Tarek> Tarek Sadek

    Tarek> Ph.D. Student
    Tarek> Department of Civil and Environmental Engineering
    Tarek> University of Melbourne
    Tarek> Parkville, Vic. 3052
    Tarek> Australia

    Tarek> Tel: +61 3 9344 7178
    Tarek> Fax: +61 3 9344 6215

    Tarek> E-mail address: sadek at aqua.civag.unimelb.edu.au
    Tarek> URL: http://www.civag.unimelb.edu.au/~sadek
    Tarek> ------

Tarek,

I think this would be relatively easy to code yourself. The
documentation for the Arc/Info GRID function SLOPE() presents a method
called the ``average maximum technique'', referring to

     Burrough, P.A., 1986.  Principles of Geographical Informatioin
     Systems for Land Resources Assessment.  Oxford University Press,
     New York, p. 50.

The method uses a 3x3 neighborhood in the DEM:

                                a b c
                                d e f
                                g h i
                                   
Slope is computed as 

          _______________
         |  dz 2    dz 2
    S =  | (--)  + (--)
        \|  dx      dy

where 

    dz   (a + 2 d + g) - (c + 2 f + i)
    -- = -----------------------------
    dx              8 cell
                          x

and


    dz   (a + 2 b + c) - (g + 2 h + i)
    -- = -----------------------------
    dy              8 cell
                          x

They also indicate that if any of the neighboring cells are missing
the value of the current cell (e) is assumed for the missing cell.
Given that, I have assumed (and I beleive I'm correct) that aspect
(as an azimuth) would be computed as 

                 dy         dz / dx
   aspect = atan(--) = atan(-------)
                 dx         dz / dy

with proper adjustments for negative angles.  

I have coded something similar in Arc/Info GRID, which uses only 4
neighbors, as in

                                  b 
                                d e f
                                  h 

and handles missing values differently.  I think either method could
be coded in r.mapcalc without difficulty, but I have not done it.

Hope this helps.

Bill

-- 
Bill Perkins
Battelle Pacific Northwest Laboratories
Environmental Technologies Division, Hydrology Group
P.O. Box 999 MSIN K9-36
Richland, Washington, USA  99352
voice: (509) 372-6131 fax: (509) 372-6328 email: wa_perkins at pnl.gov