1 #include <unistd.h>
2 #include <sys/types.h>
3 #include <stdlib.h>
4 #include <grass/gis.h>
5 #include <grass/raster.h>
6 #include <grass/glocale.h>
7
8 struct whereandwhat
9 {
10 off_t offset;
11 CELL *p;
12 };
13
recurse_cell(CELL flag,int i,int j,int nl,int ns,struct whereandwhat bas[],struct whereandwhat dir[])14 int recurse_cell(CELL flag, int i, int j, int nl, int ns,
15 struct whereandwhat bas[], struct whereandwhat dir[])
16 {
17 CELL edge;
18 int rc = 0;
19
20 if (j == 0 && j >= ns - 1)
21 return rc;
22
23 if (bas[i].p[j] != flag) {
24 rc = 1;
25 bas[i].p[j] = flag;
26 }
27
28 if (i > 0) {
29 edge = dir[i - 1].p[j - 1];
30 if (bas[i - 1].p[j - 1] == -1 && !Rast_is_c_null_value(&edge) &&
31 edge == 4)
32 rc += recurse_cell(flag, i - 1, j - 1, nl, ns, bas, dir);
33 edge = dir[i - 1].p[j];
34 if (bas[i - 1].p[j] == -1 && !Rast_is_c_null_value(&edge) && edge == 8)
35 rc += recurse_cell(flag, i - 1, j, nl, ns, bas, dir);
36 edge = dir[i - 1].p[j + 1];
37 if (bas[i - 1].p[j + 1] == -1 && !Rast_is_c_null_value(&edge) &&
38 edge == 16)
39 rc += recurse_cell(flag, i - 1, j + 1, nl, ns, bas, dir);
40
41 }
42
43 edge = dir[i].p[j - 1];
44 if (bas[i].p[j - 1] == -1 && !Rast_is_c_null_value(&edge) && edge == 2)
45 rc += recurse_cell(flag, i, j - 1, nl, ns, bas, dir);
46
47 edge = dir[i].p[j + 1];
48 if (bas[i].p[j + 1] == -1 && !Rast_is_c_null_value(&edge) && edge == 32)
49 rc += recurse_cell(flag, i, j + 1, nl, ns, bas, dir);
50
51 if (i < nl - 1) {
52 edge = dir[i + 1].p[j - 1];
53 if (bas[i + 1].p[j - 1] == -1 && !Rast_is_c_null_value(&edge) &&
54 edge == 1)
55 rc += recurse_cell(flag, i + 1, j - 1, nl, ns, bas, dir);
56 edge = dir[i + 1].p[j];
57 if (bas[i + 1].p[j] == -1 && !Rast_is_c_null_value(&edge) && edge == 128)
58 rc += recurse_cell(flag, i + 1, j, nl, ns, bas, dir);
59 edge = dir[i + 1].p[j + 1];
60 if (bas[i + 1].p[j + 1] == -1 && !Rast_is_c_null_value(&edge) &&
61 edge == 64)
62 rc += recurse_cell(flag, i + 1, j + 1, nl, ns, bas, dir);
63 }
64 return rc;
65 }
66
wtrshed(int fm,int fd,int nl,int ns,int mxbuf)67 void wtrshed(int fm, int fd, int nl, int ns, int mxbuf)
68 {
69 int pass, repeat, flag, i, j, half, bufsz;
70 int sline, nline, rdline;
71
72 struct whereandwhat hold;
73 struct whereandwhat *dir;
74 struct whereandwhat *bas;
75
76 dir = G_malloc(mxbuf * sizeof(struct whereandwhat));
77 bas = G_malloc(mxbuf * sizeof(struct whereandwhat));
78
79 bufsz = ns * sizeof(CELL);
80
81 /* adjust maxbuf to an even number */
82 half = mxbuf / 2;
83 mxbuf = 2 * half;
84
85 /* allocate buffers for drainage directions and basin areas */
86 for (i = 0; i < mxbuf; i += 1)
87 bas[i].p = (CELL *) G_calloc(ns, sizeof(CELL));
88 for (i = 0; i < mxbuf; i += 1)
89 dir[i].p = (CELL *) G_calloc(ns, sizeof(CELL));
90
91 pass = 0;
92
93 /* complete a downward pass */
94 do {
95 G_verbose_message(_("Watershed pass %d"), ++pass);
96 repeat = 0;
97
98 /* fill the buffer */
99 nline = mxbuf;
100 sline = 0;
101 rdline = 1;
102 for (i = 0; i < mxbuf; i++) {
103 bas[i].offset = dir[i].offset = (off_t) rdline *bufsz;
104
105 lseek(fm, bas[i].offset, SEEK_SET);
106 read(fm, bas[i].p, bufsz);
107
108 lseek(fd, dir[i].offset, SEEK_SET);
109 read(fd, dir[i].p, bufsz);
110
111 rdline++;
112 }
113
114 /* repeat for all subsequent rows except the first and last */
115 for (i = 1; i < nl - 1; i += 1) {
116 /* analyse one line */
117 for (j = 1; j < ns - 1; j += 1) {
118 flag = bas[sline].p[j];
119 if (flag > 0)
120 if (recurse_cell(flag, sline, j, nline, ns, bas, dir) > 0)
121 repeat = 1;
122 }
123
124 /* write one line */
125 lseek(fm, bas[sline].offset, SEEK_SET);
126 write(fm, bas[sline].p, bufsz);
127
128 /* If the bottom end of the buffers reach the bottom of the file,
129 * rotate the buffers and read new lines */
130 if (rdline < nl - 1) {
131 hold = bas[0];
132 for (j = 1; j < mxbuf; j += 1)
133 bas[j - 1] = bas[j];
134 bas[mxbuf - 1] = hold;
135
136 hold = dir[0];
137 for (j = 1; j < mxbuf; j += 1)
138 dir[j - 1] = dir[j];
139 dir[mxbuf - 1] = hold;
140
141 bas[mxbuf - 1].offset = dir[mxbuf - 1].offset =
142 (off_t) rdline *bufsz;
143
144 lseek(fm, bas[mxbuf - 1].offset, SEEK_SET);
145 read(fm, bas[mxbuf - 1].p, bufsz);
146
147 lseek(fd, dir[mxbuf - 1].offset, SEEK_SET);
148 read(fd, dir[mxbuf - 1].p, bufsz);
149
150 rdline++;
151 }
152 /* After the buffer reaches the bottom of the file, stop reading file,
153 * just advance the pointers */
154 else {
155 nline -= 1;
156 sline += 1;
157 }
158
159 }
160
161 /* fill the buffer */
162 nline = mxbuf;
163 rdline = nl - 2;
164 for (i = mxbuf - 1; i >= 0; i -= 1) {
165 bas[i].offset = dir[i].offset = (off_t) rdline *bufsz;
166
167 lseek(fm, bas[i].offset, SEEK_SET);
168 read(fm, bas[i].p, bufsz);
169
170 lseek(fd, dir[i].offset, SEEK_SET);
171 read(fd, dir[i].p, bufsz);
172
173 rdline--;
174 }
175
176 /* repeat */
177 for (i = nl - 2; i >= 1; i -= 1) {
178 /* analyse one line */
179 for (j = 1; j < ns - 1; j += 1) {
180 flag = bas[nline - 1].p[j];
181 if (flag > 0)
182 if (recurse_cell(flag, nline - 1, j, nline, ns, bas, dir)
183 > 0)
184 repeat = 1;
185 }
186
187 /* write one line */
188 lseek(fm, bas[nline - 1].offset, SEEK_SET);
189 write(fm, bas[nline - 1].p, bufsz);
190
191 /* Until the top of the buffers reach the top of the file,
192 * rotate the buffers and read new lines */
193 if (rdline >= 1) {
194 hold = bas[nline - 1];
195 for (j = nline - 1; j > 0; j -= 1)
196 bas[j] = bas[j - 1];
197 bas[0] = hold;
198
199 hold = dir[nline - 1];
200 for (j = nline - 1; j > 0; j -= 1)
201 dir[j] = dir[j - 1];
202 dir[0] = hold;
203
204 bas[0].offset = dir[0].offset = (off_t) rdline *bufsz;
205
206 lseek(fm, bas[0].offset, SEEK_SET);
207 read(fm, bas[0].p, bufsz);
208
209 lseek(fd, dir[0].offset, SEEK_SET);
210 read(fd, dir[0].p, bufsz);
211
212 rdline--;
213 }
214 /* after the buffer reaches the top of the file, just decrease the
215 * line count */
216 else
217 nline -= 1;
218
219 }
220
221 } while (repeat);
222
223 /* allocate buffers for drainage directions and basin areas */
224 for (i = 0; i < mxbuf; i += 1)
225 G_free(bas[i].p);
226 for (i = 0; i < mxbuf; i += 1)
227 G_free(dir[i].p);
228
229 G_free(dir);
230 G_free(bas);
231 }
232