1 /****************************************************************************
2 *
3 * MODULE: i.eb.eta
4 * AUTHOR(S): Yann Chemin - yann.chemin@gmail.com
5 * PURPOSE: Calculates the actual evapotranspiration for diurnal period
6 * as seen in Bastiaanssen (1995)
7 *
8 * COPYRIGHT: (C) 2002-2011 by the GRASS Development Team
9 *
10 * This program is free software under the GNU General Public
11 * License (>=v2). Read the file COPYING that comes with GRASS
12 * for details.
13 *
14 *****************************************************************************/
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <grass/gis.h>
20 #include <grass/raster.h>
21 #include <grass/glocale.h>
22
23 double et_a(double r_net_day, double evap_fr, double tempk);
24
main(int argc,char * argv[])25 int main(int argc, char *argv[])
26 {
27 int nrows, ncols;
28 int row, col;
29 struct GModule *module;
30 struct Option *input1, *input2, *input3, *output1;
31 struct History history; /*metadata */
32 char *result1; /*output raster name */
33 int infd_rnetday, infd_evapfr, infd_tempk;
34 int outfd1;
35 char *rnetday, *evapfr, *tempk;
36 void *inrast_rnetday, *inrast_evapfr, *inrast_tempk;
37
38 DCELL * outrast1;
39 G_gisinit(argv[0]);
40
41 module = G_define_module();
42 G_add_keyword(_("imagery"));
43 G_add_keyword(_("energy balance"));
44 G_add_keyword(_("actual evapotranspiration"));
45 G_add_keyword(_("SEBAL"));
46 module->description =
47 _("Actual evapotranspiration for diurnal period (Bastiaanssen, 1995).");
48
49 /* Define the different options */
50 input1 = G_define_standard_option(G_OPT_R_INPUT);
51 input1->key = "netradiationdiurnal";
52 input1->description = _("Name of the diurnal net radiation map [W/m2]");
53 input1->answer = "rnetday";
54
55 input2 = G_define_standard_option(G_OPT_R_INPUT);
56 input2->key = "evaporativefraction";
57 input2->description = _("Name of the evaporative fraction map [-]");
58 input2->answer = "evapfr";
59
60 input3 = G_define_standard_option(G_OPT_R_INPUT);
61 input3->key = "temperature";
62 input3->description = _("Name of the surface skin temperature [K]");
63 input3->answer = "tempk";
64
65 output1 = G_define_standard_option(G_OPT_R_OUTPUT);
66 output1->description =
67 _("Name of the output actual evapotranspiration layer [mm/d]");
68
69 if (G_parser(argc, argv))
70 exit(EXIT_FAILURE);
71
72 rnetday = input1->answer;
73 evapfr = input2->answer;
74 tempk = input3->answer;
75 result1 = output1->answer;
76
77 infd_rnetday = Rast_open_old(rnetday, "");
78 inrast_rnetday = Rast_allocate_d_buf();
79
80 infd_evapfr = Rast_open_old(evapfr, "");
81 inrast_evapfr = Rast_allocate_d_buf();
82
83 infd_tempk = Rast_open_old(tempk, "");
84 inrast_tempk = Rast_allocate_d_buf();
85
86 nrows = Rast_window_rows();
87 ncols = Rast_window_cols();
88 outrast1 = Rast_allocate_d_buf();
89
90 outfd1 = Rast_open_new(result1, DCELL_TYPE);
91
92 /* Process pixels */
93 for (row = 0; row < nrows; row++)
94 {
95 DCELL d;
96 DCELL d_rnetday;
97 DCELL d_evapfr;
98 DCELL d_tempk;
99 G_percent(row, nrows, 2);
100
101 /* read input maps */
102 Rast_get_d_row(infd_rnetday,inrast_rnetday,row);
103 Rast_get_d_row(infd_evapfr,inrast_evapfr,row);
104 Rast_get_d_row(infd_tempk,inrast_tempk,row);
105
106 /*process the data */
107 for (col = 0; col < ncols; col++)
108 {
109 d_rnetday = ((DCELL *) inrast_rnetday)[col];
110 d_evapfr = ((DCELL *) inrast_evapfr)[col];
111 d_tempk = ((DCELL *) inrast_tempk)[col];
112 if (Rast_is_d_null_value(&d_rnetday) ||
113 Rast_is_d_null_value(&d_evapfr) ||
114 Rast_is_d_null_value(&d_tempk))
115 Rast_set_d_null_value(&outrast1[col], 1);
116 else {
117 d = et_a(d_rnetday, d_evapfr, d_tempk);
118 outrast1[col] = d;
119 }
120 }
121 Rast_put_d_row(outfd1,outrast1);
122 }
123 G_free(inrast_rnetday);
124 G_free(inrast_evapfr);
125 G_free(inrast_tempk);
126 Rast_close(infd_rnetday);
127 Rast_close(infd_evapfr);
128 Rast_close(infd_tempk);
129 G_free(outrast1);
130 Rast_close(outfd1);
131 Rast_short_history(result1, "raster", &history);
132 Rast_command_history(&history);
133 Rast_write_history(result1, &history);
134 exit(EXIT_SUCCESS);
135 }
136
137