[GRASS-SVN] r37959 - in grass/trunk/imagery: . i.evapo.time_integration

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Jun 19 10:50:54 EDT 2009


Author: ychemin
Date: 2009-06-19 10:50:54 -0400 (Fri, 19 Jun 2009)
New Revision: 37959

Added:
   grass/trunk/imagery/i.evapo.time_integration/
   grass/trunk/imagery/i.evapo.time_integration/Makefile
   grass/trunk/imagery/i.evapo.time_integration/i.evapo.time_integration.html
   grass/trunk/imagery/i.evapo.time_integration/main.c
Modified:
   grass/trunk/imagery/Makefile
Log:
Added module for temporal integration of ET maps

Modified: grass/trunk/imagery/Makefile
===================================================================
--- grass/trunk/imagery/Makefile	2009-06-19 12:48:44 UTC (rev 37958)
+++ grass/trunk/imagery/Makefile	2009-06-19 14:50:54 UTC (rev 37959)
@@ -8,6 +8,7 @@
 	i.eb.evapfr \
 	i.eb.g0 \
 	i.eb.h_SEBAL01 \
+	i.evapo.time_integration \
 	i.emissivity \
 	i.find \
 	i.gensig \

Added: grass/trunk/imagery/i.evapo.time_integration/Makefile
===================================================================
--- grass/trunk/imagery/i.evapo.time_integration/Makefile	                        (rev 0)
+++ grass/trunk/imagery/i.evapo.time_integration/Makefile	2009-06-19 14:50:54 UTC (rev 37959)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.evapo.time_integration
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

Added: grass/trunk/imagery/i.evapo.time_integration/i.evapo.time_integration.html
===================================================================
--- grass/trunk/imagery/i.evapo.time_integration/i.evapo.time_integration.html	                        (rev 0)
+++ grass/trunk/imagery/i.evapo.time_integration/i.evapo.time_integration.html	2009-06-19 14:50:54 UTC (rev 37959)
@@ -0,0 +1,58 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.evapo.time_integration</EM> integrates ETa in time following a reference ET (typically) from a set of meteorological stations dataset.
+
+Inputs:
+- ETa images
+- ETa images DOY (Day of Year)
+- ETo images
+- ETo DOYmin as a single value 
+
+Method:
+1 - each ETa pixel is divided by the same day ETo and become ETrF
+2 - each ETrF pixel is multiplied by the ETo sum for the representative days
+3 - Sum all n temporal [ETrF*ETo_sum] pixels to make a summed(ET) in [DOYmin;DOYmax]
+
+representative days calculation:
+let assume i belongs to range [DOYmin;DOYmax]
+
+DOYbeforeETa[i] = ( DOYofETa[i] - DOYofETa[i-1] ) / 2
+DOYafterETa[i] = ( DOYofETa[i+1] - DOYofETa[i] ) / 2
+
+
+<H2>NOTES</H2>
+
+ETo images preparation:
+If you only have one meteorological station data, the easiest way is:
+
+n=0
+for ETo_val in Eto[1] Eto[2] ...
+do
+	r.mapcalc eto$n=$ETo_val 
+	`expr n = n + 1'
+done
+
+with Eto[1], Eto[2], etc being a simple copy and paste from your data file of all ETo values separated by an empty space from each other.
+
+If you have several meteorological stations data, then you need to grid them, Thiessen polygons or interpolation for each day.
+
+For multi-year calculations, just continue incrementing DOY values above 366, it will continue working, up to maximum input of 400 satellite images.
+
+<H2>TODO</H2>
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.eb.eta.html">i.eb.eta</A><br>
+<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
+<A HREF="i.evapo.SENAY.html">i.evapo.SENAY</A><br>
+<A HREF="r.surf.idw.html">r.surf.idw</A><br>
+<A HREF="r.surf.idw2.html">r.surf.idw2</A><br>
+<A HREF="r.bilinear.html">r.bilinear</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+Yann Chemin, International Rice Research Institute, The Philippines<BR>
+<p>
+<i>Last changed: $Date: 2008/08/25 06:17:43 $</i>

Added: grass/trunk/imagery/i.evapo.time_integration/main.c
===================================================================
--- grass/trunk/imagery/i.evapo.time_integration/main.c	                        (rev 0)
+++ grass/trunk/imagery/i.evapo.time_integration/main.c	2009-06-19 14:50:54 UTC (rev 37959)
@@ -0,0 +1,436 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.evapo.time_integration
+ * AUTHOR(S):    Yann Chemin - yann.chemin at gmail.com
+ * 		 Ines Cherif - icherif at yahoo.com
+ * PURPOSE:      Integrate in time the evapotranspiration from satellite,
+ *		 following a daily pattern from meteorological ETo.
+ *
+ * COPYRIGHT:    (C) 2008-09 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU Lesser General Public
+ *   	    	 License. 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>
+
+#define MAXFILES 400
+
+int main(int argc, char *argv[])
+{
+    struct Cell_head cellhd;	/*region+header info */
+    char *mapset;		/*mapset name */
+    int nrows, ncols;
+    int row, col;
+    struct GModule *module;
+    struct Option *input, *input1, *input2, *input3, *input4, *input5, *output;
+    struct History history;	/*metadata */
+    struct Colors colors;	/*Color rules */
+
+    /************************************/
+    char *name, *name1, *name2;	/*input raster name */
+    char *result;		/*output raster name */
+
+    /*File Descriptors */
+    int nfiles, nfiles1, nfiles2;
+    int infd[MAXFILES], infd1[MAXFILES], infd2[MAXFILES];
+    int outfd;
+
+    /****************************************/
+    /* Pointers for file names              */
+    char **names;
+    char **ptr;
+    char **names1;
+    char **ptr1;
+    char **names2;
+    char **ptr2;
+
+    /****************************************/
+    int DOYbeforeETa[MAXFILES], DOYafterETa[MAXFILES];
+    int bfr, aft;
+
+    /****************************************/
+    int ok;
+    int i = 0, j = 0;
+    double etodoy;		/*minimum ETo DOY */
+    double startperiod, endperiod;  /*first and last days (DOYs) of the period studied */
+    void *inrast[MAXFILES], *inrast1[MAXFILES], *inrast2[MAXFILES];
+    DCELL *outrast;
+    int data_format;		/* 0=double  1=float  2=32bit signed int  5=8bit unsigned int (ie text) */
+
+    RASTER_MAP_TYPE in_data_type[MAXFILES];	/* ETa */
+    RASTER_MAP_TYPE in_data_type1[MAXFILES];	/* DOY of ETa */
+    RASTER_MAP_TYPE in_data_type2[MAXFILES];	/* ETo */
+    RASTER_MAP_TYPE out_data_type = DCELL_TYPE;
+
+    /************************************/
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->keywords = _("evapotranspiration,temporal,integration");
+    module->description =_("Temporal integration of satellite ET actual (ETa) following the daily ET reference (ETo) from meteorological station(s)\n");
+
+    /* Define the different options */
+    input = G_define_standard_option(G_OPT_R_INPUTS);
+    input->key = _("eta");
+    input->description = _("Names of satellite ETa layers [mm/d or cm/d]");
+
+    input1 = G_define_standard_option(G_OPT_R_INPUTS);
+    input1->key = _("eta_doy");
+    input1->description =
+	_("Names of satellite ETa Day of Year (DOY) layers [0-400] [-]");
+
+    input2 = G_define_standard_option(G_OPT_R_INPUTS);
+    input2->key = _("eto");
+    input2->description =
+	_("Names of meteorological station ETo layers [0-400] [mm/d or cm/d]");
+
+    input3 = G_define_option();
+    input3->key = _("eto_doy_min");
+    input3->type = TYPE_DOUBLE;
+    input3->required = YES;
+    input3->description = _("Value of DOY for ETo first day");
+    
+    input4 = G_define_option();
+    input4->key = _("start_period");
+    input4->type = TYPE_DOUBLE;
+    input4->required = YES;
+    input4->description = _("Value of DOY for the first day of the period studied ");
+
+    input5 = G_define_option();
+    input5->key = _("end_period");
+    input5->type = TYPE_DOUBLE;
+    input5->required = YES;
+    input5->description = _("Value of DOY for the last day of the period studied ");
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+    output->description =
+	_("Name of the output temporally integrated ETa layer");
+
+    /* init nfiles */
+    nfiles = 1;
+    nfiles1 = 1;
+    nfiles2 = 1;
+
+    /********************/
+
+    if (G_parser(argc, argv))
+	exit(-1);
+
+    ok = 1;
+    names = input->answers;
+    ptr = input->answers;
+    names1 = input1->answers;
+    ptr1 = input1->answers;
+    names2 = input2->answers;
+    ptr2 = input2->answers;
+    etodoy = atof(input3->answer);
+    startperiod = atof(input4->answer);
+    endperiod = atof(input5->answer);
+    result = output->answer;
+
+    /****************************************/
+    if (endperiod<startperiod) {
+	G_fatal_error(_("The DOY for end_period can not be smaller than start_period"));
+	ok = 0;
+    }
+    if (etodoy>startperiod) {
+	G_fatal_error(_("The DOY for start_period can not be smaller than eto_doy_min"));
+	ok = 0;
+    }
+    if (G_legal_filename(result) < 0) {
+	G_fatal_error(_("<%s> is an illegal name"), result);
+	ok = 0;
+    }
+    for (; *ptr != NULL; ptr++) {
+	if (nfiles > MAXFILES)
+	    G_fatal_error(_("%s - too many ETa files. Only %d allowed"),
+			  G_program_name(), MAXFILES);
+	name = *ptr;
+	/* find map in mapset */
+	mapset = G_find_cell2(name, "");
+	if (mapset == NULL) {
+	    G_fatal_error(_("raster map <%s> not found"), name);
+	    ok = 0;
+	}
+	if (!ok)
+	    continue;
+	infd[nfiles] = G_open_cell_old(name, mapset);
+	if (infd[nfiles] < 0) {
+	    ok = 0;
+	    continue;
+	}
+	/* Allocate input buffer */
+	in_data_type[nfiles] = G_raster_map_type(name, mapset);
+	if ((infd[nfiles] = G_open_cell_old(name, mapset)) < 0)
+	    G_fatal_error(_("Unable to open raster map <%s>"), name);
+	if ((G_get_cellhd(name, mapset, &cellhd)) < 0)
+	    G_fatal_error(_("Unable to read file header of <%s>"), name);
+	inrast[nfiles] = G_allocate_raster_buf(in_data_type[nfiles]);
+	nfiles++;
+    }
+    nfiles--;
+    if (nfiles <= 1)
+	G_fatal_error(_("The min specified input map is two"));
+	
+	/****************************************/
+    ok = 1;
+    for (; *ptr1 != NULL; ptr1++) {
+	if (nfiles1 > MAXFILES)
+	    G_fatal_error(_("%s - too many ETa files. Only %d allowed"),
+			  G_program_name(), MAXFILES);
+	name1 = *ptr1;
+	/* find map in mapset */
+	mapset = G_find_cell2(name1, "");
+	if (mapset == NULL) {
+	    G_fatal_error(_("raster map <%s> not found"), name1);
+	    ok = 0;
+	}
+	if (!ok)
+	    continue;
+	infd1[nfiles1] = G_open_cell_old(name1, mapset);
+	if (infd1[nfiles1] < 0) {
+	    ok = 0;
+	    continue;
+	}
+	/* Allocate input buffer */
+	in_data_type1[nfiles1] = G_raster_map_type(name1, mapset);
+	if ((infd1[nfiles1] = G_open_cell_old(name1, mapset)) < 0)
+	    G_fatal_error(_("Unable to open raster map <%s>"), name1);
+	if ((G_get_cellhd(name1, mapset, &cellhd)) < 0)
+	    G_fatal_error(_("Unable to read file header of <%s>"), name1);
+	inrast1[nfiles1] = G_allocate_raster_buf(in_data_type1[nfiles1]);
+	nfiles1++;
+    }
+    nfiles1--;
+    if (nfiles1 <= 1)
+	G_fatal_error(_("The min specified input map is two"));
+
+
+	/****************************************/
+    if (nfiles != nfiles1)
+	G_fatal_error(_("ETa and ETa_DOY file numbers are not equal!"));
+
+	/****************************************/
+
+    ok = 1;
+    for (; *ptr2 != NULL; ptr2++) {
+	if (nfiles > MAXFILES)
+	    G_fatal_error(_("%s - too many ETa files. Only %d allowed"),
+			  G_program_name(), MAXFILES);
+	name2 = *ptr2;
+	/* find map in mapset */
+	mapset = G_find_cell2(name2, "");
+	if (mapset == NULL) {
+	    G_fatal_error(_("raster map <%s> not found"), name2);
+	    ok = 0;
+	}
+	if (!ok)
+	    continue;
+	infd2[nfiles2] = G_open_cell_old(name2, mapset);
+	if (infd2[nfiles2] < 0) {
+	    ok = 0;
+	    continue;
+	}
+	/* Allocate input buffer */
+	in_data_type2[nfiles2] = G_raster_map_type(name2, mapset);
+	if ((infd2[nfiles2] = G_open_cell_old(name2, mapset)) < 0) {
+	    G_fatal_error(_("Unable to open raster map <%s>"), name2);
+	}
+	if ((G_get_cellhd(name2, mapset, &cellhd)) < 0) {
+	    G_fatal_error(_("Unable to read file header of <%s>"), name2);
+	}
+	inrast2[nfiles2] = G_allocate_d_raster_buf();
+	nfiles2++;
+    }
+    nfiles2--;
+    if (nfiles2 <= 1)
+	G_fatal_error(_("The min specified input map is two"));
+
+    /* Allocate output buffer, use input map data_type */
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+    outrast = G_allocate_raster_buf(out_data_type);
+
+   
+    /* Create New raster files */
+    if ((outfd = G_open_raster_new(result, 1)) < 0)
+	G_fatal_error(_("Could not open <%s>"), result);
+
+    /*******************/
+    /* Process pixels */
+    double doy[MAXFILES];
+    double sum[MAXFILES];
+
+    for (row = 0; row < nrows; row++) 
+    {
+	DCELL d_out;
+	DCELL d_ETrF[MAXFILES];
+	DCELL d[MAXFILES];
+	DCELL d1[MAXFILES];
+	DCELL d2[MAXFILES];
+	G_percent(row, nrows, 2);
+
+	/* read input map */
+	for (i = 1; i <= nfiles; i++) 
+	    if ((G_get_d_raster_row(infd[i], inrast[i], row)) <	0) 
+		G_fatal_error(_("Could not read from <%s>"), name);
+	
+	for (i = 1; i <= nfiles1; i++) 
+	    if ((G_get_d_raster_row(infd1[i], inrast1[i], row)) < 0) 
+		G_fatal_error(_("Could not read from <%s>"), name1);
+
+	for (i = 1; i <= nfiles2; i++) 
+	    if ((G_get_d_raster_row (infd2[i], inrast2[i], row)) < 0) 
+		G_fatal_error(_("Could not read from <%s>"), name2);
+
+	/*process the data */
+	for (col = 0; col < ncols; col++) 
+        {
+            int	d1_null=0;
+            int	d_null=0;
+	    for (i = 1; i <= nfiles; i++) 
+            {
+		    if (G_is_d_null_value(&((DCELL *) inrast[i])[col]))
+		    	d_null=1;
+		    else
+	                d[i] = (double)((DCELL *) inrast[i])[col];
+	    }
+	    for (i = 1; i <= nfiles1; i++) 
+            {
+		    if (G_is_d_null_value(&((DCELL *) inrast1[i])[col]))
+			d1_null=1;
+		    else
+	                d1[i] = ((DCELL *) inrast1[i])[col];
+	    }
+
+	    for (i = 1; i <= nfiles2; i++) 
+		    d2[i] = ((DCELL *) inrast2[i])[col];
+
+	    /* Find out the DOY of the eto image    */
+	    for (i = 1; i <= nfiles1; i++) 
+            {
+		if ( d_null==1 || d1_null==1 )
+			G_set_d_null_value(&outrast[col],1);	
+		else
+		{
+			doy[i] = d1[i] - etodoy+1;
+			if (G_is_d_null_value(&d2[(int)doy[i]]) || d2[(int)doy[i]]==0 )
+				G_set_d_null_value(&outrast[col],1);
+			else
+				d_ETrF[i] = d[i] / d2[(int)doy[i]];
+		} 
+	    }
+
+	    for (i = 1; i <= nfiles1; i++) 
+            {
+		/* do nothing	*/
+		if ( d_null==1 || d1_null==1)
+                {
+			/*G_message("  null value ");*/
+                }
+		else
+		{
+			DOYbeforeETa[i]=0; DOYafterETa[i]=0;
+			if (i == 1)   
+				DOYbeforeETa[i] = startperiod;
+			else
+			{
+ 				int k=i-1;
+				while (d1[k]>=startperiod )
+				{
+					if (d1[k]<0)	 // case were d1[k] is null
+						k=k-1;					
+					else
+					{
+						DOYbeforeETa[i] = 1+((d1[i] + d1[k])/2.0);
+						break;
+					}			
+				}
+
+			}
+	
+			if (i == nfiles1)  
+				DOYafterETa[i] = endperiod;
+			else
+			{
+				int k=i+1;
+				while (d1[k]<=endperiod)
+				{
+					if (d1[k]<0)   // case were d1[k] is null
+						k=k+1;
+					else
+					{
+						DOYafterETa[i] = (d1[i] + d1[k]) / 2.0;
+						break;
+	   				}					
+				}
+			}
+		}	
+	    }
+
+	    sum[MAXFILES] = 0.0;
+	    for (i = 1; i <= nfiles1; i++) 
+            {
+		if(d_null==1 || d1_null==1)
+                {
+		    /* do nothing	 */
+		} 
+                else
+                {
+			if (DOYbeforeETa[i]==0 || DOYbeforeETa[i]==0 ) 	
+                            G_set_d_null_value(&outrast[col],1);
+			else 
+                        {
+				bfr = (int)DOYbeforeETa[i];
+				aft = (int)DOYafterETa[i];
+				sum[i]=0.0;
+				for (j = bfr; j < aft; j++) 
+					sum[i] += d2[(int)(j-etodoy+1)];
+			}
+		}
+	    }
+	
+	    d_out = 0.0;
+	    for (i = 1; i <= nfiles1; i++)
+            {
+		if(d_null==1 || d_null==1)
+			G_set_d_null_value(&outrast[col],1);
+		else
+                {	
+			d_out += d_ETrF[i] * sum[i];
+		     	outrast[col] = d_out;
+		}	
+	    }
+	}
+	if (G_put_raster_row(outfd, outrast, out_data_type) < 0)
+	    G_fatal_error(_("Unable to write to raster map<%s>"), result);
+    }
+
+    for (i = 1; i <= nfiles; i++) {
+	G_free(inrast[i]);
+	G_close_cell(infd[i]);
+	G_free(inrast1[i]);
+	G_close_cell(infd1[i]);
+	G_free(inrast2[i]);
+	G_close_cell(infd2[i]);
+    }
+    G_free(outrast);
+    G_close_cell(outfd);
+
+    /* Color table from 0.0 to 10.0 */
+    G_init_colors(&colors);
+    G_add_color_rule(0.0, 0, 0, 0, 10.0, 255, 255, 255, &colors);
+    /* Metadata */
+    G_short_history(result, "raster", &history);
+    G_command_history(&history);
+    G_write_history(result, &history);
+
+    exit(EXIT_SUCCESS);
+}



More information about the grass-commit mailing list