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