[GRASS-SVN] r32882 - in grass-addons/gipe: . i.albedo
i.evapo.time_integration
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Aug 19 06:43:14 EDT 2008
Author: ychemin
Date: 2008-08-19 06:43:10 -0400 (Tue, 19 Aug 2008)
New Revision: 32882
Added:
grass-addons/gipe/i.evapo.time_integration/
grass-addons/gipe/i.evapo.time_integration/Makefile
grass-addons/gipe/i.evapo.time_integration/description.html
grass-addons/gipe/i.evapo.time_integration/main.c
Modified:
grass-addons/gipe/Makefile.imagery
grass-addons/gipe/gmmenu.tcl
grass-addons/gipe/i.albedo/description.html
grass-addons/gipe/i.albedo/main.c
grass-addons/gipe/menudata.py
Log:
added i.evapo.time_integration
Modified: grass-addons/gipe/Makefile.imagery
===================================================================
--- grass-addons/gipe/Makefile.imagery 2008-08-19 10:14:07 UTC (rev 32881)
+++ grass-addons/gipe/Makefile.imagery 2008-08-19 10:43:10 UTC (rev 32882)
@@ -34,6 +34,7 @@
i.evapo.potrad \
i.evapo.PT \
i.evapo.SENAY \
+ i.evapo.time_integration \
i.evapo.TSA \
i.qc.modis \
i.latitude \
Modified: grass-addons/gipe/gmmenu.tcl
===================================================================
--- grass-addons/gipe/gmmenu.tcl 2008-08-19 10:14:07 UTC (rev 32881)
+++ grass-addons/gipe/gmmenu.tcl 2008-08-19 10:43:10 UTC (rev 32882)
@@ -600,6 +600,8 @@
{command {[G_msg "Actual ET (SEBAL)"]} {} "i.eb.eta" {} -command {execute i.eb.eta }}
{command {[G_msg "Actual ET (SENAY)"]} {} "i.evapo.SENAY" {} -command {execute i.evapo.SENAY }}
{command {[G_msg "Actual ET (TSA)"]} {} "i.evapo.TSA" {} -command {execute i.evapo.TSA }}
+ {separator}
+ {command {[G_msg "Actual ET Temporal Integration"]} {} "i.evapo.time_integration" {} -command {execute i.evapo.time_integration }}
}}
{cascad {[G_msg "Energy Balance"]} {} "" $tmenu {
{command {[G_msg "Surface roughness (generic)"]} {} "i.eb.z0m0" {} -command {execute i.eb.z0m0 }}
Modified: grass-addons/gipe/i.albedo/description.html
===================================================================
--- grass-addons/gipe/i.albedo/description.html 2008-08-19 10:14:07 UTC (rev 32881)
+++ grass-addons/gipe/i.albedo/description.html 2008-08-19 10:43:10 UTC (rev 32882)
@@ -17,7 +17,7 @@
<H2>AUTHORS</H2>
-Yann Chemin, AIT, Thailand<BR>
+Yann Chemin, International Rice Research Institute, The Philippines<BR>
<p>
Modified: grass-addons/gipe/i.albedo/main.c
===================================================================
--- grass-addons/gipe/i.albedo/main.c 2008-08-19 10:14:07 UTC (rev 32881)
+++ grass-addons/gipe/i.albedo/main.c 2008-08-19 10:43:10 UTC (rev 32882)
@@ -134,6 +134,13 @@
landsat = (flag3->answer);
aster = (flag4->answer);
+ /***************************************************/
+ if (G_legal_filename (result) < 0)
+ {
+ G_fatal_error (_("[%s] is an illegal name"), result);
+ }
+ /***************************************************/
+ /***************************************************/
for (; *ptr != NULL; ptr++)
{
if (nfiles >= MAXFILES)
@@ -146,11 +153,6 @@
G_fatal_error (_("cell file [%s] not found"), name);
ok = 0;
}
- if (G_legal_filename (result) < 0)
- {
- G_fatal_error (_("[%s] is an illegal name"), result);
- ok = 0;
- }
if (!ok){
continue;
}
Added: grass-addons/gipe/i.evapo.time_integration/Makefile
===================================================================
--- grass-addons/gipe/i.evapo.time_integration/Makefile (rev 0)
+++ grass-addons/gipe/i.evapo.time_integration/Makefile 2008-08-19 10:43:10 UTC (rev 32882)
@@ -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-addons/gipe/i.evapo.time_integration/description.html
===================================================================
--- grass-addons/gipe/i.evapo.time_integration/description.html (rev 0)
+++ grass-addons/gipe/i.evapo.time_integration/description.html 2008-08-19 10:43:10 UTC (rev 32882)
@@ -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-addons/gipe/i.evapo.time_integration/main.c
===================================================================
--- grass-addons/gipe/i.evapo.time_integration/main.c (rev 0)
+++ grass-addons/gipe/i.evapo.time_integration/main.c 2008-08-19 10:43:10 UTC (rev 32882)
@@ -0,0 +1,398 @@
+/****************************************************************************
+ *
+ * MODULE: i.evapo.time_integration
+ * AUTHOR(S): Yann Chemin - ychemin at gmail.com
+ * PURPOSE: Integrate in time the evapotranspiration from satellite,
+ * following a daily pattern from meteorological ETo.
+ *
+ * COPYRIGHT: (C) 2008 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, *output;
+
+ struct History history; //metadata
+ struct Colors colors; //Color rules
+ /************************************/
+ /* FMEO Declarations*****************/
+ 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;
+ /****************************************/
+ double DOYbeforeETa[MAXFILES],DOYafterETa[MAXFILES];
+ int bfr,aft;
+ /****************************************/
+
+ int ok;
+ int i=0,j=0;
+ double etodoy; /*minimum ETo DOY*/
+
+ 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_INPUT) ;
+ input->key =_("eta");
+ input->multiple = YES;
+ input->description= _("Names of satellite ETa layers [mm/d or cm/d]");
+
+ input1 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input1->key =_("eta_doy");
+ input1->multiple = YES;
+ input1->description= _("Names of satellite ETa Day of Year (DOY) layers [0-400] [-]");
+
+ input2 = G_define_standard_option(G_OPT_R_INPUT) ;
+ input2->key =_("eto");
+ input2->multiple = YES;
+ 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");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT) ;
+ output->description= _("Name of the output temporally integrated ETa layer");
+
+ /* FMEO 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);
+
+ result = output->answer;
+ /****************************************/
+ 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 (_("cell file [%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(_("Cannot open cell file [%s]"), name);
+ }
+ if( (G_get_cellhd(name,mapset,&cellhd)) < 0){
+ G_fatal_error(_("Cannot 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 (_("cell file [%s] not found"), name1);
+ ok = 0;
+ }
+ if (!ok){
+ continue;
+ }
+ infd1[nfiles] = G_open_cell_old (name1, mapset);
+ if (infd1[nfiles] < 0)
+ {
+ ok = 0;
+ continue;
+ }
+ /* Allocate input buffer */
+ in_data_type1[nfiles] = G_raster_map_type(name1, mapset);
+ if( (infd1[nfiles] = G_open_cell_old(name1,mapset)) < 0){
+ G_fatal_error(_("Cannot open cell file [%s]"), name1);
+ }
+ if( (G_get_cellhd(name1,mapset,&cellhd)) < 0){
+ G_fatal_error(_("Cannot read file header of [%s]"), name1);
+ }
+ inrast1[nfiles] = G_allocate_raster_buf(in_data_type[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 (_("cell file [%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(_("Cannot open cell file [%s]"), name2);
+ }
+ if( (G_get_cellhd(name2,mapset,&cellhd)) < 0){
+ G_fatal_error(_("Cannot read file header of [%s]"), name2);
+ }
+ inrast2[nfiles2] = G_allocate_raster_buf(in_data_type2[nfiles2]);
+ 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_raster_row(infd[i],inrast[i],row,in_data_type[i])) < 0){
+ G_fatal_error (_("Could not read from <%s>"),name);
+ }
+ }
+ for (i=1;i<=nfiles1;i++)
+ {
+ if( (G_get_raster_row(infd1[i],inrast1[i],row,in_data_type1[i])) < 0){
+ G_fatal_error (_("Could not read from <%s>"),name1);
+ }
+ }
+ for (i=1;i<=nfiles2;i++)
+ {
+ if( (G_get_raster_row(infd2[i],inrast2[i],row,in_data_type2[i])) < 0){
+ G_fatal_error (_("Could not read from <%s>"),name2);
+ }
+ }
+ /*process the data */
+ for (col=0; col < ncols; col++)
+ {
+ for(i=1;i<=nfiles;i++)
+ {
+ switch(in_data_type[i])
+ {
+ case CELL_TYPE:
+ d[i] = (double) ((CELL *) inrast[i])[col];
+ break;
+ case FCELL_TYPE:
+ d[i] = (double) ((FCELL *) inrast[i])[col];
+ break;
+ case DCELL_TYPE:
+ d[i] = (double) ((DCELL *) inrast[i])[col];
+ break;
+ }
+ }
+ for(i=1;i<=nfiles1;i++)
+ {
+ switch(in_data_type1[i])
+ {
+ case CELL_TYPE:
+ d1[i] = (double) ((CELL *) inrast1[i])[col];
+ break;
+ case FCELL_TYPE:
+ d1[i] = (double) ((FCELL *) inrast1[i])[col];
+ break;
+ case DCELL_TYPE:
+ d1[i] = (double) ((DCELL *) inrast1[i])[col];
+ break;
+ }
+ }
+ for(i=1;i<=nfiles2;i++)
+ {
+ switch(in_data_type2[i])
+ {
+ case CELL_TYPE:
+ d2[i] = (double) ((CELL *) inrast2[i])[col];
+ break;
+ case FCELL_TYPE:
+ d2[i] = (double) ((FCELL *) inrast2[i])[col];
+ break;
+ case DCELL_TYPE:
+ d2[i] = (double) ((DCELL *) inrast2[i])[col];
+ break;
+ }
+ }
+ /* Find out the DOY of the eto image */
+ /* DOY_eto_index = ModisDOY - etodoymin */
+ for(i=1;i<=nfiles1;i++)
+ {
+ doy[i] = d1[i] - etodoy;
+ d_ETrF[i]=d[i]/d2[(int)doy[i]];
+ }
+ for(i=1;i<=nfiles1;i++)
+ {
+ if(i==1)
+ DOYbeforeETa[i]=etodoy;
+ else
+ DOYbeforeETa[i]=(d[i]-d[i-1])/2.0;
+ if(i==nfiles1)
+ DOYafterETa[i]=etodoy+nfiles2;
+ else
+ DOYafterETa[i]=(d[i+1]-d[i])/2.0;
+ }
+ sum[MAXFILES]=0.0;
+ for(i=1;i<=nfiles1;i++)
+ {
+ bfr = (int) DOYbeforeETa[i];
+ aft = (int) DOYafterETa[i];
+ for(j=bfr;j<aft;j++)
+ {
+ sum[i]+=d2[j];
+ }
+ }
+ d_out = 0.0;
+ for(i=1;i<=nfiles1;i++)
+ {
+ 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 (_("Cannot write to <%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 1.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);
+}
Modified: grass-addons/gipe/menudata.py
===================================================================
--- grass-addons/gipe/menudata.py 2008-08-19 10:14:07 UTC (rev 32881)
+++ grass-addons/gipe/menudata.py 2008-08-19 10:43:10 UTC (rev 32882)
@@ -1510,6 +1510,11 @@
_("Actual ET (TSA) NOT WORKING YET"),
"self.RunMenuCmd",
"i.evapo.TSA"),
+ ("","","", ""),
+ (_("Actual ET Time Integration"),
+ _("Temporal Integration of Actual ET using daily ETo"),
+ "self.RunMenuCmd",
+ "i.evapo.time_integration"),
)
),
(_("Energy Balance"), (
More information about the grass-commit
mailing list