1 /*****************************************************************************
2 *
3 * MODULE: i.evapo.pt
4 * AUTHOR: Yann Chemin yann.chemin@gmail.com
5 *
6 * PURPOSE: To estimate the daily evapotranspiration by means
7 * of Prestley and Taylor method (1972).
8 *
9 * COPYRIGHT: (C) 2007-2011 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * Licence (>=2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 ***************************************************************************/
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <string.h>
20 #include <math.h>
21 #include <grass/gis.h>
22 #include <grass/raster.h>
23 #include <grass/glocale.h>
24
25 double pt_delta(double tempka);
26 double pt_ghamma(double tempka, double p_atm);
27 double pt_daily_et(double alpha_pt, double delta_pt, double ghamma_pt,
28 double rnet, double g0, double tempka);
29
main(int argc,char * argv[])30 int main(int argc, char *argv[])
31 {
32 /* buffer for input-output rasters */
33 void *inrast_TEMPKA, *inrast_PATM, *inrast_RNET, *inrast_G0;
34 DCELL *outrast;
35
36 /* pointers to input-output raster files */
37 int infd_TEMPKA, infd_PATM, infd_RNET, infd_G0;
38 int outfd;
39
40 /* names of input-output raster files */
41 char *RNET, *TEMPKA, *PATM, *G0;
42 char *ETa;
43
44 /* input-output cell values */
45 DCELL d_tempka, d_pt_patm, d_rnet, d_g0;
46 DCELL d_pt_alpha, d_pt_delta, d_pt_ghamma, d_daily_et;
47
48 /* region information and handler */
49 struct Cell_head cellhd;
50 int nrows, ncols;
51 int row, col;
52
53 /* parser stuctures definition */
54 struct GModule *module;
55 struct Option *input_RNET, *input_TEMPKA, *input_PATM, *input_G0,
56 *input_PT;
57 struct Option *output;
58 struct Flag *zero;
59 struct Colors color;
60 struct History history;
61
62 /* Initialize the GIS calls */
63 G_gisinit(argv[0]);
64
65 module = G_define_module();
66 G_add_keyword(_("imagery"));
67 G_add_keyword(_("evapotranspiration"));
68 module->description =
69 _("Computes evapotranspiration calculation "
70 "Priestley and Taylor formulation, 1972.");
71
72 /* Define different options */
73 input_RNET = G_define_standard_option(G_OPT_R_INPUT);
74 input_RNET->key = "net_radiation";
75 input_RNET->description = _("Name of input net radiation raster map [W/m2]");
76
77 input_G0 = G_define_standard_option(G_OPT_R_INPUT);
78 input_G0->key = "soil_heatflux";
79 input_G0->description = _("Name of input soil heat flux raster map [W/m2]");
80
81 input_TEMPKA = G_define_standard_option(G_OPT_R_INPUT);
82 input_TEMPKA->key = "air_temperature";
83 input_TEMPKA->description = _("Name of input air temperature raster map [K]");
84
85 input_PATM = G_define_standard_option(G_OPT_R_INPUT);
86 input_PATM->key = "atmospheric_pressure";
87 input_PATM->description = _("Name of input atmospheric pressure raster map [millibars]");
88
89 input_PT = G_define_option();
90 input_PT->key = "priestley_taylor_coeff";
91 input_PT->type = TYPE_DOUBLE;
92 input_PT->required = YES;
93 input_PT->description = _("Priestley-Taylor coefficient");
94 input_PT->answer = "1.26";
95
96 output = G_define_standard_option(G_OPT_R_OUTPUT);
97 output->description = _("Name of output evapotranspiration raster map [mm/d]");
98
99 /* Define the different flags */
100 zero = G_define_flag();
101 zero->key = 'z';
102 zero->description = _("Set negative ETa to zero");
103
104 if (G_parser(argc, argv))
105 exit(EXIT_FAILURE);
106
107 /* get entered parameters */
108 RNET = input_RNET->answer;
109 TEMPKA = input_TEMPKA->answer;
110 PATM = input_PATM->answer;
111 G0 = input_G0->answer;
112 d_pt_alpha = atof(input_PT->answer);
113
114 ETa = output->answer;
115
116 /* open pointers to input raster files */
117 infd_RNET = Rast_open_old(RNET, "");
118 infd_TEMPKA = Rast_open_old(TEMPKA, "");
119 infd_PATM = Rast_open_old(PATM, "");
120 infd_G0 = Rast_open_old(G0, "");
121
122 /* read headers of raster files */
123 Rast_get_cellhd(RNET, "", &cellhd);
124 Rast_get_cellhd(TEMPKA, "", &cellhd);
125 Rast_get_cellhd(PATM, "", &cellhd);
126 Rast_get_cellhd(G0, "", &cellhd);
127
128 /* Allocate input buffer */
129 inrast_RNET = Rast_allocate_d_buf();
130 inrast_TEMPKA = Rast_allocate_d_buf();
131 inrast_PATM = Rast_allocate_d_buf();
132 inrast_G0 = Rast_allocate_d_buf();
133
134 /* get rows and columns number of the current region */
135 nrows = Rast_window_rows();
136 ncols = Rast_window_cols();
137
138 /* allocate output buffer */
139 outrast = Rast_allocate_d_buf();
140
141 /* open pointers to output raster files */
142 outfd = Rast_open_new(ETa, DCELL_TYPE);
143
144 /* start the loop through cells */
145 for (row = 0; row < nrows; row++) {
146
147 G_percent(row, nrows, 2);
148 /* read input raster row into line buffer */
149 Rast_get_d_row(infd_RNET, inrast_RNET, row);
150 Rast_get_d_row(infd_TEMPKA, inrast_TEMPKA, row);
151 Rast_get_d_row(infd_PATM, inrast_PATM, row);
152 Rast_get_d_row(infd_G0, inrast_G0, row);
153
154 for (col = 0; col < ncols; col++) {
155 /* read current cell from line buffer */
156 d_rnet = ((DCELL *) inrast_RNET)[col];
157 d_tempka = ((DCELL *) inrast_TEMPKA)[col];
158 d_pt_patm = ((DCELL *) inrast_PATM)[col];
159 d_g0 = ((DCELL *) inrast_G0)[col];
160
161 /*Delta_pt and Ghamma_pt */
162 d_pt_delta = pt_delta(d_tempka);
163 d_pt_ghamma = pt_ghamma(d_tempka, d_pt_patm);
164
165 /*Calculate ET */
166 d_daily_et =
167 pt_daily_et(d_pt_alpha, d_pt_delta, d_pt_ghamma, d_rnet, d_g0,
168 d_tempka);
169 if (zero->answer && d_daily_et < 0)
170 d_daily_et = 0.0;
171
172 /* write calculated ETP to output line buffer */
173 outrast[col] = d_daily_et;
174 }
175
176 /* write output line buffer to output raster file */
177 Rast_put_d_row(outfd, outrast);
178 }
179 /* free buffers and close input maps */
180
181 G_free(inrast_RNET);
182 G_free(inrast_TEMPKA);
183 G_free(inrast_PATM);
184 G_free(inrast_G0);
185 Rast_close(infd_RNET);
186 Rast_close(infd_TEMPKA);
187 Rast_close(infd_PATM);
188 Rast_close(infd_G0);
189
190 /* generate color table between -20 and 20 */
191 Rast_make_rainbow_colors(&color, -20, 20);
192 Rast_write_colors(ETa, G_mapset(), &color);
193
194 Rast_short_history(ETa, "raster", &history);
195 Rast_command_history(&history);
196 Rast_write_history(ETa, &history);
197
198 /* free buffers and close output map */
199 G_free(outrast);
200 Rast_close(outfd);
201
202 return (EXIT_SUCCESS);
203 }
204