output.c 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733
  1. /* output.c (simlib), 20.nov.2002, JH */
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <grass/gis.h>
  6. #include <grass/raster.h>
  7. #include <grass/bitmap.h>
  8. #include <grass/linkm.h>
  9. #include <grass/vector.h>
  10. #include <grass/glocale.h>
  11. #include <grass/waterglobs.h>
  12. static void output_walker_as_vector(int tt, int ndigit);
  13. /* This function was added by Soeren 8. Mar 2011 */
  14. /* It replaces the site walker output implementation */
  15. /* Only the 3d coordinates of the walker are stored. */
  16. void output_walker_as_vector(int tt, int ndigit)
  17. {
  18. char buf[GNAME_MAX + 10];
  19. char *outwalk_time = NULL;
  20. double x, y, z;
  21. struct Map_info Out;
  22. struct line_pnts *Points;
  23. struct line_cats *Cats;
  24. int i;
  25. if (outwalk != NULL) {
  26. /* In case of time series we extent the output name with the time value */
  27. if (ts == 1) {
  28. G_snprintf(buf, sizeof(buf), "%s_%.*d", outwalk, ndigit, tt);
  29. outwalk_time = G_store(buf);
  30. Vect_open_new(&Out, outwalk_time, WITH_Z);
  31. G_message("Writing %i walker into vector file %s", nstack, outwalk_time);
  32. }
  33. else {
  34. Vect_open_new(&Out, outwalk, WITH_Z);
  35. G_message("Writing %i walker into vector file %s", nstack, outwalk);
  36. }
  37. Points = Vect_new_line_struct();
  38. Cats = Vect_new_cats_struct();
  39. for (i = 0; i < nstack; i++) {
  40. x = stack[i][0];
  41. y = stack[i][1];
  42. z = stack[i][2];
  43. Vect_reset_line(Points);
  44. Vect_reset_cats(Cats);
  45. Vect_cat_set(Cats, 1, i + 1);
  46. Vect_append_point(Points, x, y, z);
  47. Vect_write_line(&Out, GV_POINT, Points, Cats);
  48. }
  49. /* Close vector file */
  50. Vect_close(&Out);
  51. Vect_destroy_line_struct(Points);
  52. Vect_destroy_cats_struct(Cats);
  53. }
  54. return;
  55. }
  56. /* Soeren 8. Mar 2011 TODO:
  57. * This function needs to be refractured and splittet into smaller parts */
  58. int output_data(int tt, double ft)
  59. {
  60. FCELL *depth_cell, *disch_cell, *err_cell;
  61. FCELL *conc_cell, *flux_cell, *erdep_cell;
  62. int depth_fd, disch_fd, err_fd;
  63. int conc_fd, flux_fd, erdep_fd;
  64. int i, iarc, j;
  65. float gsmax = 0, dismax = 0., gmax = 0., ermax = -1.e+12, ermin = 1.e+12;
  66. struct Colors colors;
  67. struct History hist, hist1; /* hist2, hist3, hist4, hist5 */
  68. char *depth0 = NULL, *disch0 = NULL, *err0 = NULL;
  69. char *conc0 = NULL, *flux0 = NULL;
  70. char *erdep0 = NULL;
  71. const char *mapst = NULL;
  72. char *type;
  73. char buf[GNAME_MAX + 10];
  74. int ndigit;
  75. FCELL dat1, dat2;
  76. float a1, a2;
  77. ndigit = 2;
  78. /* more compact but harder to read:
  79. ndigit = (int)floor(log10(timesec)) + 2 */
  80. if (timesec >= 10)
  81. ndigit = 3;
  82. if (timesec >= 100)
  83. ndigit = 4;
  84. if (timesec >= 1000)
  85. ndigit = 5;
  86. if (timesec >= 10000)
  87. ndigit = 6;
  88. /* Write the output walkers */
  89. output_walker_as_vector(tt, ndigit);
  90. Rast_set_window(&cellhd);
  91. if (my != Rast_window_rows())
  92. G_fatal_error("OOPS: rows changed from %d to %d", mx,
  93. Rast_window_rows());
  94. if (mx != Rast_window_cols())
  95. G_fatal_error("OOPS: cols changed from %d to %d", my,
  96. Rast_window_cols());
  97. if (depth) {
  98. depth_cell = Rast_allocate_f_buf();
  99. if (ts == 1) {
  100. G_snprintf(buf, sizeof(buf), "%s.%.*d", depth, ndigit, tt);
  101. depth0 = G_store(buf);
  102. depth_fd = Rast_open_fp_new(depth0);
  103. }
  104. else
  105. depth_fd = Rast_open_fp_new(depth);
  106. }
  107. if (disch) {
  108. disch_cell = Rast_allocate_f_buf();
  109. if (ts == 1) {
  110. G_snprintf(buf, sizeof(buf),"%s.%.*d", disch, ndigit, tt);
  111. disch0 = G_store(buf);
  112. disch_fd = Rast_open_fp_new(disch0);
  113. }
  114. else
  115. disch_fd = Rast_open_fp_new(disch);
  116. }
  117. if (err) {
  118. err_cell = Rast_allocate_f_buf();
  119. if (ts == 1) {
  120. G_snprintf(buf, sizeof(buf), "%s.%.*d", err, ndigit, tt);
  121. err0 = G_store(buf);
  122. err_fd = Rast_open_fp_new(err0);
  123. }
  124. else
  125. err_fd = Rast_open_fp_new(err);
  126. }
  127. if (conc) {
  128. conc_cell = Rast_allocate_f_buf();
  129. if (ts == 1) {
  130. G_snprintf(buf, sizeof(buf), "%s.%.*d", conc, ndigit, tt);
  131. conc0 = G_store(buf);
  132. conc_fd = Rast_open_fp_new(conc0);
  133. }
  134. else
  135. conc_fd = Rast_open_fp_new(conc);
  136. }
  137. if (flux) {
  138. flux_cell = Rast_allocate_f_buf();
  139. if (ts == 1) {
  140. G_snprintf(buf, sizeof(buf), "%s.%.*d", flux, ndigit, tt);
  141. flux0 = G_store(buf);
  142. flux_fd = Rast_open_fp_new(flux0);
  143. }
  144. else
  145. flux_fd = Rast_open_fp_new(flux);
  146. }
  147. if (erdep) {
  148. erdep_cell = Rast_allocate_f_buf();
  149. if (ts == 1) {
  150. G_snprintf(buf, sizeof(buf), "%s.%.*d", erdep, ndigit, tt);
  151. erdep0 = G_store(buf);
  152. erdep_fd = Rast_open_fp_new(erdep0);
  153. }
  154. else
  155. erdep_fd = Rast_open_fp_new(erdep);
  156. }
  157. for (iarc = 0; iarc < my; iarc++) {
  158. i = my - iarc - 1;
  159. if (depth) {
  160. for (j = 0; j < mx; j++) {
  161. if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
  162. Rast_set_f_null_value(depth_cell + j, 1);
  163. else {
  164. a1 = pow(gama[i][j], 3. / 5.);
  165. depth_cell[j] = (FCELL) a1; /* add conv? */
  166. gmax = amax1(gmax, a1);
  167. }
  168. }
  169. Rast_put_f_row(depth_fd, depth_cell);
  170. }
  171. if (disch) {
  172. for (j = 0; j < mx; j++) {
  173. if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
  174. cchez[i][j] == UNDEF)
  175. Rast_set_f_null_value(disch_cell + j, 1);
  176. else {
  177. a2 = step * gama[i][j] * cchez[i][j]; /* cchez incl. sqrt(sinsl) */
  178. disch_cell[j] = (FCELL) a2; /* add conv? */
  179. dismax = amax1(dismax, a2);
  180. }
  181. }
  182. Rast_put_f_row(disch_fd, disch_cell);
  183. }
  184. if (err) {
  185. for (j = 0; j < mx; j++) {
  186. if (zz[i][j] == UNDEF || gammas[i][j] == UNDEF)
  187. Rast_set_f_null_value(err_cell + j, 1);
  188. else {
  189. err_cell[j] = (FCELL) gammas[i][j];
  190. gsmax = amax1(gsmax, gammas[i][j]); /* add conv? */
  191. }
  192. }
  193. Rast_put_f_row(err_fd, err_cell);
  194. }
  195. if (conc) {
  196. for (j = 0; j < mx; j++) {
  197. if (zz[i][j] == UNDEF || gama[i][j] == UNDEF)
  198. Rast_set_f_null_value(conc_cell + j, 1);
  199. else {
  200. conc_cell[j] = (FCELL) gama[i][j];
  201. /* gsmax = amax1(gsmax, gama[i][j]); */
  202. }
  203. }
  204. Rast_put_f_row(conc_fd, conc_cell);
  205. }
  206. if (flux) {
  207. for (j = 0; j < mx; j++) {
  208. if (zz[i][j] == UNDEF || gama[i][j] == UNDEF ||
  209. slope[i][j] == UNDEF)
  210. Rast_set_f_null_value(flux_cell + j, 1);
  211. else {
  212. a2 = gama[i][j] * slope[i][j];
  213. flux_cell[j] = (FCELL) a2;
  214. dismax = amax1(dismax, a2);
  215. }
  216. }
  217. Rast_put_f_row(flux_fd, flux_cell);
  218. }
  219. if (erdep) {
  220. for (j = 0; j < mx; j++) {
  221. if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
  222. Rast_set_f_null_value(erdep_cell + j, 1);
  223. else {
  224. erdep_cell[j] = (FCELL) er[i][j];
  225. ermax = amax1(ermax, er[i][j]);
  226. ermin = amin1(ermin, er[i][j]);
  227. }
  228. }
  229. Rast_put_f_row(erdep_fd, erdep_cell);
  230. }
  231. }
  232. if (depth)
  233. Rast_close(depth_fd);
  234. if (disch)
  235. Rast_close(disch_fd);
  236. if (err)
  237. Rast_close(err_fd);
  238. if (conc)
  239. Rast_close(conc_fd);
  240. if (flux)
  241. Rast_close(flux_fd);
  242. if (erdep)
  243. Rast_close(erdep_fd);
  244. if (depth) {
  245. Rast_init_colors(&colors);
  246. dat1 = (FCELL) 0.;
  247. dat2 = (FCELL) 0.001;
  248. Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
  249. &colors);
  250. dat1 = dat2;
  251. dat2 = (FCELL) 0.05;
  252. Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
  253. &colors);
  254. dat1 = dat2;
  255. dat2 = (FCELL) 0.1;
  256. Rast_add_f_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
  257. &colors);
  258. dat1 = dat2;
  259. dat2 = (FCELL) 0.5;
  260. Rast_add_f_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
  261. &colors);
  262. dat1 = dat2;
  263. dat2 = (FCELL) gmax;
  264. Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
  265. if (ts == 1) {
  266. if ((mapst = G_find_file("fcell", depth0, "")) == NULL)
  267. G_fatal_error(_("FP raster map <%s> not found"), depth0);
  268. Rast_write_colors(depth0, mapst, &colors);
  269. Rast_quantize_fp_map_range(depth0, mapst, 0., (FCELL) gmax, 0,
  270. (CELL) gmax);
  271. Rast_free_colors(&colors);
  272. }
  273. else {
  274. if ((mapst = G_find_file("fcell", depth, "")) == NULL)
  275. G_fatal_error(_("FP raster map <%s> not found"), depth);
  276. Rast_write_colors(depth, mapst, &colors);
  277. Rast_quantize_fp_map_range(depth, mapst, 0., (FCELL) gmax, 0,
  278. (CELL) gmax);
  279. Rast_free_colors(&colors);
  280. }
  281. }
  282. if (disch) {
  283. Rast_init_colors(&colors);
  284. dat1 = (FCELL) 0.;
  285. dat2 = (FCELL) 0.0005;
  286. Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
  287. &colors);
  288. dat1 = dat2;
  289. dat2 = (FCELL) 0.005;
  290. Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 0, 255, 255,
  291. &colors);
  292. dat1 = dat2;
  293. dat2 = (FCELL) 0.05;
  294. Rast_add_f_color_rule(&dat1, 0, 255, 255, &dat2, 0, 127, 255,
  295. &colors);
  296. dat1 = dat2;
  297. dat2 = (FCELL) 0.1;
  298. Rast_add_f_color_rule(&dat1, 0, 127, 255, &dat2, 0, 0, 255,
  299. &colors);
  300. dat1 = dat2;
  301. dat2 = (FCELL) dismax;
  302. Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 0, 0, &colors);
  303. if (ts == 1) {
  304. if ((mapst = G_find_file("cell", disch0, "")) == NULL)
  305. G_fatal_error(_("Raster map <%s> not found"), disch0);
  306. Rast_write_colors(disch0, mapst, &colors);
  307. Rast_quantize_fp_map_range(disch0, mapst, 0., (FCELL) dismax, 0,
  308. (CELL) dismax);
  309. Rast_free_colors(&colors);
  310. }
  311. else {
  312. if ((mapst = G_find_file("cell", disch, "")) == NULL)
  313. G_fatal_error(_("Raster map <%s> not found"), disch);
  314. Rast_write_colors(disch, mapst, &colors);
  315. Rast_quantize_fp_map_range(disch, mapst, 0., (FCELL) dismax, 0,
  316. (CELL) dismax);
  317. Rast_free_colors(&colors);
  318. }
  319. }
  320. if (flux) {
  321. Rast_init_colors(&colors);
  322. dat1 = (FCELL) 0.;
  323. dat2 = (FCELL) 0.001;
  324. Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
  325. &colors);
  326. dat1 = dat2;
  327. dat2 = (FCELL) 0.1;
  328. Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
  329. &colors);
  330. dat1 = dat2;
  331. dat2 = (FCELL) 1.;
  332. Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 191, 127, 63,
  333. &colors);
  334. dat1 = dat2;
  335. dat2 = (FCELL) dismax;
  336. Rast_add_f_color_rule(&dat1, 191, 127, 63, &dat2, 0, 0, 0,
  337. &colors);
  338. if (ts == 1) {
  339. if ((mapst = G_find_file("cell", flux0, "")) == NULL)
  340. G_fatal_error(_("Raster map <%s> not found"), flux0);
  341. Rast_write_colors(flux0, mapst, &colors);
  342. Rast_quantize_fp_map_range(flux0, mapst, 0., (FCELL) dismax, 0,
  343. (CELL) dismax);
  344. Rast_free_colors(&colors);
  345. }
  346. else {
  347. if ((mapst = G_find_file("cell", flux, "")) == NULL)
  348. G_fatal_error(_("Raster map <%s> not found"), flux);
  349. Rast_write_colors(flux, mapst, &colors);
  350. Rast_quantize_fp_map_range(flux, mapst, 0., (FCELL) dismax, 0,
  351. (CELL) dismax);
  352. Rast_free_colors(&colors);
  353. }
  354. }
  355. if (erdep) {
  356. Rast_init_colors(&colors);
  357. dat1 = (FCELL) ermax;
  358. dat2 = (FCELL) 0.1;
  359. Rast_add_f_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
  360. dat1 = dat2;
  361. dat2 = (FCELL) 0.01;
  362. Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
  363. &colors);
  364. dat1 = dat2;
  365. dat2 = (FCELL) 0.0001;
  366. Rast_add_f_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
  367. &colors);
  368. dat1 = dat2;
  369. dat2 = (FCELL) 0.;
  370. Rast_add_f_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
  371. &colors);
  372. dat1 = dat2;
  373. dat2 = (FCELL) - 0.0001;
  374. Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
  375. &colors);
  376. dat1 = dat2;
  377. dat2 = (FCELL) - 0.01;
  378. Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
  379. &colors);
  380. dat1 = dat2;
  381. dat2 = (FCELL) - 0.1;
  382. Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
  383. &colors);
  384. dat1 = dat2;
  385. dat2 = (FCELL) ermin;
  386. Rast_add_f_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
  387. &colors);
  388. if (ts == 1) {
  389. if ((mapst = G_find_file("cell", erdep0, "")) == NULL)
  390. G_fatal_error(_("Raster map <%s> not found"), erdep0);
  391. Rast_write_colors(erdep0, mapst, &colors);
  392. Rast_quantize_fp_map_range(erdep0, mapst, (FCELL) ermin,
  393. (FCELL) ermax, (CELL) ermin,
  394. (CELL) ermax);
  395. Rast_free_colors(&colors);
  396. type = "raster";
  397. Rast_short_history(erdep0, type, &hist1);
  398. Rast_append_format_history(
  399. &hist1, "The sediment flux file is %s", flux0);
  400. Rast_command_history(&hist1);
  401. Rast_write_history(erdep0, &hist1);
  402. }
  403. else {
  404. if ((mapst = G_find_file("cell", erdep, "")) == NULL)
  405. G_fatal_error(_("Raster map <%s> not found"), erdep);
  406. Rast_write_colors(erdep, mapst, &colors);
  407. Rast_quantize_fp_map_range(erdep, mapst, (FCELL) ermin,
  408. (FCELL) ermax, (CELL) ermin,
  409. (CELL) ermax);
  410. Rast_free_colors(&colors);
  411. type = "raster";
  412. Rast_short_history(erdep, type, &hist1);
  413. Rast_append_format_history(
  414. &hist1, "The sediment flux file is %s", flux);
  415. Rast_command_history(&hist1);
  416. Rast_write_history(erdep, &hist1);
  417. }
  418. }
  419. /* history section */
  420. if (depth) {
  421. type = "raster";
  422. if (ts == 0) {
  423. mapst = G_find_file("cell", depth, "");
  424. if (mapst == NULL) {
  425. G_warning(_("Raster map <%s> not found"), depth);
  426. return -1;
  427. }
  428. Rast_short_history(depth, type, &hist);
  429. }
  430. else
  431. Rast_short_history(depth0, type, &hist);
  432. Rast_append_format_history(
  433. &hist, "init.walk=%d, maxwalk=%d, remaining walkers=%d",
  434. nwalk, maxwa, nwalka);
  435. Rast_append_format_history(
  436. &hist, "duration (sec.)=%d, time-serie iteration=%d",
  437. timesec, tt);
  438. Rast_append_format_history(
  439. &hist, "written deltap=%f, mean vel.=%f",
  440. deltap, vmean);
  441. Rast_append_format_history(
  442. &hist, "mean source (si)=%e, mean infil=%e",
  443. si0, infmean);
  444. Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
  445. elevin, dxin, dyin);
  446. Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s",
  447. rain, infil, manin);
  448. Rast_command_history(&hist);
  449. if (ts == 1)
  450. Rast_write_history(depth0, &hist);
  451. else
  452. Rast_write_history(depth, &hist);
  453. }
  454. if (disch) {
  455. type = "raster";
  456. if (ts == 0) {
  457. mapst = G_find_file("cell", disch, "");
  458. if (mapst == NULL)
  459. G_fatal_error(_("Raster map <%s> not found"), disch);
  460. Rast_short_history(disch, type, &hist);
  461. }
  462. else
  463. Rast_short_history(disch0, type, &hist);
  464. Rast_append_format_history(
  465. &hist,"init.walkers=%d, maxwalk=%d, rem. walkers=%d",
  466. nwalk, maxwa, nwalka);
  467. Rast_append_format_history(
  468. &hist, "duration (sec.)=%d, time-serie iteration=%d",
  469. timesec, tt);
  470. Rast_append_format_history(
  471. &hist, "written deltap=%f, mean vel.=%f",
  472. deltap, vmean);
  473. Rast_append_format_history(
  474. &hist, "mean source (si)=%e, mean infil=%e",
  475. si0, infmean);
  476. Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
  477. elevin, dxin, dyin);
  478. Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s",
  479. rain, infil, manin);
  480. Rast_command_history(&hist);
  481. if (ts == 1)
  482. Rast_write_history(disch0, &hist);
  483. else
  484. Rast_write_history(disch, &hist);
  485. }
  486. if (flux) {
  487. type = "raster";
  488. if (ts == 0) {
  489. mapst = G_find_file("cell", flux, "");
  490. if (mapst == NULL)
  491. G_fatal_error(_("Raster map <%s> not found"), flux);
  492. Rast_short_history(flux, type, &hist);
  493. }
  494. else
  495. Rast_short_history(flux0, type, &hist);
  496. Rast_append_format_history(
  497. &hist, "init.walk=%d, maxwalk=%d, remaining walkers=%d",
  498. nwalk, maxwa, nwalka);
  499. Rast_append_format_history(
  500. &hist, "duration (sec.)=%d, time-serie iteration=%d",
  501. timesec, tt);
  502. Rast_append_format_history(
  503. &hist, "written deltap=%f, mean vel.=%f",
  504. deltap, vmean);
  505. Rast_append_format_history(
  506. &hist, "mean source (si)=%f", si0);
  507. Rast_format_history(&hist, HIST_DATSRC_1, "input files: %s %s %s",
  508. wdepth, dxin, dyin);
  509. Rast_format_history(&hist, HIST_DATSRC_2, "input files: %s %s %s %s",
  510. manin, detin, tranin, tauin);
  511. Rast_command_history(&hist);
  512. if (ts == 1)
  513. Rast_write_history(flux0, &hist);
  514. else
  515. Rast_write_history(flux, &hist);
  516. }
  517. return 1;
  518. }
  519. int output_et()
  520. {
  521. FCELL *tc_cell, *et_cell;
  522. int tc_fd, et_fd;
  523. int i, iarc, j;
  524. float etmax = -1.e+12, etmin = 1.e+12;
  525. float trc;
  526. struct Colors colors;
  527. const char *mapst = NULL;
  528. /* char buf[GNAME_MAX + 10]; */
  529. FCELL dat1, dat2;
  530. /* float a1,a2; */
  531. Rast_set_window(&cellhd);
  532. if (et) {
  533. et_cell = Rast_allocate_f_buf();
  534. /* if (ts == 1) {
  535. sprintf(buf, "%s.%.*d", et, ndigit, tt);
  536. et0 = G_store(buf);
  537. et_fd = Rast_open_fp_new(et0);
  538. }
  539. else */
  540. et_fd = Rast_open_fp_new(et);
  541. }
  542. if (tc) {
  543. tc_cell = Rast_allocate_f_buf();
  544. /* if (ts == 1) {
  545. sprintf(buf, "%s.%.*d", tc, ndigit, tt);
  546. tc0 = G_store(buf);
  547. tc_fd = Rast_open_fp_new(tc0);
  548. }
  549. else */
  550. tc_fd = Rast_open_fp_new(tc);
  551. }
  552. if (my != Rast_window_rows())
  553. G_fatal_error("OOPS: rows changed from %d to %d", mx,
  554. Rast_window_rows());
  555. if (mx != Rast_window_cols())
  556. G_fatal_error("OOPS: cols changed from %d to %d", my,
  557. Rast_window_cols());
  558. for (iarc = 0; iarc < my; iarc++) {
  559. i = my - iarc - 1;
  560. if (et) {
  561. for (j = 0; j < mx; j++) {
  562. if (zz[i][j] == UNDEF || er[i][j] == UNDEF)
  563. Rast_set_f_null_value(et_cell + j, 1);
  564. else {
  565. et_cell[j] = (FCELL) er[i][j]; /* add conv? */
  566. etmax = amax1(etmax, er[i][j]);
  567. etmin = amin1(etmin, er[i][j]);
  568. }
  569. }
  570. Rast_put_f_row(et_fd, et_cell);
  571. }
  572. if (tc) {
  573. for (j = 0; j < mx; j++) {
  574. if (zz[i][j] == UNDEF || sigma[i][j] == UNDEF ||
  575. si[i][j] == UNDEF)
  576. Rast_set_f_null_value(tc_cell + j, 1);
  577. else {
  578. if (sigma[i][j] == 0.)
  579. trc = 0.;
  580. else
  581. trc = si[i][j] / sigma[i][j];
  582. tc_cell[j] = (FCELL) trc;
  583. /* gsmax = amax1(gsmax, trc); */
  584. }
  585. }
  586. Rast_put_f_row(tc_fd, tc_cell);
  587. }
  588. }
  589. if (tc)
  590. Rast_close(tc_fd);
  591. if (et)
  592. Rast_close(et_fd);
  593. if (et) {
  594. Rast_init_colors(&colors);
  595. dat1 = (FCELL) etmax;
  596. dat2 = (FCELL) 0.1;
  597. Rast_add_f_color_rule(&dat1, 0, 0, 0, &dat2, 0, 0, 255, &colors);
  598. dat1 = dat2;
  599. dat2 = (FCELL) 0.01;
  600. Rast_add_f_color_rule(&dat1, 0, 0, 255, &dat2, 0, 191, 191,
  601. &colors);
  602. dat1 = dat2;
  603. dat2 = (FCELL) 0.0001;
  604. Rast_add_f_color_rule(&dat1, 0, 191, 191, &dat2, 170, 255, 255,
  605. &colors);
  606. dat1 = dat2;
  607. dat2 = (FCELL) 0.;
  608. Rast_add_f_color_rule(&dat1, 170, 255, 255, &dat2, 255, 255, 255,
  609. &colors);
  610. dat1 = dat2;
  611. dat2 = (FCELL) - 0.0001;
  612. Rast_add_f_color_rule(&dat1, 255, 255, 255, &dat2, 255, 255, 0,
  613. &colors);
  614. dat1 = dat2;
  615. dat2 = (FCELL) - 0.01;
  616. Rast_add_f_color_rule(&dat1, 255, 255, 0, &dat2, 255, 127, 0,
  617. &colors);
  618. dat1 = dat2;
  619. dat2 = (FCELL) - 0.1;
  620. Rast_add_f_color_rule(&dat1, 255, 127, 0, &dat2, 255, 0, 0,
  621. &colors);
  622. dat1 = dat2;
  623. dat2 = (FCELL) etmin;
  624. Rast_add_f_color_rule(&dat1, 255, 0, 0, &dat2, 255, 0, 255,
  625. &colors);
  626. /* if (ts == 1) {
  627. if ((mapst = G_find_file("cell", et0, "")) == NULL)
  628. G_fatal_error(_("Raster map <%s> not found"), et0);
  629. Rast_write_colors(et0, mapst, &colors);
  630. Rast_quantize_fp_map_range(et0, mapst, (FCELL)etmin, (FCELL)etmax, (CELL)etmin, (CELL)etmax);
  631. Rast_free_colors(&colors);
  632. }
  633. else { */
  634. if ((mapst = G_find_file("cell", et, "")) == NULL)
  635. G_fatal_error(_("Raster map <%s> not found"), et);
  636. Rast_write_colors(et, mapst, &colors);
  637. Rast_quantize_fp_map_range(et, mapst, (FCELL) etmin, (FCELL) etmax,
  638. (CELL) etmin, (CELL) etmax);
  639. Rast_free_colors(&colors);
  640. /* } */
  641. }
  642. return 1;
  643. }