source: grass/branches/releasebranch_6_4/raster/r.slope.aspect/main.c

Last change on this file was 62029, checked in by neteler, 10 years ago

user message harmonization update (partial backport of r62007, trac #2409)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 41.9 KB
Line 
1
2/****************************************************************************
3 *
4 * MODULE: r.slope.aspect
5 * AUTHOR(S): Michael Shapiro and
6 * Olga Waupotitsch (original CERL contributors),
7 * Markus Neteler <neteler itc.it>,
8 * Bernhard Reiter <bernhard intevation.de>,
9 * Brad Douglas <rez touchofmadness.com>,
10 * Glynn Clements <glynn gclements.plus.com>,
11 * Hamish Bowman <hamish_b yahoo.com>,
12 * Jachym Cepicky <jachym les-ejk.cz>,
13 * Jan-Oliver Wagner <jan intevation.de>,
14 * Radim Blazek <radim.blazek gmail.com>
15 *
16 * PURPOSE: Generates raster maps of slope, aspect, curvatures and
17 * first and second order partial derivatives from a raster map
18 * of true elevation values
19 *
20 * COPYRIGHT: (C) 1999-2011 by the GRASS Development Team
21 *
22 * This program is free software under the GNU General Public
23 * License (>=v2). Read the file COPYING that comes with GRASS
24 * for details.
25 *
26 *****************************************************************************/
27#include <stdlib.h>
28#include <string.h>
29#include <math.h>
30#include <grass/gis.h>
31#include <grass/glocale.h>
32#include "local_proto.h"
33
34/* 10/99 from GMSL, updated to new GRASS 5 code style , changed default "prec" to float */
35
36
37#define abs(x) ((x)<0?-(x):(x))
38
39
40/**************************************************************************
41 * input is from command line.
42 * arguments are elevation file, slope file, aspect file, profile curvature
43 * file and tangential curvature file
44 * elevation filename required
45 * either slope filename or aspect filename or profile curvature filename
46 * or tangential curvature filename required
47 * usage: r.slope.aspect [-av] elevation=input slope=output1 aspect=output2
48 * pcurv=output3 tcurv=output4 format=name prec=name zfactor=value
49 * min_slp_allowed=value dx=output5 dy=output6 dxx=output7
50 * dyy=output8 dxy=output9
51 * -a don't align window
52 * -q quiet
53 **************************************************************************/
54
55/* some changes made to code to retrieve correct distances when using
56 lat/lon projection. changes involve recalculating H and V. see
57 comments within code. */
58/* added colortables for topographic parameters and fixed
59 * the sign bug for the second order derivatives (jh) */
60
61
62int main(int argc, char *argv[])
63{
64 struct Categories cats;
65 int elevation_fd;
66 int aspect_fd;
67 int slope_fd;
68 int pcurv_fd;
69 int tcurv_fd;
70 int dx_fd;
71 int dy_fd;
72 int dxx_fd;
73 int dyy_fd;
74 int dxy_fd;
75 DCELL *elev_cell[3], *temp;
76 DCELL *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9;
77 DCELL tmp1, tmp2;
78 FCELL dat1, dat2;
79 void *asp_raster, *asp_ptr = NULL;
80 void *slp_raster, *slp_ptr = NULL;
81 void *pcurv_raster, *pcurv_ptr = NULL;
82 void *tcurv_raster, *tcurv_ptr = NULL;
83 void *dx_raster, *dx_ptr = NULL;
84 void *dy_raster, *dy_ptr = NULL;
85 void *dxx_raster, *dxx_ptr = NULL;
86 void *dyy_raster, *dyy_ptr = NULL;
87 void *dxy_raster, *dxy_ptr = NULL;
88 int i;
89 RASTER_MAP_TYPE out_type = 0, data_type;
90 int Wrap; /* global wraparound */
91 struct Cell_head window, cellhd;
92 struct History hist;
93 struct Colors colors;
94
95 char *elev_name;
96 char *aspect_name;
97 char *slope_name;
98 char *pcurv_name;
99 char *tcurv_name;
100 char *dx_name;
101 char *dy_name;
102 char *dxx_name;
103 char *dyy_name;
104 char *dxy_name;
105 char buf[300];
106 char *mapset;
107 int nrows, row;
108 int ncols, col;
109
110 double north, east, south, west, ns_med;
111
112 double radians_to_degrees;
113 double degrees_to_radians;
114 double H, V;
115 double dx; /* partial derivative in ew direction */
116 double dy; /* partial derivative in ns direction */
117 double dxx, dxy, dyy;
118 double s3, s4, s5, s6;
119 double pcurv, tcurv;
120 double scik1 = 100000.;
121 double zfactor;
122 double factor;
123 double aspect, min_asp = 360., max_asp = 0.;
124 double dnorm1, dx2, dy2, grad2, grad, dxy2;
125 double gradmin = 0.001;
126 double c1min = 0., c1max = 0., c2min = 0., c2max = 0.;
127
128 double answer[92];
129 double degrees;
130 double tan_ans;
131 double key;
132 double slp_in_perc, slp_in_deg;
133 double min_slp = 900., max_slp = 0., min_slp_allowed;
134 int low, hi, test = 0;
135 int deg = 0;
136 int perc = 0;
137 char *slope_fmt;
138 struct GModule *module;
139 struct
140 {
141 struct Option *elevation, *slope_fmt, *slope, *aspect, *pcurv, *tcurv,
142 *zfactor, *min_slp_allowed, *out_precision,
143 *dx, *dy, *dxx, *dyy, *dxy;
144 } parm;
145 struct
146 {
147 struct Flag *a;
148 /* please, remove before GRASS 7 released */
149 struct Flag *q;
150 } flag;
151
152 G_gisinit(argv[0]);
153
154 module = G_define_module();
155 module->keywords = _("raster, terrain");
156 module->label = _("Generates raster maps of slope, aspect, curvatures and "
157 "partial derivatives from a elevation raster map.");
158 module->description = _("Aspect is calculated counterclockwise from east.");
159
160 parm.elevation = G_define_standard_option(G_OPT_R_ELEV);
161
162 parm.slope = G_define_standard_option(G_OPT_R_OUTPUT);
163 parm.slope->key = "slope";
164 parm.slope->required = NO;
165 parm.slope->description = _("Name for output slope raster map");
166 parm.slope->guisection = _("Outputs");
167
168 parm.aspect = G_define_standard_option(G_OPT_R_OUTPUT);
169 parm.aspect->key = "aspect";
170 parm.aspect->required = NO;
171 parm.aspect->description = _("Name for output aspect raster map");
172 parm.aspect->guisection = _("Outputs");
173
174 parm.slope_fmt = G_define_option();
175 parm.slope_fmt->key = "format";
176 parm.slope_fmt->type = TYPE_STRING;
177 parm.slope_fmt->required = NO;
178 parm.slope_fmt->answer = "degrees";
179 parm.slope_fmt->options = "degrees,percent";
180 parm.slope_fmt->description = _("Format for reporting the slope");
181 parm.slope_fmt->guisection = _("Settings");
182
183 parm.out_precision = G_define_option();
184 parm.out_precision->key = "prec";
185 parm.out_precision->type = TYPE_STRING;
186 parm.out_precision->required = NO;
187 parm.out_precision->answer = "float";
188 parm.out_precision->options = "default,double,float,int";
189 parm.out_precision->description =
190 _("Type of output aspect and slope maps");
191 parm.out_precision->guisection = _("Settings");
192
193 parm.pcurv = G_define_standard_option(G_OPT_R_OUTPUT);
194 parm.pcurv->key = "pcurv";
195 parm.pcurv->required = NO;
196 parm.pcurv->description =
197 _("Name for output profile curvature raster map");
198 parm.pcurv->guisection = _("Outputs");
199
200 parm.tcurv = G_define_standard_option(G_OPT_R_OUTPUT);
201 parm.tcurv->key = "tcurv";
202 parm.tcurv->required = NO;
203 parm.tcurv->description =
204 _("Name for output tangential curvature raster map");
205 parm.tcurv->guisection = _("Outputs");
206
207 parm.dx = G_define_standard_option(G_OPT_R_OUTPUT);
208 parm.dx->key = "dx";
209 parm.dx->required = NO;
210 parm.dx->description =
211 _("Name for output first order partial derivative dx (E-W slope) raster map");
212 parm.dx->guisection = _("Outputs");
213
214 parm.dy = G_define_standard_option(G_OPT_R_OUTPUT);
215 parm.dy->key = "dy";
216 parm.dy->required = NO;
217 parm.dy->description =
218 _("Name for output first order partial derivative dy (N-S slope) raster map");
219 parm.dy->guisection = _("Outputs");
220
221 parm.dxx = G_define_standard_option(G_OPT_R_OUTPUT);
222 parm.dxx->key = "dxx";
223 parm.dxx->required = NO;
224 parm.dxx->description =
225 _("Name for output second order partial derivative dxx raster map");
226 parm.dxx->guisection = _("Outputs");
227
228 parm.dyy = G_define_standard_option(G_OPT_R_OUTPUT);
229 parm.dyy->key = "dyy";
230 parm.dyy->required = NO;
231 parm.dyy->description =
232 _("Name for output second order partial derivative dyy raster map");
233 parm.dyy->guisection = _("Outputs");
234
235 parm.dxy = G_define_standard_option(G_OPT_R_OUTPUT);
236 parm.dxy->key = "dxy";
237 parm.dxy->required = NO;
238 parm.dxy->description =
239 _("Name for output second order partial derivative dxy raster map");
240 parm.dxy->guisection = _("Outputs");
241
242 parm.zfactor = G_define_option();
243 parm.zfactor->key = "zfactor";
244 parm.zfactor->description =
245 _("Multiplicative factor to convert elevation units to meters");
246 parm.zfactor->type = TYPE_DOUBLE;
247 parm.zfactor->required = NO;
248 parm.zfactor->answer = "1.0";
249 parm.zfactor->guisection = _("Settings");
250
251 parm.min_slp_allowed = G_define_option();
252 parm.min_slp_allowed->key = "min_slp_allowed";
253 parm.min_slp_allowed->description =
254 _("Minimum slope value (in percent) for which aspect is computed");
255 parm.min_slp_allowed->type = TYPE_DOUBLE;
256 parm.min_slp_allowed->required = NO;
257 parm.min_slp_allowed->answer = "0.0";
258 parm.min_slp_allowed->guisection = _("Settings");
259
260 /* please, remove before GRASS 7 released */
261 flag.q = G_define_flag();
262 flag.q->key = 'q';
263 flag.q->description = _("Quiet");
264
265 flag.a = G_define_flag();
266 flag.a->key = 'a';
267 flag.a->description =
268 _("Do not align the current region to the elevation layer");
269 flag.a->guisection = _("Settings");
270
271
272 radians_to_degrees = 180.0 / M_PI;
273 degrees_to_radians = M_PI / 180.0;
274
275 /* INC BY ONE
276 answer[0] = 0.0;
277 answer[91] = 15000.0;
278
279 for (i = 1; i < 91; i++)
280 {
281 degrees = i - .5;
282 tan_ans = tan ( degrees / radians_to_degrees );
283 answer[i] = tan_ans * tan_ans;
284 }
285 */
286 answer[0] = 0.0;
287 answer[90] = 15000.0;
288
289 for (i = 0; i < 90; i++) {
290 degrees = i + .5;
291 tan_ans = tan(degrees / radians_to_degrees);
292 answer[i] = tan_ans * tan_ans;
293 }
294
295 if (G_parser(argc, argv))
296 exit(EXIT_FAILURE);
297
298 /* please, remove before GRASS 7 released */
299 if (flag.q->answer) {
300 putenv("GRASS_VERBOSE=0");
301 G_warning(_("The '-q' flag is superseded and will be removed "
302 "in future. Please use '--quiet' instead."));
303 }
304
305
306
307 elev_name = parm.elevation->answer;
308 slope_name = parm.slope->answer;
309 aspect_name = parm.aspect->answer;
310 pcurv_name = parm.pcurv->answer;
311 tcurv_name = parm.tcurv->answer;
312 dx_name = parm.dx->answer;
313 dy_name = parm.dy->answer;
314 dxx_name = parm.dxx->answer;
315 dyy_name = parm.dyy->answer;
316 dxy_name = parm.dxy->answer;
317
318 G_check_input_output_name(elev_name, slope_name, GR_FATAL_EXIT);
319 G_check_input_output_name(elev_name, aspect_name, GR_FATAL_EXIT);
320 G_check_input_output_name(elev_name, pcurv_name, GR_FATAL_EXIT);
321 G_check_input_output_name(elev_name, tcurv_name, GR_FATAL_EXIT);
322 G_check_input_output_name(elev_name, dx_name, GR_FATAL_EXIT);
323 G_check_input_output_name(elev_name, dy_name, GR_FATAL_EXIT);
324 G_check_input_output_name(elev_name, dxx_name, GR_FATAL_EXIT);
325 G_check_input_output_name(elev_name, dyy_name, GR_FATAL_EXIT);
326 G_check_input_output_name(elev_name, dxy_name, GR_FATAL_EXIT);
327
328 if (sscanf(parm.zfactor->answer, "%lf", &zfactor) != 1 || zfactor <= 0.0) {
329 G_fatal_error(_("%s=%s - must be a positive number"),
330 parm.zfactor->key, parm.zfactor->answer);
331 }
332
333 if (sscanf(parm.min_slp_allowed->answer, "%lf", &min_slp_allowed) != 1 ||
334 min_slp_allowed < 0.0) {
335 G_fatal_error(_("%s=%s - must be a non-negative number"),
336 parm.min_slp_allowed->key,
337 parm.min_slp_allowed->answer);
338 }
339
340 slope_fmt = parm.slope_fmt->answer;
341 if (strcmp(slope_fmt, "percent") == 0)
342 perc = 1;
343 else if (strcmp(slope_fmt, "degrees") == 0)
344 deg = 1;
345
346 if (slope_name == NULL && aspect_name == NULL
347 && pcurv_name == NULL && tcurv_name == NULL
348 && dx_name == NULL && dy_name == NULL
349 && dxx_name == NULL && dyy_name == NULL && dxy_name == NULL) {
350 G_fatal_error(_("You must specify at least one of the parameters: "
351 "<%s>, <%s>, <%s>, <%s>, <%s>, <%s>, <%s>, <%s> or <%s>"),
352 parm.slope->key, parm.aspect->key, parm.pcurv->key,
353 parm.tcurv->key, parm.dx->key, parm.dy->key,
354 parm.dxx->key, parm.dyy->key, parm.dxy->key);
355 }
356
357 /* check elevation file existence */
358 mapset = G_find_cell2(elev_name, "");
359 if (!mapset)
360 G_fatal_error(_("Raster map <%s> not found"), elev_name);
361
362 /* set the window from the header for the elevation file */
363 if (!flag.a->answer) {
364 G_get_window(&window);
365 if (G_get_cellhd(elev_name, mapset, &cellhd) >= 0) {
366 G_align_window(&window, &cellhd);
367 G_set_window(&window);
368 }
369 }
370
371 if (strcmp(parm.out_precision->answer, "double") == 0)
372 out_type = DCELL_TYPE;
373 else if (strcmp(parm.out_precision->answer, "float") == 0)
374 out_type = FCELL_TYPE;
375 else if (strcmp(parm.out_precision->answer, "int") == 0)
376 out_type = CELL_TYPE;
377 else if (strcmp(parm.out_precision->answer, "default") == 0)
378 out_type = -1;
379 else
380 G_fatal_error(_("Wrong raster type: %s"), parm.out_precision->answer);
381
382 data_type = out_type;
383 if (data_type < 0)
384 data_type = DCELL_TYPE;
385 /* data type is the type of data being processed,
386 out_type is type of map being created */
387
388 G_get_set_window(&window);
389
390 nrows = G_window_rows();
391 ncols = G_window_cols();
392
393 if (((window.west == (window.east - 360.))
394 || (window.east == (window.west - 360.))) &&
395 (G_projection() == PROJECTION_LL)) {
396 Wrap = 1;
397 ncols += 2;
398 }
399 else
400 Wrap = 0;
401
402 /* H = window.ew_res * 4 * 2/ zfactor; *//* horizontal (east-west) run
403 times 4 for weighted difference */
404 /* V = window.ns_res * 4 * 2/ zfactor; *//* vertical (north-south) run
405 times 4 for weighted difference */
406
407 /* give warning if location units are different from meters and zfactor=1 */
408 factor = G_database_units_to_meters_factor();
409 if (factor != 1.0)
410 G_warning(_("Converting units to meters, factor=%.6f"), factor);
411
412 G_begin_distance_calculations();
413 north = G_row_to_northing(0.5, &window);
414 ns_med = G_row_to_northing(1.5, &window);
415 south = G_row_to_northing(2.5, &window);
416 east = G_col_to_easting(2.5, &window);
417 west = G_col_to_easting(0.5, &window);
418 V = G_distance(east, north, east, south) * 4 / zfactor;
419 H = G_distance(east, ns_med, west, ns_med) * 4 / zfactor;
420 /* ____________________________
421 |c1 |c2 |c3 |
422 | | | |
423 | | north | |
424 | | | |
425 |________|________|________|
426 |c4 |c5 |c6 |
427 | | | |
428 | east | ns_med | west |
429 | | | |
430 |________|________|________|
431 |c7 |c8 |c9 |
432 | | | |
433 | | south | |
434 | | | |
435 |________|________|________|
436 */
437
438 /* open the elevation file for reading */
439 elevation_fd = G_open_cell_old(elev_name, mapset);
440 if (elevation_fd < 0)
441 exit(EXIT_FAILURE);
442 elev_cell[0] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
443 G_set_d_null_value(elev_cell[0], ncols);
444 elev_cell[1] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
445 G_set_d_null_value(elev_cell[1], ncols);
446 elev_cell[2] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
447 G_set_d_null_value(elev_cell[2], ncols);
448
449 if (slope_name != NULL) {
450 slope_fd = opennew(slope_name, out_type);
451 slp_raster = G_allocate_raster_buf(data_type);
452 G_set_null_value(slp_raster, G_window_cols(), data_type);
453 G_put_raster_row(slope_fd, slp_raster, data_type);
454 }
455 else {
456 slp_raster = NULL;
457 slope_fd = -1;
458 }
459
460 if (aspect_name != NULL) {
461 aspect_fd = opennew(aspect_name, out_type);
462 asp_raster = G_allocate_raster_buf(data_type);
463 G_set_null_value(asp_raster, G_window_cols(), data_type);
464 G_put_raster_row(aspect_fd, asp_raster, data_type);
465 }
466 else {
467 asp_raster = NULL;
468 aspect_fd = -1;
469 }
470
471 if (pcurv_name != NULL) {
472 pcurv_fd = opennew(pcurv_name, out_type);
473 pcurv_raster = G_allocate_raster_buf(data_type);
474 G_set_null_value(pcurv_raster, G_window_cols(), data_type);
475 G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
476 }
477 else {
478 pcurv_raster = NULL;
479 pcurv_fd = -1;
480 }
481
482 if (tcurv_name != NULL) {
483 tcurv_fd = opennew(tcurv_name, out_type);
484 tcurv_raster = G_allocate_raster_buf(data_type);
485 G_set_null_value(tcurv_raster, G_window_cols(), data_type);
486 G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
487 }
488 else {
489 tcurv_raster = NULL;
490 tcurv_fd = -1;
491 }
492
493 if (dx_name != NULL) {
494 dx_fd = opennew(dx_name, out_type);
495 dx_raster = G_allocate_raster_buf(data_type);
496 G_set_null_value(dx_raster, G_window_cols(), data_type);
497 G_put_raster_row(dx_fd, dx_raster, data_type);
498 }
499 else {
500 dx_raster = NULL;
501 dx_fd = -1;
502 }
503
504 if (dy_name != NULL) {
505 dy_fd = opennew(dy_name, out_type);
506 dy_raster = G_allocate_raster_buf(data_type);
507 G_set_null_value(dy_raster, G_window_cols(), data_type);
508 G_put_raster_row(dy_fd, dy_raster, data_type);
509 }
510 else {
511 dy_raster = NULL;
512 dy_fd = -1;
513 }
514
515 if (dxx_name != NULL) {
516 dxx_fd = opennew(dxx_name, out_type);
517 dxx_raster = G_allocate_raster_buf(data_type);
518 G_set_null_value(dxx_raster, G_window_cols(), data_type);
519 G_put_raster_row(dxx_fd, dxx_raster, data_type);
520 }
521 else {
522 dxx_raster = NULL;
523 dxx_fd = -1;
524 }
525
526 if (dyy_name != NULL) {
527 dyy_fd = opennew(dyy_name, out_type);
528 dyy_raster = G_allocate_raster_buf(data_type);
529 G_set_null_value(dyy_raster, G_window_cols(), data_type);
530 G_put_raster_row(dyy_fd, dyy_raster, data_type);
531 }
532 else {
533 dyy_raster = NULL;
534 dyy_fd = -1;
535 }
536
537 if (dxy_name != NULL) {
538 dxy_fd = opennew(dxy_name, out_type);
539 dxy_raster = G_allocate_raster_buf(data_type);
540 G_set_null_value(dxy_raster, G_window_cols(), data_type);
541 G_put_raster_row(dxy_fd, dxy_raster, data_type);
542 }
543 else {
544 dxy_raster = NULL;
545 dxy_fd = -1;
546 }
547
548 if (aspect_fd < 0 && slope_fd < 0 && pcurv_fd < 0 && tcurv_fd < 0
549 && dx_fd < 0 && dy_fd < 0 && dxx_fd < 0 && dyy_fd < 0 && dxy_fd < 0)
550 exit(EXIT_FAILURE);
551
552 if (Wrap) {
553 G_get_d_raster_row_nomask(elevation_fd, elev_cell[1] + 1, 0);
554 elev_cell[1][0] = elev_cell[1][G_window_cols() - 1];
555 elev_cell[1][G_window_cols() + 1] = elev_cell[1][2];
556 }
557 else
558 G_get_d_raster_row_nomask(elevation_fd, elev_cell[1], 0);
559
560 if (Wrap) {
561 G_get_d_raster_row_nomask(elevation_fd, elev_cell[2] + 1, 1);
562 elev_cell[2][0] = elev_cell[2][G_window_cols() - 1];
563 elev_cell[2][G_window_cols() + 1] = elev_cell[2][2];
564 }
565 else
566 G_get_d_raster_row_nomask(elevation_fd, elev_cell[2], 1);
567
568 G_verbose_message(_("Percent complete..."));
569
570 for (row = 2; row < nrows; row++) {
571 /* if projection is Lat/Lon, recalculate V and H */
572 if (G_projection() == PROJECTION_LL) {
573 north = G_row_to_northing((row - 2 + 0.5), &window);
574 ns_med = G_row_to_northing((row - 1 + 0.5), &window);
575 south = G_row_to_northing((row + 0.5), &window);
576 east = G_col_to_easting(2.5, &window);
577 west = G_col_to_easting(0.5, &window);
578 V = G_distance(east, north, east, south) * 4 / zfactor;
579 H = G_distance(east, ns_med, west, ns_med) * 4 / zfactor;
580 /* ____________________________
581 |c1 |c2 |c3 |
582 | | | |
583 | | north | |
584 | | | |
585 |________|________|________|
586 |c4 |c5 |c6 |
587 | | | |
588 | east | ns_med | west |
589 | | | |
590 |________|________|________|
591 |c7 |c8 |c9 |
592 | | | |
593 | | south | |
594 | | | |
595 |________|________|________|
596 */
597 }
598
599 G_percent(row, nrows, 2);
600 temp = elev_cell[0];
601 elev_cell[0] = elev_cell[1];
602 elev_cell[1] = elev_cell[2];
603 elev_cell[2] = temp;
604
605 if (Wrap) {
606 G_get_d_raster_row_nomask(elevation_fd, elev_cell[2] + 1, row);
607 elev_cell[2][0] = elev_cell[2][G_window_cols() - 1];
608 elev_cell[2][G_window_cols() + 1] = elev_cell[2][2];
609 }
610 else
611 G_get_d_raster_row_nomask(elevation_fd, elev_cell[2], row);
612
613 c1 = elev_cell[0];
614 c2 = c1 + 1;
615 c3 = c1 + 2;
616 c4 = elev_cell[1];
617 c5 = c4 + 1;
618 c6 = c4 + 2;
619 c7 = elev_cell[2];
620 c8 = c7 + 1;
621 c9 = c7 + 2;
622
623 if (aspect_fd >= 0) {
624 if (Wrap)
625 asp_ptr = asp_raster;
626 else
627 asp_ptr =
628 G_incr_void_ptr(asp_raster, G_raster_size(data_type));
629 }
630 if (slope_fd >= 0) {
631 if (Wrap)
632 slp_ptr = slp_raster;
633 else
634 slp_ptr =
635 G_incr_void_ptr(slp_raster, G_raster_size(data_type));
636 }
637
638 if (pcurv_fd >= 0) {
639 if (Wrap)
640 pcurv_ptr = pcurv_raster;
641 else
642 pcurv_ptr =
643 G_incr_void_ptr(pcurv_raster, G_raster_size(data_type));
644 }
645
646 if (tcurv_fd >= 0) {
647 if (Wrap)
648 tcurv_ptr = tcurv_raster;
649 else
650 tcurv_ptr =
651 G_incr_void_ptr(tcurv_raster, G_raster_size(data_type));
652 }
653
654 if (dx_fd >= 0) {
655 if (Wrap)
656 dx_ptr = dx_raster;
657 else
658 dx_ptr = G_incr_void_ptr(dx_raster, G_raster_size(data_type));
659 }
660
661 if (dy_fd >= 0) {
662 if (Wrap)
663 dy_ptr = dy_raster;
664 else
665 dy_ptr = G_incr_void_ptr(dy_raster, G_raster_size(data_type));
666 }
667
668 if (dxx_fd >= 0) {
669 if (Wrap)
670 dxx_ptr = dxx_raster;
671 else
672 dxx_ptr =
673 G_incr_void_ptr(dxx_raster, G_raster_size(data_type));
674 }
675
676 if (dyy_fd >= 0) {
677 if (Wrap)
678 dyy_ptr = dyy_raster;
679 else
680 dyy_ptr =
681 G_incr_void_ptr(dyy_raster, G_raster_size(data_type));
682 }
683
684 if (dxy_fd >= 0) {
685 if (Wrap)
686 dxy_ptr = dxy_raster;
687 else
688 dxy_ptr =
689 G_incr_void_ptr(dxy_raster, G_raster_size(data_type));
690 }
691
692
693 /*skip first cell of the row */
694
695 for (col = ncols - 2; col-- > 0;
696 c1++, c2++, c3++, c4++, c5++, c6++, c7++, c8++, c9++) {
697 /* DEBUG:
698 fprintf(stdout, "\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n",
699 *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9);
700 */
701
702 if (G_is_d_null_value(c1) || G_is_d_null_value(c2) ||
703 G_is_d_null_value(c3) || G_is_d_null_value(c4) ||
704 G_is_d_null_value(c5) || G_is_d_null_value(c6) ||
705 G_is_d_null_value(c7) || G_is_d_null_value(c8) ||
706 G_is_d_null_value(c9)) {
707 if (slope_fd > 0) {
708 G_set_null_value(slp_ptr, 1, data_type);
709 slp_ptr =
710 G_incr_void_ptr(slp_ptr, G_raster_size(data_type));
711 }
712 if (aspect_fd > 0) {
713 G_set_null_value(asp_ptr, 1, data_type);
714 asp_ptr =
715 G_incr_void_ptr(asp_ptr, G_raster_size(data_type));
716 }
717 if (pcurv_fd > 0) {
718 G_set_null_value(pcurv_ptr, 1, data_type);
719 pcurv_ptr =
720 G_incr_void_ptr(pcurv_ptr, G_raster_size(data_type));
721 }
722 if (tcurv_fd > 0) {
723 G_set_null_value(tcurv_ptr, 1, data_type);
724 tcurv_ptr =
725 G_incr_void_ptr(tcurv_ptr, G_raster_size(data_type));
726 }
727 if (dx_fd > 0) {
728 G_set_null_value(dx_ptr, 1, data_type);
729 dx_ptr =
730 G_incr_void_ptr(dx_ptr, G_raster_size(data_type));
731 }
732 if (dy_fd > 0) {
733 G_set_null_value(dy_ptr, 1, data_type);
734 dy_ptr =
735 G_incr_void_ptr(dy_ptr, G_raster_size(data_type));
736 }
737 if (dxx_fd > 0) {
738 G_set_null_value(dxx_ptr, 1, data_type);
739 dxx_ptr =
740 G_incr_void_ptr(dxx_ptr, G_raster_size(data_type));
741 }
742 if (dyy_fd > 0) {
743 G_set_null_value(dyy_ptr, 1, data_type);
744 dyy_ptr =
745 G_incr_void_ptr(dyy_ptr, G_raster_size(data_type));
746 }
747 if (dxy_fd > 0) {
748 G_set_null_value(dxy_ptr, 1, data_type);
749 dxy_ptr =
750 G_incr_void_ptr(dxy_ptr, G_raster_size(data_type));
751 }
752 continue;
753 } /* no data */
754
755 dx = ((*c1 + *c4 + *c4 + *c7) - (*c3 + *c6 + *c6 + *c9)) / H;
756 dy = ((*c7 + *c8 + *c8 + *c9) - (*c1 + *c2 + *c2 + *c3)) / V;
757
758 /* compute topographic parameters */
759 key = dx * dx + dy * dy;
760 slp_in_perc = 100 * sqrt(key);
761 slp_in_deg = atan(sqrt(key)) * radians_to_degrees;
762
763 /* now update min and max */
764 if (deg) {
765 if (min_slp > slp_in_deg)
766 min_slp = slp_in_deg;
767 if (max_slp < slp_in_deg)
768 max_slp = slp_in_deg;
769 }
770 else {
771 if (min_slp > slp_in_perc)
772 min_slp = slp_in_perc;
773 if (max_slp < slp_in_perc)
774 max_slp = slp_in_perc;
775 }
776 if (slp_in_perc < min_slp_allowed)
777 slp_in_perc = 0.;
778
779 if (deg && out_type == CELL_TYPE) {
780 /* INC BY ONE
781 low = 1;
782 hi = 91;
783 */
784 low = 0;
785 hi = 90;
786 test = 20;
787
788 while (hi >= low) {
789 if (key >= answer[test])
790 low = test + 1;
791 else if (key < answer[test - 1])
792 hi = test - 1;
793 else
794 break;
795 test = (low + hi) / 2;
796 }
797 }
798 else if (perc && out_type == CELL_TYPE)
799 /* INCR_BY_ONE */
800 /* test = slp_in_perc + 1.5; *//* All the slope categories are
801 incremented by 1 */
802 test = slp_in_perc + .5;
803
804 if (slope_fd > 0) {
805 if (data_type == CELL_TYPE)
806 *((CELL *) slp_ptr) = (CELL) test;
807 else {
808 if (deg)
809 G_set_raster_value_d(slp_ptr,
810 (DCELL) slp_in_deg, data_type);
811 else
812 G_set_raster_value_d(slp_ptr,
813 (DCELL) slp_in_perc, data_type);
814 }
815 slp_ptr = G_incr_void_ptr(slp_ptr, G_raster_size(data_type));
816 } /* computing slope */
817
818 if (aspect_fd > 0) {
819 if (key == 0.)
820 aspect = 0.;
821 else if (dx == 0) {
822 if (dy > 0)
823 aspect = 90.;
824 else
825 aspect = 270.;
826 }
827 else {
828 aspect = (atan2(dy, dx) / degrees_to_radians);
829 if ((aspect <= 0.5) && (aspect > 0) &&
830 out_type == CELL_TYPE)
831 aspect = 360.;
832 if (aspect <= 0.)
833 aspect = 360. + aspect;
834 }
835
836 /* if it's not the case that the slope for this cell
837 is below specified minimum */
838 if (!((slope_fd > 0) && (slp_in_perc < min_slp_allowed))) {
839 if (out_type == CELL_TYPE)
840 *((CELL *) asp_ptr) = (CELL) (aspect + .5);
841 else
842 G_set_raster_value_d(asp_ptr,
843 (DCELL) aspect, data_type);
844 }
845 else
846 G_set_null_value(asp_ptr, 1, data_type);
847 asp_ptr = G_incr_void_ptr(asp_ptr, G_raster_size(data_type));
848
849 /* now update min and max */
850 if (min_asp > aspect)
851 min_asp = aspect;
852 if (max_asp < aspect)
853 max_asp = aspect;
854 } /* computing aspect */
855
856 if (dx_fd > 0) {
857 if (out_type == CELL_TYPE)
858 *((CELL *) dx_ptr) = (CELL) (scik1 * dx);
859 else
860 G_set_raster_value_d(dx_ptr, (DCELL) dx, data_type);
861 dx_ptr = G_incr_void_ptr(dx_ptr, G_raster_size(data_type));
862 }
863
864 if (dy_fd > 0) {
865 if (out_type == CELL_TYPE)
866 *((CELL *) dy_ptr) = (CELL) (scik1 * dy);
867 else
868 G_set_raster_value_d(dy_ptr, (DCELL) dy, data_type);
869 dy_ptr = G_incr_void_ptr(dy_ptr, G_raster_size(data_type));
870 }
871
872 if (dxx_fd <= 0 && dxy_fd <= 0 && dyy_fd <= 0 &&
873 pcurv_fd <= 0 && tcurv_fd <= 0)
874 continue;
875
876 /* compute second order derivatives */
877 s4 = *c1 + *c3 + *c7 + *c9 - *c5 * 8.;
878 s5 = *c4 * 4. + *c6 * 4. - *c8 * 2. - *c2 * 2.;
879 s6 = *c8 * 4. + *c2 * 4. - *c4 * 2. - *c6 * 2.;
880 s3 = *c7 - *c9 + *c3 - *c1;
881
882 dxx = -(s4 + s5) / ((3. / 32.) * H * H);
883 dyy = -(s4 + s6) / ((3. / 32.) * V * V);
884 dxy = -s3 / ((1. / 16.) * H * V);
885
886 if (dxx_fd > 0) {
887 if (out_type == CELL_TYPE)
888 *((CELL *) dxx_ptr) = (CELL) (scik1 * dxx);
889 else
890 G_set_raster_value_d(dxx_ptr, (DCELL) dxx, data_type);
891 dxx_ptr = G_incr_void_ptr(dxx_ptr, G_raster_size(data_type));
892 }
893
894 if (dyy_fd > 0) {
895 if (out_type == CELL_TYPE)
896 *((CELL *) dyy_ptr) = (CELL) (scik1 * dyy);
897 else
898 G_set_raster_value_d(dyy_ptr, (DCELL) dyy, data_type);
899 dyy_ptr = G_incr_void_ptr(dyy_ptr, G_raster_size(data_type));
900 }
901
902 if (dxy_fd > 0) {
903 if (out_type == CELL_TYPE)
904 *((CELL *) dxy_ptr) = (CELL) (scik1 * dxy);
905 else
906 G_set_raster_value_d(dxy_ptr, (DCELL) dxy, data_type);
907 dxy_ptr = G_incr_void_ptr(dxy_ptr, G_raster_size(data_type));
908 }
909
910 /* compute curvature */
911 if (pcurv_fd <= 0 && tcurv_fd <= 0)
912 continue;
913
914 grad2 = key; /*dx2 + dy2 */
915 grad = sqrt(grad2);
916 if (grad <= gradmin) {
917 pcurv = 0.;
918 tcurv = 0.;
919 }
920 else {
921 dnorm1 = sqrt(grad2 + 1.);
922 dxy2 = 2. * dxy * dx * dy;
923 dx2 = dx * dx;
924 dy2 = dy * dy;
925 pcurv = (dxx * dx2 + dxy2 + dyy * dy2) /
926 (grad2 * dnorm1 * dnorm1 * dnorm1);
927 tcurv = (dxx * dy2 - dxy2 + dyy * dx2) / (grad2 * dnorm1);
928 if (c1min > pcurv)
929 c1min = pcurv;
930 if (c1max < pcurv)
931 c1max = pcurv;
932 if (c2min > tcurv)
933 c2min = tcurv;
934 if (c2max < tcurv)
935 c2max = tcurv;
936 }
937
938 if (pcurv_fd > 0) {
939 if (out_type == CELL_TYPE)
940 *((CELL *) pcurv_ptr) = (CELL) (scik1 * pcurv);
941 else
942 G_set_raster_value_d(pcurv_ptr, (DCELL) pcurv, data_type);
943 pcurv_ptr =
944 G_incr_void_ptr(pcurv_ptr, G_raster_size(data_type));
945 }
946
947 if (tcurv_fd > 0) {
948 if (out_type == CELL_TYPE)
949 *((CELL *) tcurv_ptr) = (CELL) (scik1 * tcurv);
950 else
951 G_set_raster_value_d(tcurv_ptr, (DCELL) tcurv, data_type);
952 tcurv_ptr =
953 G_incr_void_ptr(tcurv_ptr, G_raster_size(data_type));
954 }
955
956 } /* column for loop */
957
958 if (aspect_fd > 0)
959 G_put_raster_row(aspect_fd, asp_raster, data_type);
960
961 if (slope_fd > 0)
962 G_put_raster_row(slope_fd, slp_raster, data_type);
963
964 if (pcurv_fd > 0)
965 G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
966
967 if (tcurv_fd > 0)
968 G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
969
970 if (dx_fd > 0)
971 G_put_raster_row(dx_fd, dx_raster, data_type);
972
973 if (dy_fd > 0)
974 G_put_raster_row(dy_fd, dy_raster, data_type);
975
976 if (dxx_fd > 0)
977 G_put_raster_row(dxx_fd, dxx_raster, data_type);
978
979 if (dyy_fd > 0)
980 G_put_raster_row(dyy_fd, dyy_raster, data_type);
981
982 if (dxy_fd > 0)
983 G_put_raster_row(dxy_fd, dxy_raster, data_type);
984
985 } /* row loop */
986
987 G_percent(row, nrows, 2);
988
989 G_close_cell(elevation_fd);
990 G_debug(1, "Creating support files...");
991
992 G_verbose_message(_("Elevation products for mapset <%s> in <%s>"),
993 G_mapset(), G_location());
994
995 if (aspect_fd >= 0) {
996 DCELL min, max;
997 struct FPRange range;
998
999 G_set_null_value(asp_raster, G_window_cols(), data_type);
1000 G_put_raster_row(aspect_fd, asp_raster, data_type);
1001 G_close_cell(aspect_fd);
1002
1003 if (out_type != CELL_TYPE)
1004 G_quantize_fp_map_range(aspect_name, G_mapset(), 0., 360., 0,
1005 360);
1006
1007 G_read_raster_cats(aspect_name, G_mapset(), &cats);
1008 G_set_raster_cats_title
1009 ("Aspect counterclockwise in degrees from east", &cats);
1010
1011 G_verbose_message(_("Min computed aspect %.4f, max computed aspect %.4f"),
1012 min_asp, max_asp);
1013 /* the categries quant intervals are 1.0 long, plus
1014 we are using reverse order so that the label looked up
1015 for i-.5 is not the one defined for i-.5, i+.5 interval, but
1016 the one defile for i-1.5, i-.5 interval which is added later */
1017 for (i = ceil(max_asp); i >= 1; i--) {
1018 if (i == 360)
1019 sprintf(buf, "east");
1020 else if (i == 360)
1021 sprintf(buf, "east");
1022 else if (i == 45)
1023 sprintf(buf, "north ccw of east");
1024 else if (i == 90)
1025 sprintf(buf, "north");
1026 else if (i == 135)
1027 sprintf(buf, "north ccw of west");
1028 else if (i == 180)
1029 sprintf(buf, "west");
1030 else if (i == 225)
1031 sprintf(buf, "south ccw of west");
1032 else if (i == 270)
1033 sprintf(buf, "south");
1034 else if (i == 315)
1035 sprintf(buf, "south ccw of east");
1036 else
1037 sprintf(buf, "%d degree%s ccw from east", i,
1038 i == 1 ? "" : "s");
1039 if (data_type == CELL_TYPE) {
1040 G_set_cat(i, buf, &cats);
1041 continue;
1042 }
1043 tmp1 = (double)i - .5;
1044 tmp2 = (double)i + .5;
1045 G_set_d_raster_cat(&tmp1, &tmp2, buf, &cats);
1046 }
1047 if (data_type == CELL_TYPE)
1048 G_set_cat(0, "no aspect", &cats);
1049 else {
1050 tmp1 = 0.;
1051 tmp2 = .5;
1052 G_set_d_raster_cat(&tmp1, &tmp2, "no aspect", &cats);
1053 }
1054 G_write_raster_cats(aspect_name, &cats);
1055 G_free_raster_cats(&cats);
1056
1057 /* write colors for aspect file */
1058 G_init_colors(&colors);
1059 G_read_fp_range(aspect_name, G_mapset(), &range);
1060 G_get_fp_range_min_max(&range, &min, &max);
1061 G_make_aspect_fp_colors(&colors, min, max);
1062 G_write_colors(aspect_name, G_mapset(), &colors);
1063
1064 /* writing history file */
1065 G_short_history(aspect_name, "raster", &hist);
1066 G_snprintf(hist.edhist[0], RECORD_LEN, "aspect map elev = %s",
1067 elev_name);
1068 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1069 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1070 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1071 elev_name);
1072 hist.edlinecnt = 3;
1073
1074 G_command_history(&hist);
1075 G_write_history(aspect_name, &hist);
1076
1077 G_message(_("Aspect raster map <%s> complete"), aspect_name);
1078 }
1079
1080 if (slope_fd >= 0) {
1081 /* colortable for slopes */
1082 G_init_colors(&colors);
1083 G_add_color_rule(0, 255, 255, 255, 2, 255, 255, 0, &colors);
1084 G_add_color_rule(2, 255, 255, 0, 5, 0, 255, 0, &colors);
1085 G_add_color_rule(5, 0, 255, 0, 10, 0, 255, 255, &colors);
1086 G_add_color_rule(10, 0, 255, 255, 15, 0, 0, 255, &colors);
1087 G_add_color_rule(15, 0, 0, 255, 30, 255, 0, 255, &colors);
1088 G_add_color_rule(30, 255, 0, 255, 50, 255, 0, 0, &colors);
1089 G_add_color_rule(50, 255, 0, 0, 90, 0, 0, 0, &colors);
1090
1091 G_set_null_value(slp_raster, G_window_cols(), data_type);
1092 G_put_raster_row(slope_fd, slp_raster, data_type);
1093 G_close_cell(slope_fd);
1094
1095 if (out_type != CELL_TYPE) {
1096 /* INCR_BY_ONE
1097 if(deg)
1098 G_quantize_fp_map_range(slope_name, G_mapset(), 0., 90., 1, 91);
1099 else
1100 G_quantize_fp_map_range(slope_name, G_mapset(), min_slp, max_slp,
1101 (CELL) min_slp + 1, (CELL) ceil(max_slp) + 1);
1102 */
1103 G_write_colors(slope_name, G_mapset(), &colors);
1104 if (deg)
1105 G_quantize_fp_map_range(slope_name, G_mapset(), 0., 90., 0,
1106 90);
1107 else /* percent */
1108 G_quantize_fp_map_range(slope_name, G_mapset(), min_slp,
1109 max_slp, (CELL) min_slp,
1110 (CELL) ceil(max_slp));
1111 }
1112
1113 G_read_raster_cats(slope_name, G_mapset(), &cats);
1114 if (deg)
1115 G_set_raster_cats_title("slope in degrees", &cats);
1116 else if (perc)
1117 G_set_raster_cats_title("percent slope", &cats);
1118
1119 G_verbose_message(_("Min computed slope %.4f, max computed slope %.4f"),
1120 min_slp, max_slp);
1121 /* the categries quant intervals are 1.0 long, plus
1122 we are using reverse order so that the label looked up
1123 for i-.5 is not the one defined for i-.5, i+.5 interval, but
1124 the one defined for i-1.5, i-.5 interval which is added later */
1125 for (i = ceil(max_slp); i > /* INC BY ONE >= */ 0; i--) {
1126 if (deg)
1127 sprintf(buf, "%d degree%s", i, i == 1 ? "" : "s");
1128 else if (perc)
1129 sprintf(buf, "%d percent", i);
1130 if (data_type == CELL_TYPE) {
1131 /* INCR_BY_ONE
1132 G_set_cat(i+1, buf, &cats);
1133 */
1134 G_set_cat(i, buf, &cats);
1135 continue;
1136 }
1137 /* INCR_BY_ONE
1138 tmp1 = (DCELL) i+.5;
1139 tmp2 = (DCELL) i+1.5;
1140 */
1141 tmp1 = (DCELL) i - .5;
1142 tmp2 = (DCELL) i + .5;
1143 G_set_d_raster_cat(&tmp1, &tmp2, buf, &cats);
1144 }
1145 if (data_type == CELL_TYPE)
1146 G_set_cat(0, "zero slope", &cats);
1147 /* INCR_BY_ONE
1148 G_set_cat(0, "no data", &cats);
1149 */
1150 else {
1151 tmp1 = 0;
1152 tmp2 = 0.5;
1153 G_set_d_raster_cat(&tmp1, &tmp2, "zero slope", &cats);
1154 }
1155 /* INCR_BY_ONE
1156 G_set_d_raster_cat (&tmp1, &tmp1, "no data", &cats);
1157 */
1158 G_write_raster_cats(slope_name, &cats);
1159
1160 /* writing history file */
1161 G_short_history(slope_name, "raster", &hist);
1162 G_snprintf(hist.edhist[0], RECORD_LEN, "slope map elev = %s",
1163 elev_name);
1164 sprintf(hist.edhist[1], "zfactor = %.2f format = %s", zfactor,
1165 parm.slope_fmt->answer);
1166 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1167 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1168 elev_name);
1169 hist.edlinecnt = 3;
1170
1171 G_command_history(&hist);
1172 G_write_history(slope_name, &hist);
1173
1174 G_message(_("Slope raster map <%s> complete"), slope_name);
1175 }
1176
1177 /* colortable for curvatures */
1178 if (pcurv_fd >= 0 || tcurv_fd >= 0) {
1179 G_init_colors(&colors);
1180 if (c1min < c2min)
1181 dat1 = (FCELL) c1min;
1182 else
1183 dat1 = (FCELL) c2min;
1184
1185 dat2 = (FCELL) - 0.01;
1186 G_add_f_raster_color_rule(&dat1, 127, 0, 255,
1187 &dat2, 0, 0, 255, &colors);
1188 dat1 = dat2;
1189 dat2 = (FCELL) - 0.001;
1190 G_add_f_raster_color_rule(&dat1, 0, 0, 255,
1191 &dat2, 0, 127, 255, &colors);
1192 dat1 = dat2;
1193 dat2 = (FCELL) - 0.00001;
1194 G_add_f_raster_color_rule(&dat1, 0, 127, 255,
1195 &dat2, 0, 255, 255, &colors);
1196 dat1 = dat2;
1197 dat2 = (FCELL) 0.0;
1198 G_add_f_raster_color_rule(&dat1, 0, 255, 255,
1199 &dat2, 200, 255, 200, &colors);
1200 dat1 = dat2;
1201 dat2 = (FCELL) 0.00001;
1202 G_add_f_raster_color_rule(&dat1, 200, 255, 200,
1203 &dat2, 255, 255, 0, &colors);
1204 dat1 = dat2;
1205 dat2 = (FCELL) 0.001;
1206 G_add_f_raster_color_rule(&dat1, 255, 255, 0,
1207 &dat2, 255, 127, 0, &colors);
1208 dat1 = dat2;
1209 dat2 = (FCELL) 0.01;
1210 G_add_f_raster_color_rule(&dat1, 255, 127, 0,
1211 &dat2, 255, 0, 0, &colors);
1212 dat1 = dat2;
1213 if (c1max > c2max)
1214 dat2 = (FCELL) c1max;
1215 else
1216 dat2 = (FCELL) c2max;
1217
1218 G_add_f_raster_color_rule(&dat1, 255, 0, 0,
1219 &dat2, 255, 0, 200, &colors);
1220 }
1221
1222 if (pcurv_fd >= 0) {
1223 G_set_null_value(pcurv_raster, G_window_cols(), data_type);
1224 G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
1225 G_close_cell(pcurv_fd);
1226
1227 G_write_colors(pcurv_name, G_mapset(), &colors);
1228
1229 if (out_type != CELL_TYPE)
1230 G_round_fp_map(pcurv_name, G_mapset());
1231
1232 G_read_cats(pcurv_name, G_mapset(), &cats);
1233 G_set_cats_title("profile curvature", &cats);
1234 G_set_cat((CELL) 0, "no profile curve", &cats);
1235
1236 /* writing history file */
1237 G_short_history(pcurv_name, "raster", &hist);
1238 G_snprintf(hist.edhist[0], RECORD_LEN, "profile curve map elev = %s",
1239 elev_name);
1240 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1241 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1242 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1243 elev_name);
1244 hist.edlinecnt = 3;
1245
1246 G_command_history(&hist);
1247 G_write_history(pcurv_name, &hist);
1248
1249 G_message(_("Profile curve raster map <%s> complete"), pcurv_name);
1250 }
1251
1252 if (tcurv_fd >= 0) {
1253 G_set_null_value(tcurv_raster, G_window_cols(), data_type);
1254 G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
1255 G_close_cell(tcurv_fd);
1256
1257 G_write_colors(tcurv_name, G_mapset(), &colors);
1258
1259 if (out_type != CELL_TYPE)
1260 G_round_fp_map(tcurv_name, G_mapset());
1261
1262 G_read_cats(tcurv_name, G_mapset(), &cats);
1263 G_set_cats_title("tangential curvature", &cats);
1264 G_set_cat((CELL) 0, "no tangential curve", &cats);
1265
1266 /* writing history file */
1267 G_short_history(tcurv_name, "raster", &hist);
1268 G_snprintf(hist.edhist[0], RECORD_LEN, "tangential curve map elev = %s",
1269 elev_name);
1270 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1271 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1272 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1273 elev_name);
1274 hist.edlinecnt = 3;
1275
1276 G_command_history(&hist);
1277 G_write_history(tcurv_name, &hist);
1278
1279 G_message(_("Tangential curve raster map <%s> complete"), tcurv_name);
1280 }
1281
1282 if (dx_fd >= 0) {
1283 G_set_null_value(dx_raster, G_window_cols(), data_type);
1284 G_put_raster_row(dx_fd, dx_raster, data_type);
1285 G_close_cell(dx_fd);
1286
1287 if (out_type != CELL_TYPE)
1288 G_round_fp_map(dx_name, G_mapset());
1289
1290 G_read_cats(dx_name, G_mapset(), &cats);
1291 G_set_cats_title("E-W slope", &cats);
1292 G_set_cat((CELL) 0, "no E-W slope", &cats);
1293
1294 /* writing history file */
1295 G_short_history(dx_name, "raster", &hist);
1296 G_snprintf(hist.edhist[0], RECORD_LEN, "E-W slope map elev = %s",
1297 elev_name);
1298 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1299 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1300 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1301 elev_name);
1302 hist.edlinecnt = 3;
1303
1304 G_command_history(&hist);
1305 G_write_history(dx_name, &hist);
1306
1307 G_message(_("E-W slope raster map <%s> complete"), dx_name);
1308 }
1309
1310 if (dy_fd >= 0) {
1311 G_set_null_value(dy_raster, G_window_cols(), data_type);
1312 G_put_raster_row(dy_fd, dy_raster, data_type);
1313 G_close_cell(dy_fd);
1314
1315 if (out_type != CELL_TYPE)
1316 G_round_fp_map(dy_name, G_mapset());
1317
1318 G_read_cats(dy_name, G_mapset(), &cats);
1319 G_set_cats_title("N-S slope", &cats);
1320 G_set_cat((CELL) 0, "no N-S slope", &cats);
1321
1322 /* writing history file */
1323 G_short_history(dy_name, "raster", &hist);
1324 G_snprintf(hist.edhist[0], RECORD_LEN, "N-S slope map elev = %s",
1325 elev_name);
1326 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1327 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1328 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1329 elev_name);
1330 hist.edlinecnt = 3;
1331
1332 G_command_history(&hist);
1333 G_write_history(dy_name, &hist);
1334
1335 G_message(_("N-S slope raster map <%s> complete"), dy_name);
1336 }
1337
1338 if (dxx_fd >= 0) {
1339 G_set_null_value(dxx_raster, G_window_cols(), data_type);
1340 G_put_raster_row(dxx_fd, dxx_raster, data_type);
1341 G_close_cell(dxx_fd);
1342
1343 if (out_type != CELL_TYPE)
1344 G_round_fp_map(dxx_name, G_mapset());
1345
1346 G_read_cats(dxx_name, G_mapset(), &cats);
1347 G_set_cats_title("DXX", &cats);
1348 G_set_cat((CELL) 0, "DXX", &cats);
1349
1350 /* writing history file */
1351 G_short_history(dxx_name, "raster", &hist);
1352 G_snprintf(hist.edhist[0], RECORD_LEN, "DXX map elev = %s",
1353 elev_name);
1354 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1355 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1356 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1357 elev_name);
1358 hist.edlinecnt = 3;
1359
1360 G_command_history(&hist);
1361 G_write_history(dxx_name, &hist);
1362
1363 G_message(_("Dxx raster map <%s> complete"), dxx_name);
1364 }
1365
1366 if (dyy_fd >= 0) {
1367 G_set_null_value(dyy_raster, G_window_cols(), data_type);
1368 G_put_raster_row(dyy_fd, dyy_raster, data_type);
1369 G_close_cell(dyy_fd);
1370
1371 if (out_type != CELL_TYPE)
1372 G_round_fp_map(dyy_name, G_mapset());
1373
1374 G_read_cats(dyy_name, G_mapset(), &cats);
1375 G_set_cats_title("DYY", &cats);
1376 G_set_cat((CELL) 0, "DYY", &cats);
1377
1378 /* writing history file */
1379 G_short_history(dyy_name, "raster", &hist);
1380 G_snprintf(hist.edhist[0], RECORD_LEN, "DYY map elev = %s",
1381 elev_name);
1382 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1383 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1384 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1385 elev_name);
1386 hist.edlinecnt = 3;
1387
1388 G_command_history(&hist);
1389 G_write_history(dyy_name, &hist);
1390
1391 G_message(_("Dyy raster map <%s> complete"), dyy_name);
1392 }
1393
1394 if (dxy_fd >= 0) {
1395 G_set_null_value(dxy_raster, G_window_cols(), data_type);
1396 G_put_raster_row(dxy_fd, dxy_raster, data_type);
1397 G_close_cell(dxy_fd);
1398
1399 if (out_type != CELL_TYPE)
1400 G_round_fp_map(dxy_name, G_mapset());
1401
1402 G_read_cats(dxy_name, G_mapset(), &cats);
1403 G_set_cats_title("DXY", &cats);
1404 G_set_cat((CELL) 0, "DXY", &cats);
1405
1406 /* writing history file */
1407 G_short_history(dxy_name, "raster", &hist);
1408 G_snprintf(hist.edhist[0], RECORD_LEN, "DXY map elev = %s",
1409 elev_name);
1410 sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1411 sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1412 G_snprintf(hist.datsrc_1, RECORD_LEN, "raster elevation file %s",
1413 elev_name);
1414 hist.edlinecnt = 3;
1415
1416 G_command_history(&hist);
1417 G_write_history(dxy_name, &hist);
1418
1419 G_message(_("Dxy raster map <%s> complete"), dxy_name);
1420 }
1421
1422 exit(EXIT_SUCCESS);
1423}
Note: See TracBrowser for help on using the repository browser.