Added i.evapo.MH with updated imagery/Makefile and menu item

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

gui/wxpython/xml/menudata.xml

@@ -2780,6 +2780,13 @@
 	    </menuitem>
 	    <separator />
 	    <menuitem>
+	      <label>Hargreaves methods Evapotranspiration</label>
+	      <help>reference evapotranspiration (Hargreaves, modified-H, H and Samani)</help>
+	      <keywords>imagery,reference evapotranspiration,Hargreaves,Samani</keywords>
+	      <handler>OnMenuCmd</handler>
+	      <command>i.evapo.MH</command>
+	    </menuitem>
+	    <menuitem>
 	      <label>Penman-Monteith Evapotranspiration</label>
 	      <help>reference evapotranspiration (FAO 56 and Penman, 1972)</help>
 	      <keywords>imagery,reference evapotranspiration,FAO56,Penman-Monteith,HydroFOSS</keywords>

imagery/Makefile

@@ -12,6 +12,7 @@ SUBDIRS = \
 	i.eb.h_SEBAL01 \
 	i.eb.netrad \
 	i.eb.soilheatflux \
+	i.evapo.MH \
 	i.evapo.PM \
 	i.evapo.PT \
 	i.evapo.time_integration \

imagery/i.evapo.MH/Makefile

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

imagery/i.evapo.MH/i.evapo.MH.html

@@ -0,0 +1,33 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.evapo.MH</EM> Calculates the reference ET after Hargreaves (1985) and Modified Hargreaves (2001). 
+
+<H2>NOTES</H2>
+Hargreaves GL, Hargreaves GH, Riley JP, 1985. Agricultural benefits for Senegal River Basin. Journal of Irrigation and Drainange Engineering, ASCE, 111(2):113-124.
+
+Droogers P, Allen RG, 2002. Towards a simplified global reference evapotranspiration equation. Irrigation Science.
+Droogers, P., and R.G. Allen. 2002. Estimating reference evapotranspiration under inaccurate data conditions. Irrigation and Drainage Systems 16: 33-45.
+
+Hargreaves and Samani, 1985.
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.evapo.PT.html">i.evapo.PT</A><br>
+<A HREF="i.evapo.PM.html">i.evapo.PM</A><br>
+<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
+<A HREF="r.sun.html">r.sun</A><br>
+
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, GRASS Development team, 2007-2011<BR>
+
+
+<p>
+<i>Last changed: $Date$</i>

imagery/i.evapo.MH/main.c

@@ -0,0 +1,251 @@
+/*****************************************************************************
+*
+* MODULE:	i.evapo.MH
+* AUTHOR:	Yann Chemin yann.chemin@gmail.com 
+*
+* PURPOSE:	To estimate the reference evapotranspiration by means
+*		of Modified Hargreaves method (2001).
+*		Also has a switch for original Hargreaves (1985),
+*		and for Hargreaves-Samani (1985).
+*
+* 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 mh_original(double ra, double tavg, double tmax, double tmin,
+		   double p);
+double mh_eto(double ra, double tavg, double tmax, double tmin, double p);
+double mh_samani(double ra, double tavg, double tmax, double tmin);
+
+int main(int argc, char *argv[])
+{
+    /* buffer for input-output rasters */
+    void *inrast_TEMPKAVG, *inrast_TEMPKMIN, *inrast_TEMPKMAX, *inrast_RNET,
+	*inrast_P;
+    DCELL *outrast;
+
+    /* pointers to input-output raster files */
+    int infd_TEMPKAVG, infd_TEMPKMIN, infd_TEMPKMAX, infd_RNET, infd_P;
+    int outfd;
+
+    /* mapsets for input raster files */
+    char *mapset_TEMPKAVG, *mapset_TEMPKMIN, *mapset_TEMPKMAX, *mapset_RNET,
+	*mapset_P;
+
+    /* names of input-output raster files */
+    char *RNET, *TEMPKAVG, *TEMPKMIN, *TEMPKMAX, *P;
+    char *ETa;
+
+    /* input-output cell values */
+    DCELL d_tempkavg, d_tempkmin, d_tempkmax, d_rnet, d_p;
+    DCELL 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_TEMPKAVG, *input_TEMPKMIN;
+    struct Option *input_TEMPKMAX, *input_P;
+    struct Option *output;
+    struct Flag *zero, *original, *samani;
+    struct Colors color;
+    struct History history;
+
+    /* Initialize the GIS calls */
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->description =
+	_("Evapotranspiration Calculation "
+	  "Modified Hargreaves formulation, 2001."
+	  "Flag for Original Hargreaves (1985).");
+
+    /* Define different options */
+    input_RNET = G_define_standard_option(G_OPT_R_INPUT);
+    input_RNET->key = "rnetd";
+    input_RNET->key_desc = "[W/m2/d]";
+    input_RNET->description = _("Name of Diurnal Net Radiation raster map");
+
+    input_TEMPKAVG = G_define_standard_option(G_OPT_R_INPUT);
+    input_TEMPKAVG->key = "tempkavg";
+    input_TEMPKAVG->key_desc = "[C]";
+    input_TEMPKAVG->description = _("Name of avg air temperature raster map");
+
+    input_TEMPKMIN = G_define_standard_option(G_OPT_R_INPUT);
+    input_TEMPKMIN->key = "tempkmin";
+    input_TEMPKMIN->key_desc = "[C]";
+    input_TEMPKMIN->description = _("Name of min air temperature raster map");
+
+    input_TEMPKMAX = G_define_standard_option(G_OPT_R_INPUT);
+    input_TEMPKMAX->key = "TEMPKMAX";
+    input_TEMPKMAX->key_desc = "[C]";
+    input_TEMPKMAX->description = _("Name of max air temperature raster map");
+
+    input_P = G_define_standard_option(G_OPT_R_INPUT);
+    input_P->key = "p";
+    input_P->key_desc = "[mm/month]";
+    input_P->description =
+	_("Name of precipitation raster map, disabled if original Hargreaves (1985) is enabled.");
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+    output->key_desc = "[mm/d]";
+    output->description = _("Name of output Ref Evapotranspiration layer");
+
+    /* Define the different flags */
+    zero = G_define_flag();
+    zero->key = 'z';
+    zero->description = _("set negative ETa to zero");
+
+    original = G_define_flag();
+    original->key = 'h';
+    original->description = _("set to original Hargreaves (1985)");
+
+    samani = G_define_flag();
+    samani->key = 's';
+    samani->description = _("set to Hargreaves-Samani (1985)");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    /* get entered parameters */
+    RNET = input_RNET->answer;
+    TEMPKAVG = input_TEMPKAVG->answer;
+    TEMPKMIN = input_TEMPKMIN->answer;
+    TEMPKMAX = input_TEMPKMAX->answer;
+    P = input_P->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_TEMPKAVG = Rast_open_old(TEMPKAVG, "");
+    infd_TEMPKMIN = Rast_open_old(TEMPKMIN, "");
+    infd_TEMPKMAX = Rast_open_old(TEMPKMAX, "");
+    if (!original->answer) {
+        infd_P = Rast_open_old(P, "");
+    }
+    /* read headers of raster files */
+    Rast_get_cellhd(RNET, "", &cellhd);
+    Rast_get_cellhd(TEMPKAVG, "", &cellhd);
+    Rast_get_cellhd(TEMPKMIN, "", &cellhd);
+    Rast_get_cellhd(TEMPKMAX, "", &cellhd);
+    if (!original->answer) {
+	Rast_get_cellhd(P, "", &cellhd);
+    }
+    /* Allocate input buffer */
+    inrast_RNET = Rast_allocate_d_buf();
+    inrast_TEMPKAVG = Rast_allocate_d_buf();
+    inrast_TEMPKMIN = Rast_allocate_d_buf();
+    inrast_TEMPKMAX = Rast_allocate_d_buf();
+    if (!original->answer) {
+	inrast_P = 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_TEMPKAVG, inrast_TEMPKAVG, row);
+	Rast_get_d_row(infd_TEMPKMIN, inrast_TEMPKMIN, row);
+	Rast_get_d_row(infd_TEMPKMAX, inrast_TEMPKMAX, row);
+	if (!original->answer) {
+	    Rast_get_d_row(infd_P, inrast_P, row);
+	}
+	for (col = 0; col < ncols; col++) {
+	    /* read current cell from line buffer */
+            d_rnet = ((DCELL *) inrast_RNET)[col];
+            d_tempkavg = ((DCELL *) inrast_TEMPKAVG)[col];
+            d_tempkmin = ((DCELL *) inrast_TEMPKMIN)[col];
+            d_tempkmax = ((DCELL *) inrast_TEMPKMAX)[col];
+	    if (!original->answer) {
+		    d_p = ((DCELL *) inrast_P)[col];
+	    }
+	    if (Rast_is_d_null_value(&d_rnet) ||
+		Rast_is_d_null_value(&d_tempkavg) ||
+		Rast_is_d_null_value(&d_tempkmin) ||
+		Rast_is_d_null_value(&d_tempkmax) || Rast_is_d_null_value(&d_p)) {
+		Rast_set_d_null_value(&outrast[col], 1);
+	    }
+	    else {
+		if (original->answer) {
+		    d_daily_et =
+			mh_original(d_rnet, d_tempkavg, d_tempkmax,
+				    d_tempkmin, d_p);
+		}
+		else if (samani->answer) {
+		    d_daily_et =
+			mh_samani(d_rnet, d_tempkavg, d_tempkmax, d_tempkmin);
+		}
+		else {
+		    d_daily_et =
+			mh_eto(d_rnet, d_tempkavg, d_tempkmax, d_tempkmin,
+			       d_p);
+		}
+		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_TEMPKAVG);
+    G_free(inrast_TEMPKMIN);
+    G_free(inrast_TEMPKMAX);
+    if (!original->answer) {
+	G_free(inrast_P);
+    }
+    Rast_close(infd_RNET);
+    Rast_close(infd_TEMPKAVG);
+    Rast_close(infd_TEMPKMIN);
+    Rast_close(infd_TEMPKMAX);
+    if (!original->answer) {
+	Rast_close(infd_P);
+    }
+    /* 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.MH/mh_eto.c

@@ -0,0 +1,19 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+    /*Droogers and Allen, 2001. */
+    /*p is [mm/month] */
+double mh_eto(double ra, double tavg, double tmax, double tmin, double p)
+{
+    double td, result;
+
+    td = tmax - tmin;
+    if (tavg > 100.0) {
+	tavg = tavg - 273.15;	/*in case temperature is in Kelvin */
+    }
+    ra = ra * (84600.0 * 1000.0);	/*convert W -> MJ/d */
+    result =
+	0.0013 * 0.408 * ra * (tavg + 17.0) * pow((td - 0.0123 * p), 0.76);
+    return result;
+}

imagery/i.evapo.MH/mh_original.c

@@ -0,0 +1,17 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+    /*Hargreaves et al, 1985. */
+double mh_original(double ra, double tavg, double tmax, double tmin, double p)
+{
+    double td, result;
+
+    td = tmax - tmin;
+    if (tavg > 100.0) {
+	tavg = tavg - 273.15;	/*in case Temperature is in Kelvin */
+    }
+    ra = ra * (84600.0 * 1000.0);	/*convert W -> MJ/d */
+    result = 0.0023 * 0.408 * ra * (tavg + 17.8) * pow(td, 0.5);
+    return result;
+}

imagery/i.evapo.MH/mh_samani.c

@@ -0,0 +1,19 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+    /*Hargreaves-Samani, 1985. */
+double mh_samani(double ra, double tavg, double tmax, double tmin)
+{
+    double td, result;
+
+    td = tmax - tmin;
+    if (tavg > 100.0) {
+	tavg = tavg - 273.15;	/*in case Temperature is in Kelvin */
+    }
+    ra = ra * (84600.0 * 1000.0);	/* convert W -> MJ/d */
+    result =
+	0.0023 * 0.408 * ra * pow(td,
+				  0.5) * ((tmax + tmin) / 2 + 17.8) / 2.45;
+    return result;
+}