GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 
28 static 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 
45 int 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..."));
79  mapset = G_mapset();
80 
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)
434  &hist, "tension=%f, smoothing=%s",
435  params->fi * 1000. / (*dnorm), smooth);
436  else
438  &hist, "tension=%f", params->fi * 1000. / (*dnorm));
439 
441  &hist, "dnorm=%f, zmult=%f", *dnorm, params->zmult);
443  &hist, "KMAX=%d, KMIN=%d, errtotal=%f", params->kmax,
444  params->kmin, sqrt(ertot / n_points));
446  &hist, "zmin_data=%f, zmax_data=%f", zmin, zmax);
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 }
struct _Color_Info_ fixed
Definition: gis.h:679
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&#39;t touch) ...
Definition: find_file.c:247
void Rast_write_colors(const char *, const char *, struct Colors *)
Write map layer color table.
void Rast_write_history(const char *, struct History *)
Write raster history file.
Definition: history.c:158
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:75
2D/3D raster map header (used also for region)
Definition: gis.h:412
char * mcurv
Definition: interpf.h:84
FILE * Tmp_fd_xy
Definition: interpf.h:92
FILE * Tmp_fd_yy
Definition: interpf.h:92
struct _Color_Value_ low high
Definition: gis.h:633
struct _Color_Rule_ * next
Definition: gis.h:634
double amax1(double, double)
Definition: minmax.c:54
char * elev
Definition: interpf.h:84
char * pcurv
Definition: interpf.h:84
void Rast_free_history(struct History *)
Definition: history.c:317
#define NULL
Definition: ccmath.h:32
FILE * Tmp_fd_xx
Definition: interpf.h:92
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
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
FILE * Tmp_fd_dx
Definition: interpf.h:92
FCELL * Rast_allocate_f_output_buf(void)
Definition: alloc_cell.c:193
void Rast_set_output_window(struct Cell_head *)
Establishes &#39;window&#39; as the current working window for output.
char * aspect
Definition: interpf.h:84
void Rast_append_format_history(struct History *, const char *,...) __attribute__((format(printf
DCELL value
Definition: gis.h:625
Raster history info (metadata)
Definition: raster.h:180
Description of original data source (two lines)
Definition: raster.h:170
void Rast_put_f_row(int, const FCELL *)
Writes the next row for fcell file (FCELL version)
void G_fseek(FILE *, off_t, int)
Change the file position of the stream.
Definition: gis/seek.c:50
double fi
Definition: interpf.h:80
int Rast_open_fp_new(const char *)
Opens new fcell file in a database.
Definition: raster/open.c:499
float FCELL
Definition: gis.h:604
Definition: gis.h:665
void Rast_short_history(const char *, const char *, struct History *)
Initialize history structure.
Definition: history.c:226
double amin1(double, double)
Definition: minmax.c:67
char * tcurv
Definition: interpf.h:84
int cols
Number of columns for 2D data.
Definition: gis.h:431
FILE * Tmp_fd_dy
Definition: interpf.h:92
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:55
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:203
const char * G_mapset(void)
Get current mapset name.
Definition: gis/mapset.c:33
#define MULT
Definition: resout2d.c:14
int Rast_read_colors(const char *, const char *, struct Colors *)
Read color table of raster map.
FILE * Tmp_fd_z
Definition: interpf.h:92
void G_warning(const char *,...) __attribute__((format(printf
int CELL
Definition: gis.h:602
unsigned char red
Definition: gis.h:626
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
#define _(str)
Definition: glocale.h:10
unsigned char blu
Definition: gis.h:628
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:35
double zmult
Definition: interpf.h:70
struct _Color_Info_ modular
Definition: gis.h:680
void Rast_format_history(struct History *, int, const char *,...) __attribute__((format(printf
unsigned char grn
Definition: gis.h:627
const char * name
Definition: named_colr.c:7
void void G_verbose_message(const char *,...) __attribute__((format(printf
int rows
Number of rows for 2D data.
Definition: gis.h:427
void Rast_init_colors(struct Colors *)
Initialize color structure.
Definition: color_init.c:25
struct _Color_Rule_ * prev
Definition: gis.h:635
void Rast_close(int)
Close a raster map.
Definition: raster/close.c:99
char * slope
Definition: interpf.h:84
struct _Color_Rule_ * rules
Definition: gis.h:640