GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
resout2d.c
Go to the documentation of this file.
1
2/*-
3 * Written by H. Mitasova, I. Kosinovsky, D. Gerdes Summer 1993
4 * University of Illinois
5 * US Army Construction Engineering Research Lab
6 * Copyright 1993, H. Mitasova (University of Illinois),
7 * I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
8 *
9 * modified by McCauley in August 1995
10 * modified by Mitasova in August 1995
11 *
12 */
13
14#define MULT 100000
15
16#include <stdio.h>
17#include <math.h>
18
19#include <grass/gis.h>
20#include <grass/raster.h>
21#include <grass/bitmap.h>
22#include <grass/linkm.h>
23#include <grass/interpf.h>
24#include <grass/glocale.h>
25
26/* output cell maps for elevation, aspect, slope and curvatures */
27
28static void do_history(const char *name, const char *input,
29 const struct interp_params *params)
30{
31 struct History hist;
32
33 Rast_short_history(name, "raster", &hist);
34 if (params->elev)
35 Rast_append_format_history(&hist, "The elevation map is %s",
36 params->elev);
37
38 Rast_format_history(&hist, HIST_DATSRC_1, "raster map %s", input);
39
40 Rast_write_history(name, &hist);
41
42 Rast_free_history(&hist);
43}
44
45int IL_resample_output_2d(struct interp_params *params, double zmin, double zmax, /* min,max input z-values */
46 double zminac, double zmaxac, /* min,max interpolated values */
47 double c1min, double c1max, double c2min, double c2max, double gmin, double gmax, double ertot, /* total interplating func. error */
48 char *input, /* input file name */
49 double *dnorm, struct Cell_head *outhd, /* Region with desired resolution */
50 struct Cell_head *winhd, /* Current region */
51 char *smooth, int n_points)
52
53/*
54 * Creates output files as well as history files and color tables for
55 * them.
56 */
57{
58 FCELL *cell1; /* cell buffer */
59 int cf1 = 0, cf2 = 0, cf3 = 0, cf4 = 0, cf5 = 0, cf6 = 0; /* cell file descriptors */
60 int nrows, ncols; /* current region rows and columns */
61 int i; /* loop counter */
62 const char *mapset;
63 float dat1, dat2;
64 struct Colors colors, colors2;
65 double value1, value2;
66 struct History hist;
67 struct _Color_Rule_ *rule;
68 const char *maps;
69 int cond1, cond2;
70 CELL val1, val2;
71
72 cond2 = ((params->pcurv != NULL) ||
73 (params->tcurv != NULL) || (params->mcurv != NULL));
74 cond1 = ((params->slope != NULL) || (params->aspect != NULL) || cond2);
75
76 /* change region to output cell file region */
77 G_verbose_message(_("Temporarily changing the region to desired resolution..."));
78 Rast_set_output_window(outhd);
79 mapset = G_mapset();
80
81 cell1 = Rast_allocate_f_output_buf();
82
83 if (params->elev)
84 cf1 = Rast_open_fp_new(params->elev);
85
86 if (params->slope)
87 cf2 = Rast_open_fp_new(params->slope);
88
89 if (params->aspect)
90 cf3 = Rast_open_fp_new(params->aspect);
91
92 if (params->pcurv)
93 cf4 = Rast_open_fp_new(params->pcurv);
94
95 if (params->tcurv)
96 cf5 = Rast_open_fp_new(params->tcurv);
97
98 if (params->mcurv)
99 cf6 = Rast_open_fp_new(params->mcurv);
100
101 nrows = outhd->rows;
102 if (nrows != params->nsizr) {
103 G_warning(_("First change your rows number(%d) to %d"),
104 nrows, params->nsizr);
105 return -1;
106 }
107
108 ncols = outhd->cols;
109 if (ncols != params->nsizc) {
110 G_warning(_("First change your columns number(%d) to %d"),
111 ncols, params->nsizr);
112 return -1;
113 }
114
115 if (params->elev != NULL) {
116 G_fseek(params->Tmp_fd_z, 0L, 0); /* seek to the beginning */
117 for (i = 0; i < params->nsizr; i++) {
118 /* seek to the right row */
119 G_fseek(params->Tmp_fd_z, (off_t) (params->nsizr - 1 - i) *
120 params->nsizc * sizeof(FCELL), 0);
121 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_z);
122 Rast_put_f_row(cf1, cell1);
123 }
124 }
125
126 if (params->slope != NULL) {
127 G_fseek(params->Tmp_fd_dx, 0L, 0); /* seek to the beginning */
128 for (i = 0; i < params->nsizr; i++) {
129 /* seek to the right row */
130 G_fseek(params->Tmp_fd_dx, (off_t) (params->nsizr - 1 - i) *
131 params->nsizc * sizeof(FCELL), 0);
132 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dx);
133 /*
134 * for (ii==0;ii<params->nsizc;ii++) { fprintf(stderr,"ii=%d ",ii);
135 * fprintf(stderr,"%f ",cell1[ii]); }
136 * fprintf(stderr,"params->nsizc=%d \n",params->nsizc);
137 */
138 Rast_put_f_row(cf2, cell1);
139 }
140 }
141
142 if (params->aspect != NULL) {
143 G_fseek(params->Tmp_fd_dy, 0L, 0); /* seek to the beginning */
144 for (i = 0; i < params->nsizr; i++) {
145 /* seek to the right row */
146 G_fseek(params->Tmp_fd_dy, (off_t) (params->nsizr - 1 - i) *
147 params->nsizc * sizeof(FCELL), 0);
148 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_dy);
149 Rast_put_f_row(cf3, cell1);
150 }
151 }
152
153 if (params->pcurv != NULL) {
154 G_fseek(params->Tmp_fd_xx, 0L, 0); /* seek to the beginning */
155 for (i = 0; i < params->nsizr; i++) {
156 /* seek to the right row */
157 G_fseek(params->Tmp_fd_xx, (off_t) (params->nsizr - 1 - i) *
158 params->nsizc * sizeof(FCELL), 0);
159 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xx);
160 Rast_put_f_row(cf4, cell1);
161 }
162 }
163
164 if (params->tcurv != NULL) {
165 G_fseek(params->Tmp_fd_yy, 0L, 0); /* seek to the beginning */
166 for (i = 0; i < params->nsizr; i++) {
167 /* seek to the right row */
168 G_fseek(params->Tmp_fd_yy, (off_t) (params->nsizr - 1 - i) *
169 params->nsizc * sizeof(FCELL), 0);
170 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_yy);
171 Rast_put_f_row(cf5, cell1);
172 }
173 }
174
175 if (params->mcurv != NULL) {
176 G_fseek(params->Tmp_fd_xy, 0L, 0); /* seek to the beginning */
177 for (i = 0; i < params->nsizr; i++) {
178 /* seek to the right row */
179 G_fseek(params->Tmp_fd_xy, (off_t) (params->nsizr - 1 - i) *
180 params->nsizc * sizeof(FCELL), 0);
181 fread(cell1, sizeof(FCELL), params->nsizc, params->Tmp_fd_xy);
182 Rast_put_f_row(cf6, cell1);
183 }
184 }
185
186 if (cf1)
187 Rast_close(cf1);
188 if (cf2)
189 Rast_close(cf2);
190 if (cf3)
191 Rast_close(cf3);
192 if (cf4)
193 Rast_close(cf4);
194 if (cf5)
195 Rast_close(cf5);
196 if (cf6)
197 Rast_close(cf6);
198
199 /* write colormaps and history for output cell files */
200 /* colortable for elevations */
201 maps = G_find_file("cell", input, "");
202
203 if (params->elev != NULL) {
204 if (maps == NULL) {
205 G_warning(_("Raster map <%s> not found"), input);
206 return -1;
207 }
208 Rast_init_colors(&colors2);
209 /*
210 * Rast_mark_colors_as_fp(&colors2);
211 */
212
213 if (Rast_read_colors(input, maps, &colors) >= 0) {
214 if (colors.modular.rules) {
215 rule = colors.modular.rules;
216
217 while (rule->next)
218 rule = rule->next;
219
220 for (; rule; rule = rule->prev) {
221 value1 = rule->low.value * params->zmult;
222 value2 = rule->high.value * params->zmult;
223 Rast_add_modular_d_color_rule(&value1, rule->low.red,
224 rule->low.grn,
225 rule->low.blu, &value2,
226 rule->high.red,
227 rule->high.grn,
228 rule->high.blu,
229 &colors2);
230 }
231 }
232
233 if (colors.fixed.rules) {
234 rule = colors.fixed.rules;
235
236 while (rule->next)
237 rule = rule->next;
238
239 for (; rule; rule = rule->prev) {
240 value1 = rule->low.value * params->zmult;
241 value2 = rule->high.value * params->zmult;
242 Rast_add_d_color_rule(&value1, rule->low.red,
243 rule->low.grn, rule->low.blu,
244 &value2, rule->high.red,
245 rule->high.grn, rule->high.blu,
246 &colors2);
247 }
248 }
249
250 maps = NULL;
251 maps = G_find_file("cell", params->elev, "");
252 if (maps == NULL) {
253 G_warning(_("Raster map <%s> not found"), params->elev);
254 return -1;
255 }
256
257 Rast_write_colors(params->elev, maps, &colors2);
258 Rast_quantize_fp_map_range(params->elev, mapset,
259 zminac - 0.5, zmaxac + 0.5,
260 (CELL) (zminac - 0.5),
261 (CELL) (zmaxac + 0.5));
262 }
263 else
264 G_warning(_("No color table for input raster map -- will not create color table"));
265 }
266
267 /* colortable for slopes */
268 if (cond1 & (!params->deriv)) {
269 Rast_init_colors(&colors);
270 val1 = 0;
271 val2 = 2;
272 Rast_add_c_color_rule(&val1, 255, 255, 255, &val2, 255, 255, 0, &colors);
273 val1 = 2;
274 val2 = 5;
275 Rast_add_c_color_rule(&val1, 255, 255, 0, &val2, 0, 255, 0, &colors);
276 val1 = 5;
277 val2 = 10;
278 Rast_add_c_color_rule(&val1, 0, 255, 0, &val2, 0, 255, 255, &colors);
279 val1 = 10;
280 val2 = 15;
281 Rast_add_c_color_rule(&val1, 0, 255, 255, &val2, 0, 0, 255, &colors);
282 val1 = 15;
283 val2 = 30;
284 Rast_add_c_color_rule(&val1, 0, 0, 255, &val2, 255, 0, 255, &colors);
285 val1 = 30;
286 val2 = 50;
287 Rast_add_c_color_rule(&val1, 255, 0, 255, &val2, 255, 0, 0, &colors);
288 val1 = 50;
289 val2 = 90;
290 Rast_add_c_color_rule(&val1, 255, 0, 0, &val2, 0, 0, 0, &colors);
291
292 if (params->slope != NULL) {
293 maps = NULL;
294 maps = G_find_file("cell", params->slope, "");
295 if (maps == NULL) {
296 G_warning(_("Raster map <%s> not found"), params->slope);
297 return -1;
298 }
299 Rast_write_colors(params->slope, maps, &colors);
300 Rast_quantize_fp_map_range(params->slope, mapset, 0., 90., 0, 90);
301
302 do_history(params->slope, input, params);
303 }
304
305 /* colortable for aspect */
306 Rast_init_colors(&colors);
307 val1 = 0;
308 val2 = 0;
309 Rast_add_c_color_rule(&val1, 255, 255, 255, &val2, 255, 255, 255, &colors);
310 val1 = 1;
311 val2 = 90;
312 Rast_add_c_color_rule(&val1, 255, 255, 0, &val2, 0, 255, 0, &colors);
313 val1 = 90;
314 val2 = 180;
315 Rast_add_c_color_rule(&val1, 0, 255, 0, &val2, 0, 255, 255, &colors);
316 val1 = 180;
317 val2 = 270;
318 Rast_add_c_color_rule(&val1, 0, 255, 255, &val2, 255, 0, 0, &colors);
319 val1 = 270;
320 val2 = 360;
321 Rast_add_c_color_rule(&val1, 255, 0, 0, &val2, 255, 255, 0, &colors);
322
323 if (params->aspect != NULL) {
324 maps = NULL;
325 maps = G_find_file("cell", params->aspect, "");
326 if (maps == NULL) {
327 G_warning(_("Raster map <%s> not found"), params->aspect);
328 return -1;
329 }
330 Rast_write_colors(params->aspect, maps, &colors);
331 Rast_quantize_fp_map_range(params->aspect, mapset, 0., 360., 0, 360);
332
333 do_history(params->aspect, input, params);
334 }
335
336 /* colortable for curvatures */
337 if (cond2) {
338 Rast_init_colors(&colors);
339
340 dat1 = (FCELL) amin1(c1min, c2min);
341 dat2 = (FCELL) - 0.01;
342
343 Rast_add_f_color_rule(&dat1, 50, 0, 155,
344 &dat2, 0, 0, 255, &colors);
345 dat1 = dat2;
346 dat2 = (FCELL) - 0.001;
347 Rast_add_f_color_rule(&dat1, 0, 0, 255,
348 &dat2, 0, 127, 255, &colors);
349 dat1 = dat2;
350 dat2 = (FCELL) - 0.00001;
351 Rast_add_f_color_rule(&dat1, 0, 127, 255,
352 &dat2, 0, 255, 255, &colors);
353 dat1 = dat2;
354 dat2 = (FCELL) 0.00;
355 Rast_add_f_color_rule(&dat1, 0, 255, 255,
356 &dat2, 200, 255, 200, &colors);
357 dat1 = dat2;
358 dat2 = (FCELL) 0.00001;
359 Rast_add_f_color_rule(&dat1, 200, 255, 200,
360 &dat2, 255, 255, 0, &colors);
361 dat1 = dat2;
362 dat2 = (FCELL) 0.001;
363 Rast_add_f_color_rule(&dat1, 255, 255, 0,
364 &dat2, 255, 127, 0, &colors);
365 dat1 = dat2;
366 dat2 = (FCELL) 0.01;
367 Rast_add_f_color_rule(&dat1, 255, 127, 0,
368 &dat2, 255, 0, 0, &colors);
369 dat1 = dat2;
370 dat2 = (FCELL) amax1(c1max, c2max);
371 Rast_add_f_color_rule(&dat1, 255, 0, 0,
372 &dat2, 155, 0, 20, &colors);
373 maps = NULL;
374 if (params->pcurv != NULL) {
375 maps = G_find_file("cell", params->pcurv, "");
376 if (maps == NULL) {
377 G_warning(_("Raster map <%s> not found"), params->pcurv);
378 return -1;
379 }
380 Rast_write_colors(params->pcurv, maps, &colors);
381
382 fprintf(stderr, "color map written\n");
383
384 Rast_quantize_fp_map_range(params->pcurv, mapset,
385 dat1, dat2,
386 (CELL) (dat1 * MULT),
387 (CELL) (dat2 * MULT));
388 do_history(params->pcurv, input, params);
389 }
390
391 if (params->tcurv != NULL) {
392 maps = NULL;
393 maps = G_find_file("cell", params->tcurv, "");
394 if (maps == NULL) {
395 G_warning(_("Raster map <%s> not found"), params->tcurv);
396 return -1;
397 }
398 Rast_write_colors(params->tcurv, maps, &colors);
399 Rast_quantize_fp_map_range(params->tcurv, mapset,
400 dat1, dat2, (CELL) (dat1 * MULT),
401 (CELL) (dat2 * MULT));
402
403 do_history(params->tcurv, input, params);
404 }
405
406 if (params->mcurv != NULL) {
407 maps = NULL;
408 maps = G_find_file("cell", params->mcurv, "");
409 if (maps == NULL) {
410 G_warning(_("Raster map <%s> not found"), params->mcurv);
411 return -1;
412 }
413 Rast_write_colors(params->mcurv, maps, &colors);
414 Rast_quantize_fp_map_range(params->mcurv, mapset,
415 dat1, dat2,
416 (CELL) (dat1 * MULT),
417 (CELL) (dat2 * MULT));
418
419 do_history(params->mcurv, input, params);
420 }
421 }
422 }
423
424 if (params->elev != NULL) {
425 if (!G_find_file2("cell", params->elev, "")) {
426 G_warning(_("Raster map <%s> not found"), params->elev);
427 return -1;
428 }
429
430 Rast_short_history(params->elev, "raster", &hist);
431
432 if (smooth != NULL)
433 Rast_append_format_history(
434 &hist, "tension=%f, smoothing=%s",
435 params->fi * 1000. / (*dnorm), smooth);
436 else
437 Rast_append_format_history(
438 &hist, "tension=%f", params->fi * 1000. / (*dnorm));
439
440 Rast_append_format_history(
441 &hist, "dnorm=%f, zmult=%f", *dnorm, params->zmult);
442 Rast_append_format_history(
443 &hist, "KMAX=%d, KMIN=%d, errtotal=%f", params->kmax,
444 params->kmin, sqrt(ertot / n_points));
445 Rast_append_format_history(
446 &hist, "zmin_data=%f, zmax_data=%f", zmin, zmax);
447 Rast_append_format_history(
448 &hist, "zmin_int=%f, zmax_int=%f", zminac, zmaxac);
449
450 Rast_format_history(&hist, HIST_DATSRC_1, "raster map %s", input);
451
452 Rast_write_history(params->elev, &hist);
453
454 Rast_free_history(&hist);
455 }
456
457/* change region to initial region */
458 G_verbose_message(_("Changing the region back to initial..."));
459 Rast_set_output_window(winhd);
460
461 return 1;
462}
#define NULL
Definition: ccmath.h:32
const char * G_find_file2(const char *element, const char *name, const char *mapset)
Searches for a file from the mapset search list or in a specified mapset. (look but don't touch)
Definition: find_file.c:249
const char * G_find_file(const char *element, char *name, const char *mapset)
Searches for a file from the mapset search list or in a specified mapset.
Definition: find_file.c:203
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition: gis/error.c:109
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
void G_fseek(FILE *fp, off_t offset, int whence)
Change the file position of the stream.
Definition: gis/seek.c:50
double amin1(double, double)
Definition: minmax.c:67
double amax1(double, double)
Definition: minmax.c:54
const char * G_mapset(void)
Get current mapset name.
Definition: mapset.c:33
const char * name
Definition: named_colr.c:7
int IL_resample_output_2d(struct interp_params *params, double zmin, double zmax, double zminac, double zmaxac, double c1min, double c1max, double c2min, double c2max, double gmin, double gmax, double ertot, char *input, double *dnorm, struct Cell_head *outhd, struct Cell_head *winhd, char *smooth, int n_points)
Definition: resout2d.c:45
#define MULT
Definition: resout2d.c:14
double zmult
Definition: interpf.h:70
FILE * Tmp_fd_xx
Definition: interpf.h:93
FILE * Tmp_fd_xy
Definition: interpf.h:94
char * pcurv
Definition: interpf.h:85
FILE * Tmp_fd_yy
Definition: interpf.h:94
double fi
Definition: interpf.h:80
FILE * Tmp_fd_dx
Definition: interpf.h:92
FILE * Tmp_fd_z
Definition: interpf.h:92
char * tcurv
Definition: interpf.h:85
char * mcurv
Definition: interpf.h:85
FILE * Tmp_fd_dy
Definition: interpf.h:93
char * aspect
Definition: interpf.h:84
char * elev
Definition: interpf.h:84
char * slope
Definition: interpf.h:84