brewer.lsc

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

# Script to create a slope-aspect map using the recommendations of Brewer 
# and Marlow (1993) http://www.personal.psu.edu/cab38/Terrain/AutoCarto.html
# --------------------------------------------------------------------------

baseDir = "c:\research\reliefPaper\data\";
object = open(baseDir&"MtRainier.srf");
#object = open(baseDir&"Loess.srf");
#object = open(baseDir&"Snowflake.srf");
#object = open(baseDir&"Hemisphere500.srf");

# Calculate slope and aspect
slope  = surfparam(object,"slope",3,1);
aspect = surfparam(object,"aspect",3,1);

# 8 aspect classes - 10 is N, 20 is NE, 30 is E etc.
aspect = ifelse(aspect>337,10,
          ifelse(aspect>292, 80,
           ifelse(aspect>247, 70,
            ifelse(aspect>202, 60,
             ifelse(aspect>157, 50,
              ifelse(aspect>112, 40,
               ifelse(aspect>67, 30,
                ifelse(aspect>22, 20,
                 10))))))));

# 4 slope categories based on <5%, 5-20%, 20-40% and >40%
slope = ifelse(slope>21.8,3,
         ifelse(slope>11.3,2,
           ifelse(slope>2.9,1,
             0)));

# Create new surface based on slope and aspect categories.
brewer = new(object);
brewer = ifelse(slope==0,0,slope+aspect);

# Set up colourtable.
colouredit(brewer,"rules"," 0 161 161 161 (D),"&

                          "11 152 181 129 (D),"&
                          "21 114 168 144 (D),"&
                          "31 124 142 173 (D),"&
                          "41 140 117 160 (D),"&
                          "51 180 123 161 (D),"&
                          "61 203 139 143 (D),"&
                          "71 197 165 138 (D),"&
                          "81 189 191 137 (D),"&

                          "12 141 196  88 (D),"&
                          "22  61 171 113 (D),"&
                          "32  80 120 182 (D),"&
                          "42 119  71 157 (D),"&
                          "52 192  77 156 (D),"&
                          "62 231 111 122 (D),"&
                          "72 226 166 108 (D),"&
                          "82 214 219  94 (D),"&

                          "13 132 214   0 (D),"&
                          "23   0 171  68 (D),"&
                          "33   0 104 192 (D),"&
                          "43 108   0 163 (D),"&
                          "53 202   0 156 (D),"&
                          "63 255  85 104 (D),"&
                          "73 255 171  71 (D),"&
                          "83 244 250   0 (D)");

#Save the reclassified surface
save(brewer,baseDir&"BrewerRainier.srf");