1
2 /****************************************************************************
3 *
4 * MODULE: d.profile
5 * AUTHOR(S): Dave Johnson (original contributor)
6 * DBA Systems, Inc. 10560 Arrowhead Drive Fairfax, VA 22030
7 * Markus Neteler <neteler itc.it>,
8 * Bernhard Reiter <bernhard intevation.de>,
9 * Huidae Cho <grass4u gmail.com>,
10 * Eric G. Miller <egm2 jps.net>,
11 * Glynn Clements <glynn gclements.plus.com>
12 * PURPOSE: user chooses transects path, and profile of raster data drawn
13 * COPYRIGHT: (C) 1999-2007 by the GRASS Development Team
14 *
15 * This program is free software under the GNU General Public
16 * License (>=v2). Read the file COPYING that comes with GRASS
17 * for details.
18 *
19 *****************************************************************************/
20
21
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <math.h>
25 #include <grass/gis.h>
26 #include <grass/raster.h>
27 #include <grass/display.h>
28 #include <grass/glocale.h>
29
30 static char *mapname;
31 static double min, max;
32
33 struct point
34 {
35 double x, y;
36 double d;
37 };
38
get_region_range(int fd)39 static void get_region_range(int fd)
40 {
41 DCELL *buf = Rast_allocate_d_buf();
42
43 int nrows = Rast_window_rows();
44 int ncols = Rast_window_cols();
45 int row, col;
46
47 min = 1e300;
48 max = -1e300;
49
50 for (row = 0; row < nrows; row++) {
51 Rast_get_d_row(fd, buf, row);
52 for (col = 0; col < ncols; col++) {
53 if (min > buf[col])
54 min = buf[col];
55 if (max < buf[col])
56 max = buf[col];
57 }
58 }
59 }
60
get_map_range(void)61 static void get_map_range(void)
62 {
63 if (Rast_map_type(mapname, "") == CELL_TYPE) {
64 struct Range range;
65 CELL xmin, xmax;
66
67 if (Rast_read_range(mapname, "", &range) <= 0)
68 G_fatal_error(_("Unable to read range for %s"), mapname);
69
70 Rast_get_range_min_max(&range, &xmin, &xmax);
71
72 max = xmax;
73 min = xmin;
74 }
75 else {
76 struct FPRange fprange;
77
78 if (Rast_read_fp_range(mapname, "", &fprange) <= 0)
79 G_fatal_error(_("Unable to read FP range for %s"), mapname);
80
81 Rast_get_fp_range_min_max(&fprange, &min, &max);
82 }
83 }
84
plot_axes(void)85 static void plot_axes(void)
86 {
87 char str[64];
88 double scale;
89 double t, b, l, r;
90
91 D_use_color(D_translate_color("red"));
92
93 D_begin();
94 D_move_abs(0, 1);
95 D_cont_abs(0, 0);
96 D_cont_abs(1, 0);
97 D_end();
98 D_stroke();
99
100 D_use_color(D_translate_color(DEFAULT_FG_COLOR));
101
102 /* set text size for y-axis labels */
103 scale = fabs(D_get_u_to_d_yconv());
104 D_text_size(scale * 0.04, scale * 0.05);
105
106 /* plot y-axis label (bottom) */
107 sprintf(str, "%.1f", min);
108 D_get_text_box(str, &t, &b, &l, &r);
109 D_pos_abs(-0.02 - (r - l), 0 - (t - b) / 2);
110 D_text(str);
111
112 /* plot y-axis label (top) */
113 sprintf(str, "%.1f", max);
114 D_get_text_box(str, &t, &b, &l, &r);
115 D_pos_abs(-0.02 - (r - l), 1 - (t - b) / 2);
116 D_text(str);
117 }
118
get_cell(DCELL * result,int fd,double x,double y)119 static int get_cell(DCELL *result, int fd, double x, double y)
120 {
121 static DCELL *row1, *row2;
122 static int cur_row = -1;
123 static int row, col;
124 DCELL *tmp;
125
126 if (!row1) {
127 row1 = Rast_allocate_d_buf();
128 row2 = Rast_allocate_d_buf();
129 }
130
131 col = (int)floor(x - 0.5);
132 row = (int)floor(y - 0.5);
133 x -= col + 0.5;
134 y -= row + 0.5;
135
136 if (row < 0 || row + 1 >= Rast_window_rows() ||
137 col < 0 || col + 1 >= Rast_window_cols()) {
138 Rast_set_d_null_value(result, 1);
139 return 0;
140 }
141
142 if (cur_row != row) {
143 if (cur_row == row + 1) {
144 tmp = row1; row1 = row2; row2 = tmp;
145 Rast_get_d_row(fd, row1, row);
146 }
147 else if (cur_row == row - 1) {
148 tmp = row1; row1 = row2; row2 = tmp;
149 Rast_get_d_row(fd, row2, row + 1);
150 }
151 else {
152 Rast_get_d_row(fd, row1, row);
153 Rast_get_d_row(fd, row2, row + 1);
154 }
155 cur_row = row;
156 }
157
158 if (Rast_is_d_null_value(&row1[col]) || Rast_is_d_null_value(&row1[col+1]) ||
159 Rast_is_d_null_value(&row2[col]) || Rast_is_d_null_value(&row2[col+1])) {
160 Rast_set_d_null_value(result, 1);
161 return 0;
162 }
163
164 *result = Rast_interp_bilinear(x, y,
165 row1[col], row1[col+1],
166 row2[col], row2[col+1]);
167
168 return 1;
169 }
170
main(int argc,char ** argv)171 int main(int argc, char **argv)
172 {
173 struct GModule *module;
174 struct Option *map, *profile;
175 struct Flag *stored;
176 struct Cell_head window;
177 struct point *points = NULL;
178 int num_points, max_points = 0;
179 double length;
180 double t, b, l, r;
181 int fd;
182 int i;
183 double sx;
184 int last;
185
186 /* Initialize the GIS calls */
187 G_gisinit(argv[0]);
188
189 /* Set description */
190 module = G_define_module();
191 G_add_keyword(_("display"));
192 G_add_keyword(_("profile"));
193 G_add_keyword(_("raster"));
194 G_add_keyword(_("plot"));
195 module->description = _("Plots profile of a transect.");
196
197 /* set up command line */
198 map = G_define_standard_option(G_OPT_R_MAP);
199 map->description = _("Raster map to be profiled");
200
201 profile = G_define_standard_option(G_OPT_M_COORDS);
202 profile->required = YES;
203 profile->multiple = YES;
204 profile->description = _("Profile coordinate pairs");
205
206 stored = G_define_flag();
207 stored->key = 'r';
208 stored->description = _("Use map's range recorded range");
209
210 if (G_parser(argc, argv))
211 exit(EXIT_FAILURE);
212
213 mapname = map->answer;
214
215 fd = Rast_open_old(mapname, "");
216
217 if (stored->answer)
218 get_map_range();
219 else
220 get_region_range(fd);
221
222 G_get_window(&window);
223
224 num_points = 0;
225 length = 0;
226 for (i = 0; profile->answers[i]; i += 2) {
227 struct point *p;
228 double x, y;
229
230 if (num_points >= max_points) {
231 max_points = num_points + 100;
232 points = G_realloc(points, max_points * sizeof(struct point));
233 }
234
235 p = &points[num_points];
236
237 G_scan_easting( profile->answers[i+0], &x, G_projection());
238 G_scan_northing(profile->answers[i+1], &y, G_projection());
239
240 p->x = Rast_easting_to_col (x, &window);
241 p->y = Rast_northing_to_row(y, &window);
242
243 if (num_points > 0) {
244 const struct point *prev = &points[num_points-1];
245 double dx = fabs(p->x - prev->x);
246 double dy = fabs(p->y - prev->y);
247 double d = sqrt(dx * dx + dy * dy);
248 length += d;
249 p->d = length;
250 }
251
252 num_points++;
253 }
254 points[0].d = 0;
255
256 if (num_points < 2)
257 G_fatal_error(_("At least two points are required"));
258
259 /* establish connection with graphics driver */
260 D_open_driver();
261
262 D_setup2(1, 0, 1.05, -0.05, -0.15, 1.05);
263
264 plot_axes();
265
266 D_use_color(D_translate_color(DEFAULT_FG_COLOR));
267
268 D_get_src(&t, &b, &l, &r);
269 t -= 0.1 * (t - b);
270 b += 0.1 * (t - b);
271 l += 0.1 * (r - l);
272 r -= 0.1 * (r - l);
273
274 D_begin();
275
276 i = 0;
277 last = 0;
278 for (sx = 0; sx < 1; sx += D_get_d_to_u_xconv()) {
279 double d = length * (sx - l);
280 const struct point *p, *next;
281 double k, sy, x, y;
282 DCELL v;
283
284 for (;;) {
285 p = &points[i];
286 next = &points[i + 1];
287 k = (d - p->d) / (next->d - p->d);
288 if (k < 1)
289 break;
290 i++;
291 }
292
293 x = p->x * (1 - k) + next->x * k;
294 y = p->y * (1 - k) + next->y * k;
295
296 if (!get_cell(&v, fd, x, y)) {
297 last = 0;
298 continue;
299 }
300
301 sy = (v - min) / (max - min);
302
303 if (last)
304 D_cont_abs(sx, sy);
305 else
306 D_move_abs(sx, sy);
307
308 last = 1;
309 }
310
311 D_end();
312 D_stroke();
313
314 D_close_driver();
315
316 exit(EXIT_SUCCESS);
317 }
318
319