2
/****************************************************************************
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_nospam yahoo.com>,
12
* Jachym Cepicky <jachym les-ejk.cz>,
13
* Jan-Oliver Wagner <jan intevation.de>,
14
* Radim Blazek <radim.blazek gmail.com>
15
* PURPOSE: generates raster maps of slope, aspect, curvatures and
16
* first and second order partial derivatives from a raster map
17
* of true elevation values
18
* COPYRIGHT: (C) 1999-2006 by the GRASS Development Team
20
* This program is free software under the GNU General Public
21
* License (>=v2). Read the file COPYING that comes with GRASS
24
*****************************************************************************/
28
#include <grass/gis.h>
29
#include <grass/glocale.h>
30
#include "local_proto.h"
32
/* 10/99 from GMSL, updated to new GRASS 5 code style , changed default "prec" to float */
35
#define abs(x) ((x)<0?-(x):(x))
38
/**************************************************************************
39
* input is from command line.
40
* arguments are elevation file, slope file, aspect file, profile curvature
41
* file and tangential curvature file
42
* elevation filename required
43
* either slope filename or aspect filename or profile curvature filename
44
* or tangential curvature filename required
45
* usage: r.slope.aspect [-av] elevation=input slope=output1 aspect=output2
46
* pcurv=output3 tcurv=output4 format=name prec=name zfactor=value
47
* min_slp_allowed=value dx=output5 dy=output6 dxx=output7
48
* dyy=output8 dxy=output9
49
* -a don't align window
51
**************************************************************************/
53
/* some changes made to code to retrieve correct distances when using
54
lat/lon projection. changes involve recalculating H and V. see
55
comments within code. */
56
/* added colortables for topographic parameters and fixed
57
* the sign bug for the second order derivatives (jh) */
60
int r_slope_aspect(int argc, char *argv[])
62
struct Categories cats;
73
DCELL *elev_cell[3], *temp;
74
DCELL *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9;
77
void *asp_raster, *asp_ptr = NULL;
78
void *slp_raster, *slp_ptr = NULL;
79
void *pcurv_raster, *pcurv_ptr = NULL;
80
void *tcurv_raster, *tcurv_ptr = NULL;
81
void *dx_raster, *dx_ptr = NULL;
82
void *dy_raster, *dy_ptr = NULL;
83
void *dxx_raster, *dxx_ptr = NULL;
84
void *dyy_raster, *dyy_ptr = NULL;
85
void *dxy_raster, *dxy_ptr = NULL;
87
RASTER_MAP_TYPE out_type = 0, data_type;
88
int Wrap; /* global wraparound */
89
struct Cell_head window, cellhd;
108
double north, east, south, west, ns_med;
110
double radians_to_degrees;
111
double degrees_to_radians;
113
double dx; /* partial derivative in ew direction */
114
double dy; /* partial derivative in ns direction */
115
double dxx, dxy, dyy;
116
double s3, s4, s5, s6;
118
double scik1 = 100000.;
121
double aspect, min_asp = 360., max_asp = 0.;
122
double dnorm1, dx2, dy2, grad2, grad, dxy2;
123
double gradmin = 0.001;
124
double c1min = 0., c1max = 0., c2min = 0., c2max = 0.;
130
double slp_in_perc, slp_in_deg;
131
double min_slp = 900., max_slp = 0., min_slp_allowed;
132
int low, hi, test = 0;
136
struct GModule *module;
139
struct Option *elevation, *slope_fmt, *slope, *aspect, *pcurv, *tcurv,
140
*zfactor, *min_slp_allowed, *out_precision,
141
*dx, *dy, *dxx, *dyy, *dxy;
146
/* please, remove before GRASS 7 released */
152
module = G_define_module();
153
module->keywords = _("raster");
154
module->description =
155
_("Generates raster map layers of slope, aspect, curvatures and "
156
"partial derivatives from a raster map layer of true elevation "
157
"values. Aspect is calculated counterclockwise from east.");
159
parm.elevation = G_define_option();
160
parm.elevation->key = "elevation";
161
parm.elevation->type = TYPE_STRING;
162
parm.elevation->required = YES;
163
parm.elevation->gisprompt = "old,cell,raster";
164
parm.elevation->description = _("Raster elevation file name");
166
parm.slope = G_define_option();
167
parm.slope->key = "slope";
168
parm.slope->type = TYPE_STRING;
169
parm.slope->required = NO;
170
parm.slope->answer = NULL;
171
parm.slope->gisprompt = "new,cell,raster";
172
parm.slope->description = _("Output slope filename");
174
parm.aspect = G_define_option();
175
parm.aspect->key = "aspect";
176
parm.aspect->type = TYPE_STRING;
177
parm.aspect->required = NO;
178
parm.aspect->answer = NULL;
179
parm.aspect->gisprompt = "new,cell,raster";
180
parm.aspect->description = _("Output aspect filename");
182
parm.slope_fmt = G_define_option();
183
parm.slope_fmt->key = "format";
184
parm.slope_fmt->type = TYPE_STRING;
185
parm.slope_fmt->required = NO;
186
parm.slope_fmt->answer = "degrees";
187
parm.slope_fmt->options = "degrees,percent";
188
parm.slope_fmt->description = _("Format for reporting the slope");
189
parm.slope_fmt->guisection = _("Settings");
191
parm.out_precision = G_define_option();
192
parm.out_precision->key = "prec";
193
parm.out_precision->type = TYPE_STRING;
194
parm.out_precision->required = NO;
195
parm.out_precision->answer = "float";
196
parm.out_precision->options = "default,double,float,int";
197
parm.out_precision->description =
198
_("Type of output aspect and slope maps");
199
parm.out_precision->guisection = _("Settings");
201
parm.pcurv = G_define_option();
202
parm.pcurv->key = "pcurv";
203
parm.pcurv->type = TYPE_STRING;
204
parm.pcurv->required = NO;
205
parm.pcurv->answer = NULL;
206
parm.pcurv->gisprompt = "new,cell,raster";
207
parm.pcurv->description = _("Output profile curvature filename");
208
parm.pcurv->guisection = _("Advanced");
210
parm.tcurv = G_define_option();
211
parm.tcurv->key = "tcurv";
212
parm.tcurv->type = TYPE_STRING;
213
parm.tcurv->required = NO;
214
parm.tcurv->answer = NULL;
215
parm.tcurv->gisprompt = "new,cell,raster";
216
parm.tcurv->description = _("Output tangential curvature filename");
217
parm.tcurv->guisection = _("Advanced");
219
parm.dx = G_define_option();
221
parm.dx->type = TYPE_STRING;
222
parm.dx->required = NO;
223
parm.dx->answer = NULL;
224
parm.dx->gisprompt = "new,cell,raster";
225
parm.dx->description =
226
_("Output first order partial derivative dx (E-W slope) filename");
227
parm.dx->guisection = _("Advanced");
229
parm.dy = G_define_option();
231
parm.dy->type = TYPE_STRING;
232
parm.dy->required = NO;
233
parm.dy->answer = NULL;
234
parm.dy->gisprompt = "new,cell,raster";
235
parm.dy->description =
236
_("Output first order partial derivative dy (N-S slope) filename");
237
parm.dy->guisection = _("Advanced");
239
parm.dxx = G_define_option();
240
parm.dxx->key = "dxx";
241
parm.dxx->type = TYPE_STRING;
242
parm.dxx->required = NO;
243
parm.dxx->answer = NULL;
244
parm.dxx->gisprompt = "new,cell,raster";
245
parm.dxx->description =
246
_("Output second order partial derivative dxx filename");
247
parm.dxx->guisection = _("Advanced");
249
parm.dyy = G_define_option();
250
parm.dyy->key = "dyy";
251
parm.dyy->type = TYPE_STRING;
252
parm.dyy->required = NO;
253
parm.dyy->answer = NULL;
254
parm.dyy->gisprompt = "new,cell,raster";
255
parm.dyy->description =
256
_("Output second order partial derivative dyy filename");
257
parm.dyy->guisection = _("Advanced");
259
parm.dxy = G_define_option();
260
parm.dxy->key = "dxy";
261
parm.dxy->type = TYPE_STRING;
262
parm.dxy->required = NO;
263
parm.dxy->answer = NULL;
264
parm.dxy->gisprompt = "new,cell,raster";
265
parm.dxy->description =
266
_("Output second order partial derivative dxy filename");
267
parm.dxy->guisection = _("Advanced");
269
parm.zfactor = G_define_option();
270
parm.zfactor->key = "zfactor";
271
parm.zfactor->description =
272
_("Multiplicative factor to convert elevation units to meters");
273
parm.zfactor->type = TYPE_DOUBLE;
274
parm.zfactor->required = NO;
275
parm.zfactor->answer = "1.0";
276
parm.zfactor->guisection = _("Settings");
278
parm.min_slp_allowed = G_define_option();
279
parm.min_slp_allowed->key = "min_slp_allowed";
280
parm.min_slp_allowed->description =
281
_("Minimum slope val. (in percent) for which aspect is computed");
282
parm.min_slp_allowed->type = TYPE_DOUBLE;
283
parm.min_slp_allowed->required = NO;
284
parm.min_slp_allowed->answer = "0.0";
285
parm.min_slp_allowed->guisection = _("Settings");
287
/* please, remove before GRASS 7 released */
288
flag.q = G_define_flag();
290
flag.q->description = _("Quiet");
292
flag.a = G_define_flag();
294
flag.a->description =
295
_("Do not align the current region to the elevation layer");
296
flag.a->guisection = _("Settings");
299
radians_to_degrees = 180.0 / M_PI;
300
degrees_to_radians = M_PI / 180.0;
304
answer[91] = 15000.0;
306
for (i = 1; i < 91; i++)
309
tan_ans = tan ( degrees / radians_to_degrees );
310
answer[i] = tan_ans * tan_ans;
314
answer[90] = 15000.0;
316
for (i = 0; i < 90; i++) {
318
tan_ans = tan(degrees / radians_to_degrees);
319
answer[i] = tan_ans * tan_ans;
322
if (G_parser(argc, argv))
325
/* please, remove before GRASS 7 released */
326
if (flag.q->answer) {
327
putenv("GRASS_VERBOSE=0");
328
G_warning(_("The '-q' flag is superseded and will be removed "
329
"in future. Please use '--quiet' instead."));
334
elev_name = parm.elevation->answer;
335
slope_name = parm.slope->answer;
336
aspect_name = parm.aspect->answer;
337
pcurv_name = parm.pcurv->answer;
338
tcurv_name = parm.tcurv->answer;
339
dx_name = parm.dx->answer;
340
dy_name = parm.dy->answer;
341
dxx_name = parm.dxx->answer;
342
dyy_name = parm.dyy->answer;
343
dxy_name = parm.dxy->answer;
345
G_check_input_output_name(elev_name, slope_name, GR_FATAL_EXIT);
346
G_check_input_output_name(elev_name, aspect_name, GR_FATAL_EXIT);
347
G_check_input_output_name(elev_name, pcurv_name, GR_FATAL_EXIT);
348
G_check_input_output_name(elev_name, tcurv_name, GR_FATAL_EXIT);
349
G_check_input_output_name(elev_name, dx_name, GR_FATAL_EXIT);
350
G_check_input_output_name(elev_name, dy_name, GR_FATAL_EXIT);
351
G_check_input_output_name(elev_name, dxx_name, GR_FATAL_EXIT);
352
G_check_input_output_name(elev_name, dyy_name, GR_FATAL_EXIT);
353
G_check_input_output_name(elev_name, dxy_name, GR_FATAL_EXIT);
355
if (sscanf(parm.zfactor->answer, "%lf", &zfactor) != 1 || zfactor <= 0.0) {
356
G_warning("%s=%s - must be a postive number", parm.zfactor->key,
357
parm.zfactor->answer);
362
if (sscanf(parm.min_slp_allowed->answer, "%lf", &min_slp_allowed) != 1 ||
363
min_slp_allowed < 0.0) {
364
G_warning("%s=%s - must be a non-negative number",
365
parm.min_slp_allowed->key, parm.min_slp_allowed->answer);
370
slope_fmt = parm.slope_fmt->answer;
371
if (strcmp(slope_fmt, "percent") == 0)
373
else if (strcmp(slope_fmt, "degrees") == 0)
376
if (slope_name == NULL && aspect_name == NULL
377
&& pcurv_name == NULL && tcurv_name == NULL
378
&& dx_name == NULL && dy_name == NULL
379
&& dxx_name == NULL && dyy_name == NULL && dxy_name == NULL) {
380
G_warning("You must specify at least one of the parameters:"
381
"\n<%s>, <%s>, <%s>, <%s>, <%s>, <%s>, <%s>, <%s> or <%s>.\n",
382
parm.slope->key, parm.aspect->key, parm.pcurv->key,
383
parm.tcurv->key, parm.dx->key, parm.dy->key,
384
parm.dxx->key, parm.dyy->key, parm.dxy->key);
389
/* check elevation file existence */
390
mapset = G_find_cell2(elev_name, "");
392
G_fatal_error(_("Raster map <%s> not found"), elev_name);
394
/* set the window from the header for the elevation file */
395
if (!flag.a->answer) {
396
G_get_window(&window);
397
if (G_get_cellhd(elev_name, mapset, &cellhd) >= 0) {
398
G_align_window(&window, &cellhd);
399
G_set_window(&window);
403
if (strcmp(parm.out_precision->answer, "double") == 0)
404
out_type = DCELL_TYPE;
405
else if (strcmp(parm.out_precision->answer, "float") == 0)
406
out_type = FCELL_TYPE;
407
else if (strcmp(parm.out_precision->answer, "int") == 0)
408
out_type = CELL_TYPE;
409
else if (strcmp(parm.out_precision->answer, "default") == 0)
412
G_fatal_error(_("wrong type: %s"), parm.out_precision->answer);
414
data_type = out_type;
416
data_type = DCELL_TYPE;
417
/* data type is the type of data being processed,
418
out_type is type of map being created */
420
G_get_set_window(&window);
422
nrows = G_window_rows();
423
ncols = G_window_cols();
425
if (((window.west == (window.east - 360.))
426
|| (window.east == (window.west - 360.))) &&
427
(G_projection() == PROJECTION_LL)) {
434
/* H = window.ew_res * 4 * 2/ zfactor; *//* horizontal (east-west) run
435
times 4 for weighted difference */
436
/* V = window.ns_res * 4 * 2/ zfactor; *//* vertical (north-south) run
437
times 4 for weighted difference */
439
/* give warning if location units are different from meters and zfactor=1 */
440
factor = G_database_units_to_meters_factor();
442
G_warning("Converting units to meters, factor=%.6f", factor);
444
G_begin_distance_calculations();
445
north = G_row_to_northing(0.5, &window);
446
ns_med = G_row_to_northing(1.5, &window);
447
south = G_row_to_northing(2.5, &window);
448
east = G_col_to_easting(2.5, &window);
449
west = G_col_to_easting(0.5, &window);
450
V = G_distance(east, north, east, south) * 4 / zfactor;
451
H = G_distance(east, ns_med, west, ns_med) * 4 / zfactor;
452
/* ____________________________
457
|________|________|________|
460
| east | ns_med | west |
462
|________|________|________|
467
|________|________|________|
470
/* open the elevation file for reading */
471
elevation_fd = G_open_cell_old(elev_name, mapset);
472
if (elevation_fd < 0)
474
elev_cell[0] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
475
G_set_d_null_value(elev_cell[0], ncols);
476
elev_cell[1] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
477
G_set_d_null_value(elev_cell[1], ncols);
478
elev_cell[2] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
479
G_set_d_null_value(elev_cell[2], ncols);
481
if (slope_name != NULL) {
482
slope_fd = opennew(slope_name, out_type);
483
slp_raster = G_allocate_raster_buf(data_type);
484
G_set_null_value(slp_raster, G_window_cols(), data_type);
485
G_put_raster_row(slope_fd, slp_raster, data_type);
492
if (aspect_name != NULL) {
493
aspect_fd = opennew(aspect_name, out_type);
494
asp_raster = G_allocate_raster_buf(data_type);
495
G_set_null_value(asp_raster, G_window_cols(), data_type);
496
G_put_raster_row(aspect_fd, asp_raster, data_type);
503
if (pcurv_name != NULL) {
504
pcurv_fd = opennew(pcurv_name, out_type);
505
pcurv_raster = G_allocate_raster_buf(data_type);
506
G_set_null_value(pcurv_raster, G_window_cols(), data_type);
507
G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
514
if (tcurv_name != NULL) {
515
tcurv_fd = opennew(tcurv_name, out_type);
516
tcurv_raster = G_allocate_raster_buf(data_type);
517
G_set_null_value(tcurv_raster, G_window_cols(), data_type);
518
G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
525
if (dx_name != NULL) {
526
dx_fd = opennew(dx_name, out_type);
527
dx_raster = G_allocate_raster_buf(data_type);
528
G_set_null_value(dx_raster, G_window_cols(), data_type);
529
G_put_raster_row(dx_fd, dx_raster, data_type);
536
if (dy_name != NULL) {
537
dy_fd = opennew(dy_name, out_type);
538
dy_raster = G_allocate_raster_buf(data_type);
539
G_set_null_value(dy_raster, G_window_cols(), data_type);
540
G_put_raster_row(dy_fd, dy_raster, data_type);
547
if (dxx_name != NULL) {
548
dxx_fd = opennew(dxx_name, out_type);
549
dxx_raster = G_allocate_raster_buf(data_type);
550
G_set_null_value(dxx_raster, G_window_cols(), data_type);
551
G_put_raster_row(dxx_fd, dxx_raster, data_type);
558
if (dyy_name != NULL) {
559
dyy_fd = opennew(dyy_name, out_type);
560
dyy_raster = G_allocate_raster_buf(data_type);
561
G_set_null_value(dyy_raster, G_window_cols(), data_type);
562
G_put_raster_row(dyy_fd, dyy_raster, data_type);
569
if (dxy_name != NULL) {
570
dxy_fd = opennew(dxy_name, out_type);
571
dxy_raster = G_allocate_raster_buf(data_type);
572
G_set_null_value(dxy_raster, G_window_cols(), data_type);
573
G_put_raster_row(dxy_fd, dxy_raster, data_type);
580
if (aspect_fd < 0 && slope_fd < 0 && pcurv_fd < 0 && tcurv_fd < 0
581
&& dx_fd < 0 && dy_fd < 0 && dxx_fd < 0 && dyy_fd < 0 && dxy_fd < 0)
585
G_get_d_raster_row_nomask(elevation_fd, elev_cell[1] + 1, 0);
586
elev_cell[1][0] = elev_cell[1][G_window_cols() - 1];
587
elev_cell[1][G_window_cols() + 1] = elev_cell[1][2];
590
G_get_d_raster_row_nomask(elevation_fd, elev_cell[1], 0);
593
G_get_d_raster_row_nomask(elevation_fd, elev_cell[2] + 1, 1);
594
elev_cell[2][0] = elev_cell[2][G_window_cols() - 1];
595
elev_cell[2][G_window_cols() + 1] = elev_cell[2][2];
598
G_get_d_raster_row_nomask(elevation_fd, elev_cell[2], 1);
600
G_message(_("Percent complete: "));
601
for (row = 2; row < nrows; row++) {
602
/* if projection is Lat/Lon, recalculate V and H */
603
if (G_projection() == PROJECTION_LL) {
604
north = G_row_to_northing((row - 2 + 0.5), &window);
605
ns_med = G_row_to_northing((row - 1 + 0.5), &window);
606
south = G_row_to_northing((row + 0.5), &window);
607
east = G_col_to_easting(2.5, &window);
608
west = G_col_to_easting(0.5, &window);
609
V = G_distance(east, north, east, south) * 4 / zfactor;
610
H = G_distance(east, ns_med, west, ns_med) * 4 / zfactor;
611
/* ____________________________
616
|________|________|________|
619
| east | ns_med | west |
621
|________|________|________|
626
|________|________|________|
630
G_percent(row, nrows, 2);
632
elev_cell[0] = elev_cell[1];
633
elev_cell[1] = elev_cell[2];
637
G_get_d_raster_row_nomask(elevation_fd, elev_cell[2] + 1, row);
638
elev_cell[2][0] = elev_cell[2][G_window_cols() - 1];
639
elev_cell[2][G_window_cols() + 1] = elev_cell[2][2];
642
G_get_d_raster_row_nomask(elevation_fd, elev_cell[2], row);
654
if (aspect_fd >= 0) {
656
asp_ptr = asp_raster;
659
G_incr_void_ptr(asp_raster, G_raster_size(data_type));
663
slp_ptr = slp_raster;
666
G_incr_void_ptr(slp_raster, G_raster_size(data_type));
671
pcurv_ptr = pcurv_raster;
674
G_incr_void_ptr(pcurv_raster, G_raster_size(data_type));
679
tcurv_ptr = tcurv_raster;
682
G_incr_void_ptr(tcurv_raster, G_raster_size(data_type));
689
dx_ptr = G_incr_void_ptr(dx_raster, G_raster_size(data_type));
696
dy_ptr = G_incr_void_ptr(dy_raster, G_raster_size(data_type));
701
dxx_ptr = dxx_raster;
704
G_incr_void_ptr(dxx_raster, G_raster_size(data_type));
709
dyy_ptr = dyy_raster;
712
G_incr_void_ptr(dyy_raster, G_raster_size(data_type));
717
dxy_ptr = dxy_raster;
720
G_incr_void_ptr(dxy_raster, G_raster_size(data_type));
724
/*skip first cell of the row */
726
for (col = ncols - 2; col-- > 0;
727
c1++, c2++, c3++, c4++, c5++, c6++, c7++, c8++, c9++) {
729
fprintf(stdout, "\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n%.0f %.0f %.0f\n",
730
*c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9);
733
if (G_is_d_null_value(c1) || G_is_d_null_value(c2) ||
734
G_is_d_null_value(c3) || G_is_d_null_value(c4) ||
735
G_is_d_null_value(c5) || G_is_d_null_value(c6) ||
736
G_is_d_null_value(c7) || G_is_d_null_value(c8) ||
737
G_is_d_null_value(c9)) {
739
G_set_null_value(slp_ptr, 1, data_type);
741
G_incr_void_ptr(slp_ptr, G_raster_size(data_type));
744
G_set_null_value(asp_ptr, 1, data_type);
746
G_incr_void_ptr(asp_ptr, G_raster_size(data_type));
749
G_set_null_value(pcurv_ptr, 1, data_type);
751
G_incr_void_ptr(pcurv_ptr, G_raster_size(data_type));
754
G_set_null_value(tcurv_ptr, 1, data_type);
756
G_incr_void_ptr(tcurv_ptr, G_raster_size(data_type));
759
G_set_null_value(dx_ptr, 1, data_type);
761
G_incr_void_ptr(dx_ptr, G_raster_size(data_type));
764
G_set_null_value(dy_ptr, 1, data_type);
766
G_incr_void_ptr(dy_ptr, G_raster_size(data_type));
769
G_set_null_value(dxx_ptr, 1, data_type);
771
G_incr_void_ptr(dxx_ptr, G_raster_size(data_type));
774
G_set_null_value(dyy_ptr, 1, data_type);
776
G_incr_void_ptr(dyy_ptr, G_raster_size(data_type));
779
G_set_null_value(dxy_ptr, 1, data_type);
781
G_incr_void_ptr(dxy_ptr, G_raster_size(data_type));
786
dx = ((*c1 + *c4 + *c4 + *c7) - (*c3 + *c6 + *c6 + *c9)) / H;
787
dy = ((*c7 + *c8 + *c8 + *c9) - (*c1 + *c2 + *c2 + *c3)) / V;
789
/* compute topographic parameters */
790
key = dx * dx + dy * dy;
791
slp_in_perc = 100 * sqrt(key);
792
slp_in_deg = atan(sqrt(key)) * radians_to_degrees;
794
/* now update min and max */
796
if (min_slp > slp_in_deg)
797
min_slp = slp_in_deg;
798
if (max_slp < slp_in_deg)
799
max_slp = slp_in_deg;
802
if (min_slp > slp_in_perc)
803
min_slp = slp_in_perc;
804
if (max_slp < slp_in_perc)
805
max_slp = slp_in_perc;
807
if (slp_in_perc < min_slp_allowed)
810
if (deg && out_type == CELL_TYPE) {
820
if (key >= answer[test])
822
else if (key < answer[test - 1])
826
test = (low + hi) / 2;
829
else if (perc && out_type == CELL_TYPE)
831
/* test = slp_in_perc + 1.5; *//* All the slope categories are
833
test = slp_in_perc + .5;
836
if (data_type == CELL_TYPE)
837
*((CELL *) slp_ptr) = (CELL) test;
840
G_set_raster_value_d(slp_ptr,
841
(DCELL) slp_in_deg, data_type);
843
G_set_raster_value_d(slp_ptr,
844
(DCELL) slp_in_perc, data_type);
846
slp_ptr = G_incr_void_ptr(slp_ptr, G_raster_size(data_type));
847
} /* computing slope */
859
aspect = (atan2(dy, dx) / degrees_to_radians);
860
if ((aspect <= 0.5) && (aspect > 0) &&
861
out_type == CELL_TYPE)
864
aspect = 360. + aspect;
867
/* if it's not the case that the slope for this cell
868
is below specified minimum */
869
if (!((slope_fd > 0) && (slp_in_perc < min_slp_allowed))) {
870
if (out_type == CELL_TYPE)
871
*((CELL *) asp_ptr) = (CELL) (aspect + .5);
873
G_set_raster_value_d(asp_ptr,
874
(DCELL) aspect, data_type);
877
G_set_null_value(asp_ptr, 1, data_type);
878
asp_ptr = G_incr_void_ptr(asp_ptr, G_raster_size(data_type));
880
/* now update min and max */
881
if (min_asp > aspect)
883
if (max_asp < aspect)
885
} /* computing aspect */
888
if (out_type == CELL_TYPE)
889
*((CELL *) dx_ptr) = (CELL) (scik1 * dx);
891
G_set_raster_value_d(dx_ptr, (DCELL) dx, data_type);
892
dx_ptr = G_incr_void_ptr(dx_ptr, G_raster_size(data_type));
896
if (out_type == CELL_TYPE)
897
*((CELL *) dy_ptr) = (CELL) (scik1 * dy);
899
G_set_raster_value_d(dy_ptr, (DCELL) dy, data_type);
900
dy_ptr = G_incr_void_ptr(dy_ptr, G_raster_size(data_type));
903
if (dxx_fd <= 0 && dxy_fd <= 0 && dyy_fd <= 0 &&
904
pcurv_fd <= 0 && tcurv_fd <= 0)
907
/* compute second order derivatives */
908
s4 = *c1 + *c3 + *c7 + *c9 - *c5 * 8.;
909
s5 = *c4 * 4. + *c6 * 4. - *c8 * 2. - *c2 * 2.;
910
s6 = *c8 * 4. + *c2 * 4. - *c4 * 2. - *c6 * 2.;
911
s3 = *c7 - *c9 + *c3 - *c1;
913
dxx = -(s4 + s5) / ((3. / 32.) * H * H);
914
dyy = -(s4 + s6) / ((3. / 32.) * V * V);
915
dxy = -s3 / ((1. / 16.) * H * V);
918
if (out_type == CELL_TYPE)
919
*((CELL *) dxx_ptr) = (CELL) (scik1 * dxx);
921
G_set_raster_value_d(dxx_ptr, (DCELL) dxx, data_type);
922
dxx_ptr = G_incr_void_ptr(dxx_ptr, G_raster_size(data_type));
926
if (out_type == CELL_TYPE)
927
*((CELL *) dyy_ptr) = (CELL) (scik1 * dyy);
929
G_set_raster_value_d(dyy_ptr, (DCELL) dyy, data_type);
930
dyy_ptr = G_incr_void_ptr(dyy_ptr, G_raster_size(data_type));
934
if (out_type == CELL_TYPE)
935
*((CELL *) dxy_ptr) = (CELL) (scik1 * dxy);
937
G_set_raster_value_d(dxy_ptr, (DCELL) dxy, data_type);
938
dxy_ptr = G_incr_void_ptr(dxy_ptr, G_raster_size(data_type));
941
/* compute curvature */
942
if (pcurv_fd <= 0 && tcurv_fd <= 0)
945
grad2 = key; /*dx2 + dy2 */
947
if (grad <= gradmin) {
952
dnorm1 = sqrt(grad2 + 1.);
953
dxy2 = 2. * dxy * dx * dy;
956
pcurv = (dxx * dx2 + dxy2 + dyy * dy2) /
957
(grad2 * dnorm1 * dnorm1 * dnorm1);
958
tcurv = (dxx * dy2 - dxy2 + dyy * dx2) / (grad2 * dnorm1);
970
if (out_type == CELL_TYPE)
971
*((CELL *) pcurv_ptr) = (CELL) (scik1 * pcurv);
973
G_set_raster_value_d(pcurv_ptr, (DCELL) pcurv, data_type);
975
G_incr_void_ptr(pcurv_ptr, G_raster_size(data_type));
979
if (out_type == CELL_TYPE)
980
*((CELL *) tcurv_ptr) = (CELL) (scik1 * tcurv);
982
G_set_raster_value_d(tcurv_ptr, (DCELL) tcurv, data_type);
984
G_incr_void_ptr(tcurv_ptr, G_raster_size(data_type));
987
} /* column for loop */
990
G_put_raster_row(aspect_fd, asp_raster, data_type);
993
G_put_raster_row(slope_fd, slp_raster, data_type);
996
G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
999
G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
1002
G_put_raster_row(dx_fd, dx_raster, data_type);
1005
G_put_raster_row(dy_fd, dy_raster, data_type);
1008
G_put_raster_row(dxx_fd, dxx_raster, data_type);
1011
G_put_raster_row(dyy_fd, dyy_raster, data_type);
1014
G_put_raster_row(dxy_fd, dxy_raster, data_type);
1018
G_percent(row, nrows, 2);
1020
G_close_cell(elevation_fd);
1021
G_message(_("Creating support files..."));
1023
G_message(_("Elevation products for mapset [%s] in [%s]"),
1024
G_mapset(), G_location());
1026
if (aspect_fd >= 0) {
1028
struct FPRange range;
1030
G_set_null_value(asp_raster, G_window_cols(), data_type);
1031
G_put_raster_row(aspect_fd, asp_raster, data_type);
1032
G_close_cell(aspect_fd);
1034
if (out_type != CELL_TYPE)
1035
G_quantize_fp_map_range(aspect_name, G_mapset(), 0., 360., 0,
1038
G_read_raster_cats(aspect_name, G_mapset(), &cats);
1039
G_set_raster_cats_title
1040
("Aspect counterclockwise in degrees from east", &cats);
1042
G_message(_("Min computed aspect %.4f, max computed aspect %.4f"),
1044
/* the categries quant intervals are 1.0 long, plus
1045
we are using reverse order so that the label looked up
1046
for i-.5 is not the one defined for i-.5, i+.5 interval, but
1047
the one defile for i-1.5, i-.5 interval which is added later */
1048
for (i = ceil(max_asp); i >= 1; i--) {
1050
sprintf(buf, "east");
1052
sprintf(buf, "east");
1054
sprintf(buf, "north ccw of east");
1056
sprintf(buf, "north");
1058
sprintf(buf, "north ccw of west");
1060
sprintf(buf, "west");
1062
sprintf(buf, "south ccw of west");
1064
sprintf(buf, "south");
1066
sprintf(buf, "south ccw of east");
1068
sprintf(buf, "%d degree%s ccw from east", i,
1070
if (data_type == CELL_TYPE) {
1071
G_set_cat(i, buf, &cats);
1074
tmp1 = (double)i - .5;
1075
tmp2 = (double)i + .5;
1076
G_set_d_raster_cat(&tmp1, &tmp2, buf, &cats);
1078
if (data_type == CELL_TYPE)
1079
G_set_cat(0, "no aspect", &cats);
1083
G_set_d_raster_cat(&tmp1, &tmp2, "no aspect", &cats);
1085
G_write_raster_cats(aspect_name, &cats);
1086
G_free_raster_cats(&cats);
1088
/* write colors for aspect file */
1089
G_init_colors(&colors);
1090
G_read_fp_range(aspect_name, G_mapset(), &range);
1091
G_get_fp_range_min_max(&range, &min, &max);
1092
G_make_aspect_fp_colors(&colors, min, max);
1093
G_write_colors(aspect_name, G_mapset(), &colors);
1095
/* writing history file */
1096
G_short_history(aspect_name, "raster", &hist);
1097
sprintf(hist.edhist[0], "aspect map elev = %s", elev_name);
1098
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1099
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1100
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1102
G_write_history(aspect_name, &hist);
1104
G_message(_("ASPECT [%s] COMPLETE"), aspect_name);
1107
if (slope_fd >= 0) {
1108
/* colortable for slopes */
1109
G_init_colors(&colors);
1110
G_add_color_rule(0, 255, 255, 255, 2, 255, 255, 0, &colors);
1111
G_add_color_rule(2, 255, 255, 0, 5, 0, 255, 0, &colors);
1112
G_add_color_rule(5, 0, 255, 0, 10, 0, 255, 255, &colors);
1113
G_add_color_rule(10, 0, 255, 255, 15, 0, 0, 255, &colors);
1114
G_add_color_rule(15, 0, 0, 255, 30, 255, 0, 255, &colors);
1115
G_add_color_rule(30, 255, 0, 255, 50, 255, 0, 0, &colors);
1116
G_add_color_rule(50, 255, 0, 0, 90, 0, 0, 0, &colors);
1118
G_set_null_value(slp_raster, G_window_cols(), data_type);
1119
G_put_raster_row(slope_fd, slp_raster, data_type);
1120
G_close_cell(slope_fd);
1122
if (out_type != CELL_TYPE) {
1125
G_quantize_fp_map_range(slope_name, G_mapset(), 0., 90., 1, 91);
1127
G_quantize_fp_map_range(slope_name, G_mapset(), min_slp, max_slp,
1128
(CELL) min_slp + 1, (CELL) ceil(max_slp) + 1);
1130
G_write_colors(slope_name, G_mapset(), &colors);
1132
G_quantize_fp_map_range(slope_name, G_mapset(), 0., 90., 0,
1135
G_quantize_fp_map_range(slope_name, G_mapset(), min_slp,
1136
max_slp, (CELL) min_slp,
1137
(CELL) ceil(max_slp));
1140
G_read_raster_cats(slope_name, G_mapset(), &cats);
1142
G_set_raster_cats_title("slope in degrees", &cats);
1144
G_set_raster_cats_title("percent slope", &cats);
1146
G_message(_("Min computed slope %.4f, max computed slope %.4f"),
1148
/* the categries quant intervals are 1.0 long, plus
1149
we are using reverse order so that the label looked up
1150
for i-.5 is not the one defined for i-.5, i+.5 interval, but
1151
the one defined for i-1.5, i-.5 interval which is added later */
1152
for (i = ceil(max_slp); i > /* INC BY ONE >= */ 0; i--) {
1154
sprintf(buf, "%d degree%s", i, i == 1 ? "" : "s");
1156
sprintf(buf, "%d percent", i);
1157
if (data_type == CELL_TYPE) {
1159
G_set_cat(i+1, buf, &cats);
1161
G_set_cat(i, buf, &cats);
1165
tmp1 = (DCELL) i+.5;
1166
tmp2 = (DCELL) i+1.5;
1168
tmp1 = (DCELL) i - .5;
1169
tmp2 = (DCELL) i + .5;
1170
G_set_d_raster_cat(&tmp1, &tmp2, buf, &cats);
1172
if (data_type == CELL_TYPE)
1173
G_set_cat(0, "zero slope", &cats);
1175
G_set_cat(0, "no data", &cats);
1180
G_set_d_raster_cat(&tmp1, &tmp2, "zero slope", &cats);
1183
G_set_d_raster_cat (&tmp1, &tmp1, "no data", &cats);
1185
G_write_raster_cats(slope_name, &cats);
1187
/* writing history file */
1188
G_short_history(slope_name, "raster", &hist);
1189
sprintf(hist.edhist[0], "slope map elev = %s", elev_name);
1190
sprintf(hist.edhist[1], "zfactor = %.2f format = %s", zfactor,
1191
parm.slope_fmt->answer);
1192
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1193
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1195
G_write_history(slope_name, &hist);
1197
G_message(_("SLOPE [%s] COMPLETE"), slope_name);
1200
/* colortable for curvatures */
1201
if (pcurv_fd >= 0 || tcurv_fd >= 0) {
1202
G_init_colors(&colors);
1204
dat1 = (FCELL) c1min;
1206
dat1 = (FCELL) c2min;
1208
dat2 = (FCELL) - 0.01;
1209
G_add_f_raster_color_rule(&dat1, 127, 0, 255,
1210
&dat2, 0, 0, 255, &colors);
1212
dat2 = (FCELL) - 0.001;
1213
G_add_f_raster_color_rule(&dat1, 0, 0, 255,
1214
&dat2, 0, 127, 255, &colors);
1216
dat2 = (FCELL) - 0.00001;
1217
G_add_f_raster_color_rule(&dat1, 0, 127, 255,
1218
&dat2, 0, 255, 255, &colors);
1221
G_add_f_raster_color_rule(&dat1, 0, 255, 255,
1222
&dat2, 200, 255, 200, &colors);
1224
dat2 = (FCELL) 0.00001;
1225
G_add_f_raster_color_rule(&dat1, 200, 255, 200,
1226
&dat2, 255, 255, 0, &colors);
1228
dat2 = (FCELL) 0.001;
1229
G_add_f_raster_color_rule(&dat1, 255, 255, 0,
1230
&dat2, 255, 127, 0, &colors);
1232
dat2 = (FCELL) 0.01;
1233
G_add_f_raster_color_rule(&dat1, 255, 127, 0,
1234
&dat2, 255, 0, 0, &colors);
1237
dat2 = (FCELL) c1max;
1239
dat2 = (FCELL) c2max;
1241
G_add_f_raster_color_rule(&dat1, 255, 0, 0,
1242
&dat2, 255, 0, 200, &colors);
1245
if (pcurv_fd >= 0) {
1246
G_set_null_value(pcurv_raster, G_window_cols(), data_type);
1247
G_put_raster_row(pcurv_fd, pcurv_raster, data_type);
1248
G_close_cell(pcurv_fd);
1250
G_write_colors(pcurv_name, G_mapset(), &colors);
1252
if (out_type != CELL_TYPE)
1253
G_round_fp_map(pcurv_name, G_mapset());
1255
G_read_cats(pcurv_name, G_mapset(), &cats);
1256
G_set_cats_title("profile curvature", &cats);
1257
G_set_cat((CELL) 0, "no profile curve", &cats);
1259
/* writing history file */
1260
G_short_history(pcurv_name, "raster", &hist);
1261
sprintf(hist.edhist[0], "profile curve map elev = %s", elev_name);
1262
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1263
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1264
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1266
G_write_history(pcurv_name, &hist);
1268
G_message(_("PROFILE CURVE [%s] COMPLETE"), pcurv_name);
1271
if (tcurv_fd >= 0) {
1272
G_set_null_value(tcurv_raster, G_window_cols(), data_type);
1273
G_put_raster_row(tcurv_fd, tcurv_raster, data_type);
1274
G_close_cell(tcurv_fd);
1276
G_write_colors(tcurv_name, G_mapset(), &colors);
1278
if (out_type != CELL_TYPE)
1279
G_round_fp_map(tcurv_name, G_mapset());
1281
G_read_cats(tcurv_name, G_mapset(), &cats);
1282
G_set_cats_title("tangential curvature", &cats);
1283
G_set_cat((CELL) 0, "no tangential curve", &cats);
1285
/* writing history file */
1286
G_short_history(tcurv_name, "raster", &hist);
1287
sprintf(hist.edhist[0], "tangential curve map elev = %s", elev_name);
1288
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1289
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1290
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1292
G_write_history(tcurv_name, &hist);
1294
G_message(_("TANGENTIAL CURVE [%s] COMPLETE"), tcurv_name);
1298
G_set_null_value(dx_raster, G_window_cols(), data_type);
1299
G_put_raster_row(dx_fd, dx_raster, data_type);
1300
G_close_cell(dx_fd);
1302
if (out_type != CELL_TYPE)
1303
G_round_fp_map(dx_name, G_mapset());
1305
G_read_cats(dx_name, G_mapset(), &cats);
1306
G_set_cats_title("E-W slope", &cats);
1307
G_set_cat((CELL) 0, "no E-W slope", &cats);
1309
/* writing history file */
1310
G_short_history(dx_name, "raster", &hist);
1311
sprintf(hist.edhist[0], "E-W slope map elev = %s", elev_name);
1312
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1313
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1314
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1316
G_write_history(dx_name, &hist);
1318
G_message(_("E-W SLOPE [%s] COMPLETE"), dx_name);
1322
G_set_null_value(dy_raster, G_window_cols(), data_type);
1323
G_put_raster_row(dy_fd, dy_raster, data_type);
1324
G_close_cell(dy_fd);
1326
if (out_type != CELL_TYPE)
1327
G_round_fp_map(dy_name, G_mapset());
1329
G_read_cats(dy_name, G_mapset(), &cats);
1330
G_set_cats_title("N-S slope", &cats);
1331
G_set_cat((CELL) 0, "no N-S slope", &cats);
1333
/* writing history file */
1334
G_short_history(dy_name, "raster", &hist);
1335
sprintf(hist.edhist[0], "N-S slope map elev = %s", elev_name);
1336
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1337
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1338
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1340
G_write_history(dy_name, &hist);
1342
G_message(_("N-S SLOPE [%s] COMPLETE"), dy_name);
1346
G_set_null_value(dxx_raster, G_window_cols(), data_type);
1347
G_put_raster_row(dxx_fd, dxx_raster, data_type);
1348
G_close_cell(dxx_fd);
1350
if (out_type != CELL_TYPE)
1351
G_round_fp_map(dxx_name, G_mapset());
1353
G_read_cats(dxx_name, G_mapset(), &cats);
1354
G_set_cats_title("DXX", &cats);
1355
G_set_cat((CELL) 0, "DXX", &cats);
1357
/* writing history file */
1358
G_short_history(dxx_name, "raster", &hist);
1359
sprintf(hist.edhist[0], "DXX map elev = %s", elev_name);
1360
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1361
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1362
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1364
G_write_history(dxx_name, &hist);
1366
G_message(_("DXX [%s] COMPLETE"), dxx_name);
1370
G_set_null_value(dyy_raster, G_window_cols(), data_type);
1371
G_put_raster_row(dyy_fd, dyy_raster, data_type);
1372
G_close_cell(dyy_fd);
1374
if (out_type != CELL_TYPE)
1375
G_round_fp_map(dyy_name, G_mapset());
1377
G_read_cats(dyy_name, G_mapset(), &cats);
1378
G_set_cats_title("DYY", &cats);
1379
G_set_cat((CELL) 0, "DYY", &cats);
1381
/* writing history file */
1382
G_short_history(dyy_name, "raster", &hist);
1383
sprintf(hist.edhist[0], "DYY map elev = %s", elev_name);
1384
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1385
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1386
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1388
G_write_history(dyy_name, &hist);
1390
G_message(_("DYY [%s] COMPLETE"), dyy_name);
1394
G_set_null_value(dxy_raster, G_window_cols(), data_type);
1395
G_put_raster_row(dxy_fd, dxy_raster, data_type);
1396
G_close_cell(dxy_fd);
1398
if (out_type != CELL_TYPE)
1399
G_round_fp_map(dxy_name, G_mapset());
1401
G_read_cats(dxy_name, G_mapset(), &cats);
1402
G_set_cats_title("DXY", &cats);
1403
G_set_cat((CELL) 0, "DXY", &cats);
1405
/* writing history file */
1406
G_short_history(dxy_name, "raster", &hist);
1407
sprintf(hist.edhist[0], "DXY map elev = %s", elev_name);
1408
sprintf(hist.edhist[1], "zfactor = %.2f", zfactor);
1409
sprintf(hist.edhist[2], "min_slp_allowed = %f", min_slp_allowed);
1410
sprintf(hist.datsrc_1, "raster elevation file %s", elev_name);
1412
G_write_history(dxy_name, &hist);
1414
G_message(_("DXY [%s] COMPLETE"), dxy_name);
added
removed
swig/perl/R_slope_aspect/r_slope_aspect/r_slope_aspect.c
