Added i.evapo.PT with imagery/Makefile modif and GUI menu item

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@44977 15284696-431f-4ddb-bdfa-cd5b030d7da7

gui/wxpython/xml/menudata.xml

@@ -2786,6 +2786,13 @@
 	      <handler>OnMenuCmd</handler>
 	      <command>i.evapo.PM</command>
 	    </menuitem>
+	    <menuitem>
+	      <label>Priestley-Taylor Evapotranspiration</label>
+	      <help>reference evapotranspiration (Priestley and Taylor, 1972)</help>
+	      <keywords>imagery,reference evapotranspiration,Priestley-Taylor</keywords>
+	      <handler>OnMenuCmd</handler>
+	      <command>i.evapo.PT</command>
+	    </menuitem>
 	  </items>
 	</menu>
 	<separator />

imagery/Makefile

@@ -13,6 +13,7 @@ SUBDIRS = \
 	i.eb.netrad \
 	i.eb.soilheatflux \
 	i.evapo.PM \
+	i.evapo.PT \
 	i.evapo.time_integration \
 	i.emissivity \
 	i.find \

imagery/i.evapo.PT/Makefile

@@ -0,0 +1,16 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.evapo.PT
+
+LIBES = $(VECTLIB) $(RASTERLIB) $(GISLIB)
+DEPENDENCIES= $(VECTDEP) $(GISDEP) $(DISPLAYDEP) $(RASTERDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

imagery/i.evapo.PT/i.evapo.PT.html

@@ -0,0 +1,36 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.evapo.PT</EM> Calculates the diurnal evapotranspiration after Prestley and Taylor (1972). 
+The Priestley-Taylor model (Priestley and Taylor, 1972) is a modification of Penman’s more theoretical equation.
+
+<H2>NOTES</H2>
+RNETD optional output from i.evapo.potrad is giving good results as input for net radiation in this module.
+
+Alpha values:
+1.32 for estimates from vegetated areas as a result of the increase in surface roughness (Morton, 1983; Brutsaert and Stricker, 1979)
+1.26 is applicable in humid climates (De Bruin and Keijman, 1979; Stewart and Rouse, 1976; Shuttleworth and Calder, 1979), and temperate hardwood swamps (Munro, 1979)
+1.74 has been recommended for estimating potential evapotranspiration in more arid regions (ASCE, 1990). This Worked well in Greece with University of Thessaloniki.
+
+Alpha values extracted from:
+http://www.civil.uwaterloo.ca/Watflood/Manual/02_03_1.htm
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.evapo.PM.html">i.evapo.PM</A><br>
+<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
+<A HREF="i.eb.netrad.html">i.eb.netrad</A><br>
+<A HREF="i.eb.g0.html">i.eb.g0</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, GRASS Development team, 2007-08<BR>
+
+
+<p>
+<i>Last changed: $Date$</i>

imagery/i.evapo.PT/main.c

@@ -0,0 +1,215 @@
+/*****************************************************************************
+*
+* MODULE:	i.evapo.PT
+* AUTHOR:	Yann Chemin [email protected] 
+*
+* PURPOSE:	To estimate the daily evapotranspiration by means
+*		of Prestley and Taylor method (1972).
+*
+* COPYRIGHT:	(C) 2007-2011 by the GRASS Development Team
+*
+*		This program is free software under the GNU General Public
+*		Licence (>=2). Read the file COPYING that comes with GRASS
+*		for details.
+*
+***************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+double pt_delta(double tempka);
+double pt_ghamma(double tempka, double p_atm);
+double pt_daily_et(double alpha_pt, double delta_pt, double ghamma_pt,
+		   double rnet, double g0, double tempka);
+
+int main(int argc, char *argv[])
+{
+    /* buffer for input-output rasters */
+    void *inrast_TEMPKA, *inrast_PATM, *inrast_RNET, *inrast_G0;
+    DCELL *outrast;
+
+    /* pointers to input-output raster files */
+    int infd_TEMPKA, infd_PATM, infd_RNET, infd_G0;
+    int outfd;
+
+    /* mapsets for input raster files */
+    char *mapset_TEMPKA, *mapset_PATM, *mapset_RNET, *mapset_G0;
+
+    /* names of input-output raster files */
+    char *RNET, *TEMPKA, *PATM, *G0;
+    char *ETa;
+
+    /* input-output cell values */
+    DCELL d_tempka, d_pt_patm, d_rnet, d_g0;
+    DCELL d_pt_alpha, d_pt_delta, d_pt_ghamma, d_daily_et;
+
+    /* region informations and handler */
+    struct Cell_head cellhd;
+    int nrows, ncols;
+    int row, col;
+
+    /* parser stuctures definition */
+    struct GModule *module;
+    struct Option *input_RNET, *input_TEMPKA, *input_PATM, *input_G0,
+	*input_PT;
+    struct Option *output;
+    struct Flag *flag1, *zero;
+    struct Colors color;
+    struct History history;
+
+    /* Initialize the GIS calls */
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->description =
+	_("Evapotranspiration Calculation "
+	  "Prestley and Taylor formulation, 1972.");
+
+    /* Define different options */
+    input_RNET = G_define_standard_option(G_OPT_R_INPUT);
+    input_RNET->key = "RNET";
+    input_RNET->key_desc = "[W/m2]";
+    input_RNET->description = _("Name of Net Radiation raster map");
+
+    input_G0 = G_define_standard_option(G_OPT_R_INPUT);
+    input_G0->key = "G0";
+    input_G0->key_desc = "[W/m2]";
+    input_G0->description = _("Name of Soil Heat Flux raster map");
+
+    input_TEMPKA = G_define_standard_option(G_OPT_R_INPUT);
+    input_TEMPKA->key = "TEMPKA";
+    input_TEMPKA->key_desc = "[K]";
+    input_TEMPKA->description = _("Name of air temperature raster map");
+
+    input_PATM = G_define_standard_option(G_OPT_R_INPUT);
+    input_PATM->key = "PATM";
+    input_PATM->key_desc = "[millibars]";
+    input_PATM->description = _("Name of Atmospheric Pressure raster map");
+
+    input_PT = G_define_option();
+    input_PT->key = "PT";
+    input_PT->key_desc = "[-]";
+    input_PT->type = TYPE_DOUBLE;
+    input_PT->required = YES;
+    input_PT->gisprompt = "old,cell,raster";
+    input_PT->description = _("Prestley-Taylor Coefficient");
+    input_PT->answer = "1.26";
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+    output->key_desc = "[mm/d]";
+    output->description = _("Name of output Evapotranspiration layer");
+
+    /* Define the different flags */
+    zero = G_define_flag();
+    zero->key = 'z';
+    zero->description = _("set negative ETa to zero");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    /* get entered parameters */
+    RNET = input_RNET->answer;
+    TEMPKA = input_TEMPKA->answer;
+    PATM = input_PATM->answer;
+    G0 = input_G0->answer;
+    d_pt_alpha = atof(input_PT->answer);
+
+    ETa = output->answer;
+
+    /* check legal output name */
+    if (G_legal_filename(ETa) < 0)
+	G_fatal_error(_("[%s] is an illegal name"), ETa);
+
+    /* open pointers to input raster files */
+    infd_RNET = Rast_open_old(RNET, "");
+    infd_TEMPKA = Rast_open_old(TEMPKA, "");
+    infd_PATM = Rast_open_old(PATM, "");
+    infd_G0 = Rast_open_old(G0, "");
+
+    /* read headers of raster files */
+    Rast_get_cellhd(RNET, "", &cellhd);
+    Rast_get_cellhd(TEMPKA, "", &cellhd);
+    Rast_get_cellhd(PATM, "", &cellhd);
+    Rast_get_cellhd(G0, "", &cellhd);
+
+    /* Allocate input buffer */
+    inrast_RNET = Rast_allocate_d_buf();
+    inrast_TEMPKA = Rast_allocate_d_buf();
+    inrast_PATM = Rast_allocate_d_buf();
+    inrast_G0 = Rast_allocate_d_buf();
+
+    /* get rows and columns number of the current region */
+    nrows = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    /* allocate output buffer */
+    outrast = Rast_allocate_d_buf();
+
+    /* open pointers to output raster files */
+    outfd = Rast_open_new(ETa, DCELL_TYPE);
+
+    /* start the loop through cells */
+    for (row = 0; row < nrows; row++) {
+
+	G_percent(row, nrows, 2);
+	/* read input raster row into line buffer */
+	Rast_get_d_row(infd_RNET, inrast_RNET, row);
+	Rast_get_d_row(infd_TEMPKA, inrast_TEMPKA, row);
+	Rast_get_d_row(infd_PATM, inrast_PATM, row);
+	Rast_get_d_row(infd_G0, inrast_G0, row);
+
+	for (col = 0; col < ncols; col++) {
+	    /* read current cell from line buffer */
+            d_rnet = ((DCELL *) inrast_RNET)[col];
+            d_tempka = ((DCELL *) inrast_TEMPKA)[col];
+            d_pt_patm = ((DCELL *) inrast_PATM)[col];
+            d_g0 = ((DCELL *) inrast_G0)[col];
+
+	    /*Delta_pt and Ghamma_pt */
+	    d_pt_delta = pt_delta(d_tempka);
+	    d_pt_ghamma = pt_ghamma(d_tempka, d_pt_patm);
+
+	    /*Calculate ET */
+	    d_daily_et =
+		pt_daily_et(d_pt_alpha, d_pt_delta, d_pt_ghamma, d_rnet, d_g0,
+			    d_tempka);
+	    if (zero->answer && d_daily_et < 0)
+		d_daily_et = 0.0;
+
+	    /* write calculated ETP to output line buffer */
+	    outrast[col] = d_daily_et;
+	}
+
+	/* write output line buffer to output raster file */
+	Rast_put_d_row(outfd, outrast);
+    }
+    /* free buffers and close input maps */
+
+    G_free(inrast_RNET);
+    G_free(inrast_TEMPKA);
+    G_free(inrast_PATM);
+    G_free(inrast_G0);
+    Rast_close(infd_RNET);
+    Rast_close(infd_TEMPKA);
+    Rast_close(infd_PATM);
+    Rast_close(infd_G0);
+
+    /* generate color table between -20 and 20 */
+    Rast_make_rainbow_colors(&color, -20, 20);
+    Rast_write_colors(ETa, G_mapset(), &color);
+
+    Rast_short_history(ETa, "raster", &history);
+    Rast_command_history(&history);
+    Rast_write_history(ETa, &history);
+
+    /* free buffers and close output map */
+    G_free(outrast);
+    Rast_close(outfd);
+
+    return (EXIT_SUCCESS);
+}

imagery/i.evapo.PT/pt_daily_et.c

@@ -0,0 +1,24 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+/*Prestely and Taylor, 1972. */
+double pt_daily_et(double alpha_pt, double delta_pt, double ghamma_pt,
+		   double rnet, double g0, double tempka)
+{
+    double result, latentHv, t_celsius;
+    double roh_w = 1004.15;	/*mass density of water */
+    double vap_slope_ratio;
+
+    /*Latent Heat of vaporization (W/m2/d) */
+    t_celsius = tempka - 273.15;
+    latentHv = 86400 / ((2.501 - 0.002361 * t_celsius) * pow(10, 6));
+
+    /* Ratio of slope of saturation-vapour pressure Vs Temperature */
+    /* ghamma_pt = psychrometric constant */
+    vap_slope_ratio = delta_pt / (delta_pt + ghamma_pt);
+
+    /*(Rn-g0)/latentHv returns [-] */
+    result = (alpha_pt / roh_w) * vap_slope_ratio * (rnet - g0) / latentHv;
+    return result;
+}

imagery/i.evapo.PT/pt_delta.c

@@ -0,0 +1,17 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+/* Prestely and Taylor, 1972. 
+ * INPUT in Kelvin
+ */
+double pt_delta(double tempka)
+{
+    double a, b, result;
+
+    tempka -= 273.15;		/*Celsius */
+    b = tempka + 237.3;
+    a = (17.27 * tempka) / b;
+    result = 2504.0 * exp(a) / pow(b, 2);
+    return result;
+}

imagery/i.evapo.PT/pt_ghamma.c

@@ -0,0 +1,17 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+/* Prestely and Taylor, 1972. 
+ * INPUT in Kelvin
+ */
+double pt_ghamma(double tempka, double patm_pt)
+{
+    double a, result;
+    double Cp = 1004.16;
+
+    tempka -= 273.15;
+    a = 0.622 * pow(10, 7) * (2.501 - (2.361 * pow(10, -3) * tempka));
+    result = Cp * patm_pt / a;
+    return result;
+}