Added soil heat flux module

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

imagery/Makefile

@@ -5,6 +5,7 @@ SUBDIRS = \
 	i.atcorr \
 	i.cca \
 	i.cluster \
+	i.eb.g0 \
 	i.emissivity \
 	i.find \
 	i.gensig \

imagery/i.eb.g0/Makefile

@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.eb.g0
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

imagery/i.eb.g0/g0.c

@@ -0,0 +1,28 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+ double g_0(double bbalb, double ndvi, double tempk, double rnet,
+	     double time, int roerink) 
+{
+    double a, b, result;
+    double r0_coef;
+    if (time <= 9.0 || time > 15.0) 
+	r0_coef = 1.1;
+    else if (time > 9.0 && time <= 11.0)
+	r0_coef = 1.0;
+    else if (time > 11.0 && time <= 13.0)
+	r0_coef = 0.9;
+    else if (time > 13.0 && time <= 15.0) 
+	r0_coef = 1.0;
+    a = (0.0032 * (bbalb / r0_coef) +
+	  0.0062 * (bbalb / r0_coef) * (bbalb / r0_coef));
+    b = (1 - 0.978 * pow(ndvi, 4));
+    /* Spain (Bastiaanssen, 1995) */ 
+    result = (rnet * (tempk - 273.15) / bbalb) * a * b;
+    /* HAPEX-Sahel (Roerink, 1995) */ 
+    if (roerink)
+	result = result * 1.430 - 0.0845;
+    return result;
+}
+
+

imagery/i.eb.g0/i.eb.g0.html

@@ -0,0 +1,37 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.eb.g0</EM> calculates the soil heat flux approximation (g0) after [1]. Main reference for implementation is [3].
+It takes input of Albedo, NDVI, Surface Skin temperature, Net Radiation (see r.sun), time of satellite overpass, and a flag for the Roerink empirical modification from the HAPEX-Sahel experiment. 
+<H2>NOTES</H2>
+
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="r.sun.html">r.sun</A><br>
+<A HREF="i.albedo.html">i.albedo</A><br>
+<A HREF="i.emissivity.html">i.emissivity</A><br>
+<A HREF="i.eb.h_SEBAL01.html">i.eb.h_SEBAL01</A><br>
+<A HREF="i.eb.evapfr.html">i.eb.evapfr</A><br>
+</em>
+
+<H2>REFERENCES</H2>
+
+  <p>[1] Bastiaanssen, W.G.M., 1995.
+  Estimation of Land surface parameters by remote sensing under clear-sky conditions. PhD thesis, Wageningen University, Wageningen, The Netherlands.
+
+  <p>[2] Chemin Y., Alexandridis T.A., 2001. Improving spatial resolution of ET seasonal for irrigated rice in Zhanghe, China. Asian Journal of Geoinformatics. 5(1):3-11,2004. 
+
+  <p>[3] Alexandridis T.K., Cherif I., Chemin Y., Silleos N.G., Stavrinos E., Zalidis G.C. Integrated methodology for estimating water use in Mediterranean agricultural areas. Remote Sensing. -(-):,2009. (submitted))
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, Asian Institute of Technology, Thailand<BR>
+
+
+<p>
+<i>Last changed: $Date: 2006/10/08 21:30:42 $</i>

imagery/i.eb.g0/main.c

@@ -0,0 +1,193 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.eb.g0
+ * AUTHOR(S):    Yann Chemin - [email protected]
+ * PURPOSE:      Calculates an approximation of soil heat flux
+ *               as seen in Bastiaanssen (1995) using time of
+ *               satellite overpass.
+ *
+ * COPYRIGHT:    (C) 2006-2007 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *   	    	 License (>=v2). Read the file COPYING that comes with GRASS
+ *   	    	 for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+double g_0(double bbalb, double ndvi, double tempk, double rnet,
+	     double time, int roerink);
+int main(int argc, char *argv[]) 
+{
+    int nrows, ncols;
+    int row, col;
+    int roerink = 0;		/*Roerink Flag for HAPEX-Sahel conditions */
+    struct GModule *module;
+    struct Option *input1, *input2, *input3, *input4, *input5, *output1;
+    struct Flag *flag1;
+    struct History history;	/*metadata */
+    struct Colors colors;	/*metadata */
+    char *result;		/*output raster name */
+    int infd_albedo, infd_ndvi, infd_tempk, infd_rnet, infd_time;
+    int outfd;
+    char *albedo, *ndvi, *tempk, *rnet, *time;
+    void *inrast_albedo, *inrast_ndvi, *inrast_tempk, *inrast_rnet,
+	*inrast_time;
+    DCELL * outrast;
+
+    /************************************/ 
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    module->keywords = _("soil heat flux, energy balance, SEBAL");
+    module->description =
+	_("soil heat flux approximation (Bastiaanssen, 1995)");
+
+    /* Define the different options */ 
+    input1 = G_define_standard_option(G_OPT_R_INPUT);
+    input1->key = "albedo";
+    input1->description = _("Name of the Albedo map [0.0;1.0]");
+    input1->answer = "albedo";
+
+    input2 = G_define_standard_option(G_OPT_R_INPUT);
+    input2->key = "ndvi";
+    input2->description = _("Name of the ndvi map [-1.0;+1.0]");
+    input2->answer = "ndvi";
+
+    input3 = G_define_standard_option(G_OPT_R_INPUT);
+    input3->key = "tempk";
+    input3->description =
+	_("Name of the Surface temperature map [degree Kelvin]");
+    input3->answer = "tempk";
+
+    input4 = G_define_standard_option(G_OPT_R_INPUT);
+    input4->key = "rnet";
+    input4->description = _("Name of the Net Radiation map [W/m2]");
+    input4->answer = "rnet";
+
+    input5 = G_define_standard_option(G_OPT_R_INPUT);
+    input5->key = "time";
+    input5->description =
+	_("Name of the time of satellite overpass map [local UTC]");
+    input5->answer = "time";
+
+    output1 = G_define_standard_option(G_OPT_R_OUTPUT);
+    output1->key = "g0";
+    output1->description = _("Name of the output g0 layer");
+    output1->answer = "g0";
+
+    flag1 = G_define_flag();
+    flag1->key = 'r';
+    flag1->description =
+	_("HAPEX-Sahel empirical correction (Roerink, 1995)");
+
+    /********************/ 
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    albedo = input1->answer;
+    ndvi = input2->answer;
+    tempk = input3->answer;
+    rnet = input4->answer;
+    time = input5->answer;
+    result = output1->answer;
+    roerink = flag1->answer;
+
+    if ((infd_albedo = G_open_cell_old(albedo, "")) < 0)
+	G_fatal_error(_("Unable to open raster map <%s>"), albedo);
+    inrast_albedo = G_allocate_d_raster_buf();
+
+    if ((infd_ndvi = G_open_cell_old(ndvi, "")) < 0)
+	G_fatal_error(_("Unable to open raster map <%s>"), ndvi);
+    inrast_ndvi = G_allocate_d_raster_buf();
+
+    if ((infd_tempk = G_open_cell_old(tempk, "")) < 0)
+	G_fatal_error(_("Unable to open raster map <%s>"), tempk);
+    inrast_tempk = G_allocate_d_raster_buf();
+
+    if ((infd_rnet = G_open_cell_old(rnet, "")) < 0)
+	G_fatal_error(_("Unable to open raster map [%s]"), rnet);
+    inrast_rnet = G_allocate_d_raster_buf();
+
+    if ((infd_time = G_open_cell_old(time, "")) < 0)
+	G_fatal_error(_("Unable to open raster map [%s]"), time);
+    inrast_time = G_allocate_d_raster_buf();
+
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+    outrast = G_allocate_d_raster_buf();
+
+    /* Create New raster files */ 
+    if ((outfd = G_open_raster_new(result,DCELL_TYPE)) < 0)
+	G_fatal_error(_("Unable to open raster map<%s>"), result);
+    /* Process pixels */ 
+    for (row = 0; row < nrows; row++)
+    {
+	DCELL d;
+	DCELL d_albedo;
+	DCELL d_ndvi;
+	DCELL d_tempk;
+	DCELL d_rnet;
+	DCELL d_time;
+	G_percent(row, nrows, 2);	
+        /* read input maps */ 
+        if (G_get_d_raster_row(infd_albedo, inrast_albedo, row) < 0)
+	    G_fatal_error(_("Unable to read from <%s>"), albedo);
+	if (G_get_d_raster_row(infd_ndvi, inrast_ndvi, row)<0)
+	    G_fatal_error(_("Unable to read from <%s>"), ndvi);
+	if (G_get_d_raster_row(infd_tempk, inrast_tempk, row)< 0)
+	    G_fatal_error(_("Unable to read from <%s>"), tempk);
+	if (G_get_d_raster_row(infd_rnet, inrast_rnet, row)<0)
+	    G_fatal_error(_("Unable to read from <%s>"), rnet);
+	if (G_get_d_raster_row(infd_time, inrast_time, row)<0)
+	    G_fatal_error(_("Unable to read from <%s>"), time);
+        /*process the data */ 
+        for (col = 0; col < ncols; col++)
+        {
+            d_albedo = ((DCELL *) inrast_albedo)[col];
+            d_ndvi   = ((DCELL *) inrast_ndvi)[col];
+            d_tempk  = ((DCELL *) inrast_tempk)[col];
+            d_rnet   = ((DCELL *) inrast_rnet)[col];
+            d_time   = ((DCELL *) inrast_time)[col];
+	    if (G_is_d_null_value(&d_albedo) || 
+                 G_is_d_null_value(&d_ndvi) ||
+                 G_is_d_null_value(&d_tempk) ||
+		 G_is_d_null_value(&d_rnet) || 
+                 G_is_d_null_value(&d_time)) {
+		G_set_d_null_value(&outrast[col], 1);
+	    }
+	    else {
+                /* calculate soil heat flux         */ 
+                d=g_0(d_albedo,d_ndvi,d_tempk,d_rnet,d_time,roerink);
+		outrast[col] = d;
+	    }
+	}
+	if (G_put_d_raster_row(outfd, outrast) < 0)
+	    G_fatal_error(_("Unable to write to output raster map"));
+    }
+    G_free(inrast_albedo);
+    G_free(inrast_ndvi);
+    G_free(inrast_tempk);
+    G_free(inrast_rnet);
+    G_free(inrast_time);
+    G_close_cell(infd_albedo);
+    G_close_cell(infd_ndvi);
+    G_close_cell(infd_tempk);
+    G_close_cell(infd_rnet);
+    G_close_cell(infd_time);
+    G_free(outrast);
+    G_close_cell(outfd);
+
+    /* Colors in grey shade */ 
+    G_init_colors(&colors);
+    G_add_color_rule(0.0, 0, 0, 0, 200.0, 255, 255, 255, &colors);
+    G_short_history(result, "raster", &history);
+    G_command_history(&history);
+    G_write_history(result, &history);
+    exit(EXIT_SUCCESS);
+}
+