[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