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