1 #include <string.h>
2 #include <grass/gis.h>
3 #include <grass/raster.h>
4 #include <grass/glocale.h>
5 #include "global.h"
6
read_cells(void)7 void read_cells(void)
8 {
9 int fd, i, j;
10 RASTER_MAP_TYPE data_type;
11 CELL *ccell = NULL;
12 FCELL *fcell = NULL;
13 struct Cell_head inhead;
14 char buf_wrns[32], buf_wrew[32], buf_mrns[32], buf_mrew[32];
15
16 fd = Rast_open_old(input, "");
17
18 data_type = Rast_get_map_type(fd);
19 Rast_get_cellhd(input, "", &inhead);
20
21 if (data_type == CELL_TYPE)
22 ccell = (CELL *) G_malloc(sizeof(CELL) * window.cols);
23 else if (data_type == FCELL_TYPE)
24 fcell = (FCELL *) G_malloc(sizeof(FCELL) * window.cols);
25
26 cell = (DCELL **) G_malloc(sizeof(DCELL *) * window.rows);
27 atb = (DCELL **) G_malloc(sizeof(DCELL *) * window.rows);
28 a = (DCELL **) G_malloc(sizeof(DCELL *) * window.rows);
29
30 if (window.ew_res < inhead.ew_res || window.ns_res < inhead.ns_res) {
31 G_format_resolution(window.ew_res, buf_wrew, G_projection());
32 G_format_resolution(window.ns_res, buf_wrns, G_projection());
33 G_format_resolution(inhead.ew_res, buf_mrew, G_projection());
34 G_format_resolution(inhead.ns_res, buf_mrns, G_projection());
35 G_fatal_error(_("The current region resolution [%s x %s] is finer "
36 "than the input map's resolution [%s x %s]. "
37 "The current region resolution must be identical "
38 "to, or coarser than, the input map's resolution."),
39 buf_wrew, buf_wrns, buf_mrew, buf_mrns);
40 }
41
42 G_message(_("Reading elevation map..."));
43
44 for (i = 0; i < window.rows; i++) {
45 G_percent(i, window.rows, 2);
46
47 cell[i] = (DCELL *) G_malloc(sizeof(DCELL) * window.cols);
48 atb[i] = (DCELL *) G_malloc(sizeof(DCELL) * window.cols);
49 a[i] = (DCELL *) G_malloc(sizeof(DCELL) * window.cols);
50
51 if (data_type == CELL_TYPE) {
52 Rast_get_c_row(fd, ccell, i);
53 for (j = 0; j < window.cols; j++) {
54 if (Rast_is_c_null_value(&ccell[j]))
55 Rast_set_d_null_value(&cell[i][j], 1);
56 else
57 cell[i][j] = (DCELL) ccell[j];
58 }
59 }
60 else if (data_type == FCELL_TYPE) {
61 Rast_get_f_row(fd, fcell, i);
62 for (j = 0; j < window.cols; j++) {
63 if (Rast_is_f_null_value(&fcell[j]))
64 Rast_set_d_null_value(&cell[i][j], 1);
65 else
66 cell[i][j] = (DCELL) fcell[j];
67 }
68 }
69 else
70 Rast_get_d_row(fd, cell[i], i);
71 }
72 if (data_type == CELL_TYPE)
73 G_free(ccell);
74 else if (data_type == FCELL_TYPE)
75 G_free(fcell);
76 G_percent(i, window.rows, 2);
77 Rast_close(fd);
78 }
79
write_cells(void)80 void write_cells(void)
81 {
82 int fd, i;
83 struct History history;
84
85 fd = Rast_open_new(output, DCELL_TYPE);
86
87 G_message(_("Writing topographic index map..."));
88
89 for (i = 0; i < window.rows; i++) {
90 G_percent(i, window.rows, 2);
91 Rast_put_d_row(fd, atb[i]);
92 }
93 G_percent(i, window.rows, 2);
94 Rast_close(fd);
95
96 Rast_short_history(output, "raster", &history);
97 Rast_set_history(&history, HIST_DATSRC_1, input);
98 Rast_write_history(output, &history);
99 }
100