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