[GRASS-SVN] r45119 - grass-addons/imagery/gipe/i.evapo.senay
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jan 20 17:36:08 EST 2011
Author: ychemin
Date: 2011-01-20 14:36:08 -0800 (Thu, 20 Jan 2011)
New Revision: 45119
Added:
grass-addons/imagery/gipe/i.evapo.senay/i.evapo.senay.html
Removed:
grass-addons/imagery/gipe/i.evapo.senay/description.html
Modified:
grass-addons/imagery/gipe/i.evapo.senay/Makefile
grass-addons/imagery/gipe/i.evapo.senay/main.c
Log:
Converted i.evapo.senay to G7, added ET Potential map input option
Modified: grass-addons/imagery/gipe/i.evapo.senay/Makefile
===================================================================
--- grass-addons/imagery/gipe/i.evapo.senay/Makefile 2011-01-20 22:31:13 UTC (rev 45118)
+++ grass-addons/imagery/gipe/i.evapo.senay/Makefile 2011-01-20 22:36:08 UTC (rev 45119)
@@ -2,8 +2,8 @@
PGM = i.evapo.senay
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Deleted: grass-addons/imagery/gipe/i.evapo.senay/description.html
===================================================================
--- grass-addons/imagery/gipe/i.evapo.senay/description.html 2011-01-20 22:31:13 UTC (rev 45118)
+++ grass-addons/imagery/gipe/i.evapo.senay/description.html 2011-01-20 22:36:08 UTC (rev 45119)
@@ -1,39 +0,0 @@
-<H2>DESCRIPTION</H2>
-
-<EM>i.evapo.SENAY</EM> Calculates the diurnal actual evapotranspiration after Senay (2007). This is converting all Net radiation from the diurnal period into ET, then uses Senay equation for evaporative fraction.
-
-It takes input maps of Albedo, surface skin temperature, latitude, day of year, single-way transmissivity and takes input value of the density of fresh water.
-
-DEM is used for calculating min and max temperature for Senay equation.
-
-The "-s" flag permits output map of evaporative fraction from Senay.
-
-
-<H2>NOTES</H2>
-If you are trying to map irrigated crops, and you know there is at least one crop pixel that is evapotranspiring at maximum (ETa=ETpot), then read this.
-
-i.evapo.SENAY is highly dependent on the wet pixel being the lowest temperature in the crop pixels to work for non water stressed crops, force it that way, even if it breaks non crop areas. I suggest you reduce your region to the irrigation system boundaries, checking that it includes a bit of dry area for the hot/dry pixel.
-
-Since it is a direct relationship to LST, evaporative fraction can be very sensitive to the kind of pixel sample you feed it with.
-
-<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.eb.eta.html">i.eb.eta</A><br>
-<A HREF="i.eb.evapfr.html">i.eb.evapfr</A><br>
-<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
-</em>
-
-
-<H2>AUTHORS</H2>
-
-Yann Chemin, International Rice Research Institute, The Philippines.<BR>
-
-
-<p>
-<i>Last changed: $Date$</i>
Copied: grass-addons/imagery/gipe/i.evapo.senay/i.evapo.senay.html (from rev 45115, grass-addons/imagery/gipe/i.evapo.senay/description.html)
===================================================================
--- grass-addons/imagery/gipe/i.evapo.senay/i.evapo.senay.html (rev 0)
+++ grass-addons/imagery/gipe/i.evapo.senay/i.evapo.senay.html 2011-01-20 22:36:08 UTC (rev 45119)
@@ -0,0 +1,39 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.evapo.SENAY</EM> Calculates the diurnal actual evapotranspiration after Senay (2007). This is converting all Net radiation from the diurnal period into ET, then uses Senay equation for evaporative fraction.
+
+It takes input maps of Albedo, surface skin temperature, latitude, day of year, single-way transmissivity and takes input value of the density of fresh water.
+
+DEM is used for calculating min and max temperature for Senay equation.
+
+The "-s" flag permits output map of evaporative fraction from Senay.
+
+
+<H2>NOTES</H2>
+If you are trying to map irrigated crops, and you know there is at least one crop pixel that is evapotranspiring at maximum (ETa=ETpot), then read this.
+
+i.evapo.SENAY is highly dependent on the wet pixel being the lowest temperature in the crop pixels to work for non water stressed crops, force it that way, even if it breaks non crop areas. I suggest you reduce your region to the irrigation system boundaries, checking that it includes a bit of dry area for the hot/dry pixel.
+
+Since it is a direct relationship to LST, evaporative fraction can be very sensitive to the kind of pixel sample you feed it with.
+
+<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.eb.eta.html">i.eb.eta</A><br>
+<A HREF="i.eb.evapfr.html">i.eb.evapfr</A><br>
+<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, International Rice Research Institute, The Philippines.<BR>
+
+
+<p>
+<i>Last changed: $Date$</i>
Modified: grass-addons/imagery/gipe/i.evapo.senay/main.c
===================================================================
--- grass-addons/imagery/gipe/i.evapo.senay/main.c 2011-01-20 22:31:13 UTC (rev 45118)
+++ grass-addons/imagery/gipe/i.evapo.senay/main.c 2011-01-20 22:36:08 UTC (rev 45119)
@@ -6,7 +6,7 @@
* PURPOSE: Creates a map of actual evapotranspiration following
* the method of Senay (2007).
*
- * COPYRIGHT: (C) 2008 by the GRASS Development Team
+ * COPYRIGHT: (C) 2008-2011 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
@@ -19,464 +19,284 @@
#include <stdlib.h>
#include <string.h>
#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
-double solar_day(double lat, double doy, double tsw);
+double solar_day(double lat, double doy, double tsw);
double solar_day_3d(double lat, double doy, double tsw, double slope,
double aspect);
double r_net_day(double bbalb, double solar, double tsw);
-
double r_net_day_bandara98(double surface_albedo, double solar_day,
double apparent_atm_emissivity,
double surface_emissivity,
double air_temperature);
double et_pot_day(double rnetd, double tempk, double roh_w);
-
double evapfr_senay(double t_hot, double t_cold, double tempk);
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 *input1, *input2, *input3, *input4;
-
struct Option *input5, *input6, *input7, *input8, *input9;
-
- struct Option *input10, *input11, *input12, *input13;
-
+ struct Option *input10, *input11, *input12, *input13, *input14;
struct Option *output1, *output2;
-
- struct Flag *flag2, *flag3;
-
+ struct Flag *flag1, *flag2, *flag3;
struct History history; /*metadata */
-
struct Colors colors; /*Color rules */
+ CELL val1,val2; /*Color ranges */
-
-
- /************************************/
- /* FMEO Declarations**************** */
char *name; /*input raster name */
-
char *result1, *result2; /*output raster name */
-
-
- /*File Descriptors */
int infd_albedo, infd_tempk, infd_dem;
-
int infd_lat, infd_doy, infd_tsw;
-
int infd_slope, infd_aspect;
-
int infd_tair, infd_e0;
-
- int infd_ndvi;
-
+ int infd_ndvi, infd_etpotd;
int outfd1, outfd2;
-
char *albedo, *tempk, *dem, *lat, *doy, *tsw;
-
- char *slope, *aspect, *tair, *e0, *ndvi;
-
+ char *slope, *aspect, *tair, *e0, *ndvi, *etpotd;
double roh_w, e_atm;
-
int i = 0, j = 0;
-
void *inrast_albedo, *inrast_tempk;
-
void *inrast_dem, *inrast_lat;
-
void *inrast_doy, *inrast_tsw;
-
void *inrast_slope, *inrast_aspect;
-
void *inrast_tair, *inrast_e0;
+ void *inrast_ndvi, *inrast_etpotd;
- void *inrast_ndvi;
-
DCELL * outrast1, *outrast2;
- RASTER_MAP_TYPE data_type_output = DCELL_TYPE;
- RASTER_MAP_TYPE data_type_albedo;
- RASTER_MAP_TYPE data_type_tempk;
- RASTER_MAP_TYPE data_type_dem;
- RASTER_MAP_TYPE data_type_lat;
- RASTER_MAP_TYPE data_type_doy;
- RASTER_MAP_TYPE data_type_tsw;
- RASTER_MAP_TYPE data_type_slope;
- RASTER_MAP_TYPE data_type_aspect;
- RASTER_MAP_TYPE data_type_tair;
- RASTER_MAP_TYPE data_type_e0;
- RASTER_MAP_TYPE data_type_ndvi;
-
- /********************************/
- /* Stats for Senay equation */
+ /********************************/
+ /* Stats for Senay equation */
double t0dem_min = 400.0, t0dem_max = 200.0;
-
double tempk_min = 400.0, tempk_max = 200.0;
+ /********************************/
-
+ G_gisinit(argv[0]);
- /********************************/
- G_gisinit(argv[0]);
module = G_define_module();
- module->keywords = _("Actual ET, evapotranspiration, Senay");
+ G_add_keyword(_("Actual ET"));
+ G_add_keyword(_("Evapotranspiration"));
+ G_add_keyword(_("Senay"));
module->description =
_("Actual evapotranspiration, method after Senay (2007)");
- /* Define the different options */
- input1 = G_define_standard_option(G_OPT_R_INPUT);
- input1->key = _("albedo");
- input1->description = _("Name of the Albedo layer [0.0-1.0]");
+ /* Define the different options */
+ input1 = G_define_standard_option(G_OPT_R_INPUT);
+ input1->key = _("temperature");
+ input1->description = _("Name of the temperature layer [Degree Kelvin]");
+
input2 = G_define_standard_option(G_OPT_R_INPUT);
- input2->key = _("tempk");
- input2->description = _("Name of the temperature layer [Degree Kelvin]");
+ input2->key = _("albedo");
+ input2->description = _("Name of the Albedo layer [0.0-1.0]");
+
input3 = G_define_standard_option(G_OPT_R_INPUT);
- input3->key = _("dem");
+ input3->key = _("elevation");
input3->description = _("Name of the elevation layer [m]");
+
input4 = G_define_standard_option(G_OPT_R_INPUT);
- input4->key = _("lat");
+ input4->key = _("latitude");
+ input4->required = NO;
input4->description = _("Name of the degree latitude layer [dd.ddd]");
+ input4->guisection = _("Optional");
+
input5 = G_define_standard_option(G_OPT_R_INPUT);
- input5->key = _("doy");
+ input5->key = _("dayofyear");
+ input5->required = NO;
input5->description = _("Name of the Day of Year layer [0.0-366.0]");
+ input5->guisection = _("Optional");
+
input6 = G_define_standard_option(G_OPT_R_INPUT);
- input6->key = _("tsw");
+ input6->key = _("transmissivitysingleway");
+ input6->required = NO;
input6->description =
_("Name of the single-way transmissivity layer [0.0-1.0]");
+ input6->guisection = _("Optional");
+
input7 = G_define_option();
- input7->key = _("roh_w");
+ input7->key = _("waterdensity");
input7->type = TYPE_DOUBLE;
- input7->required = YES;
+ input7->required = NO;
input7->description =
_("Value of the density of fresh water ~[1000-1020]");
input7->answer = _("1005.0");
+ input7->guisection = _("Optional");
+
input8 = G_define_standard_option(G_OPT_R_INPUT);
input8->key = _("slope");
input8->required = NO;
input8->description = _("Name of the Slope layer ~[0-90]");
input8->guisection = _("Optional");
+
input9 = G_define_standard_option(G_OPT_R_INPUT);
input9->key = _("aspect");
input9->required = NO;
input9->description = _("Name of the Aspect layer ~[0-360]");
input9->guisection = _("Optional");
+
input10 = G_define_option();
- input10->key = _("e_atm");
+ input10->key = _("atmosphericemissivity");
input10->type = TYPE_DOUBLE;
input10->required = NO;
input10->description =
_("Value of the apparent atmospheric emissivity (Bandara, 1998 used 0.845 for Sri Lanka)");
input10->guisection = _("Optional");
+
input11 = G_define_standard_option(G_OPT_R_INPUT);
- input11->key = _("t_air");
+ input11->key = _("airtemperature");
input11->required = NO;
input11->description =
_("Name of the Air Temperature layer [Kelvin], use with -b");
input11->guisection = _("Optional");
+
input12 = G_define_standard_option(G_OPT_R_INPUT);
- input12->key = _("e0");
+ input12->key = _("surfaceemissivity");
input12->required = NO;
input12->description =
_("Name of the Surface Emissivity layer [-], use with -b");
input12->guisection = _("Optional");
+
input13 = G_define_standard_option(G_OPT_R_INPUT);
input13->key = _("ndvi");
input13->description = _("Name of the NDVI layer [-]");
+
+ input14 = G_define_standard_option(G_OPT_R_INPUT);
+ input14->key = _("diurnaletpotential");
+ input14->required = NO;
+ input14->description = _("Name of the ET Potential layer [mm/day]");
+ input14->guisection = _("Optional");
+
output1 = G_define_standard_option(G_OPT_R_OUTPUT);
output1->description = _("Name of the output Actual ET layer");
+
output2 = G_define_standard_option(G_OPT_R_OUTPUT);
output2->key = _("evapfr");
output2->required = NO;
output2->description =
_("Name of the output evaporative fraction layer");
output2->guisection = _("Optional");
+
+ flag1 = G_define_flag();
+ flag1->key = 'e';
+ flag1->description = _("ET Potential Map as input (By-Pass creation of one)");
+
flag2 = G_define_flag();
flag2->key = 'd';
flag2->description = _("Slope/Aspect correction");
+
flag3 = G_define_flag();
flag3->key = 'b';
flag3->description =
_("Net Radiation Bandara (1998), generic Longwave calculation, need apparent atmospheric emissivity, Air temperature and surface emissivity inputs");
-
- /********************/
- if (G_parser(argc, argv))
+ /********************/
+ if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- albedo = input1->answer;
- tempk = input2->answer;
+
+ tempk = input1->answer;
+ albedo = input2->answer;
+ ndvi = input13->answer;
dem = input3->answer;
- lat = input4->answer;
- doy = input5->answer;
- tsw = input6->answer;
- roh_w = atof(input7->answer);
- slope = input8->answer;
- aspect = input9->answer;
- if (input10->answer) {
- e_atm = atof(input10->answer);
+ if(flag1->answer)
+ etpotd = input14->answer;
+ else{
+ lat = input4->answer;
+ doy = input5->answer;
+ tsw = input6->answer;
+ roh_w = atof(input7->answer);
+ if (flag2->answer) {
+ slope = input8->answer;
+ aspect = input9->answer;
+ }
+ if (input10->answer) {
+ e_atm = atof(input10->answer);
+ }
+ if (flag3->answer) {
+ tair = input11->answer;
+ e0 = input12->answer;
+ }
}
- tair = input11->answer;
- e0 = input12->answer;
- ndvi = input13->answer;
result1 = output1->answer;
result2 = output2->answer;
+ /***************************************************/
+ infd_albedo = Rast_open_old(albedo, "");
+ inrast_albedo = Rast_allocate_d_buf();
+ infd_tempk = Rast_open_old(tempk, "");
+ inrast_tempk = Rast_allocate_d_buf();
+ infd_dem = Rast_open_old(dem, "");
+ inrast_dem = Rast_allocate_d_buf();
+ infd_ndvi = Rast_open_old(ndvi, "");
+ inrast_ndvi = Rast_allocate_d_buf();
- /***************************************************/
- mapset = G_find_cell2(albedo, "");
- if (mapset == NULL) {
- G_fatal_error(_("cell file [%s] not found"), albedo);
+ if (flag1->answer) {
+ infd_etpotd = Rast_open_old(etpotd, "");
+ inrast_etpotd = Rast_allocate_d_buf();
+ }
+ else {
+ infd_lat = Rast_open_old(lat, "");
+ inrast_lat = Rast_allocate_d_buf();
+ infd_doy = Rast_open_old(doy, "");
+ inrast_doy = Rast_allocate_d_buf();
+ infd_tsw = Rast_open_old(tsw, "");
+ inrast_tsw = Rast_allocate_d_buf();
+ if (flag2->answer) {
+ infd_slope = Rast_open_old(slope, "");
+ inrast_slope = Rast_allocate_d_buf();
+ infd_aspect = Rast_open_old(aspect, "");
+ inrast_aspect = Rast_allocate_d_buf();
+ }
+ if (flag3->answer) {
+ infd_tair = Rast_open_old(tair, "");
+ inrast_tair = Rast_allocate_d_buf();
+ infd_e0 = Rast_open_old(e0, "");
+ inrast_e0 = Rast_allocate_d_buf();
+ }
}
- data_type_albedo = G_raster_map_type(albedo, mapset);
- if ((infd_albedo = G_open_cell_old(albedo, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), albedo);
- if (G_get_cellhd(albedo, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s])"), albedo);
- inrast_albedo = G_allocate_raster_buf(data_type_albedo);
+ /***************************************************/
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ outrast1 = Rast_allocate_d_buf();
+ if (result2)
+ outrast2 = Rast_allocate_d_buf();
-
- /***************************************************/
- mapset = G_find_cell2(tempk, "");
- if (mapset == NULL) {
- G_fatal_error(_("cell file [%s] not found"), tempk);
- }
- data_type_tempk = G_raster_map_type(tempk, mapset);
- if ((infd_tempk = G_open_cell_old(tempk, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), tempk);
- if (G_get_cellhd(tempk, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), tempk);
- inrast_tempk = G_allocate_raster_buf(data_type_tempk);
+ /* Create New raster files */
+ outfd1 = Rast_open_new(result1, DCELL_TYPE);
+ if (result2)
+ outfd2 = Rast_open_new(result2, DCELL_TYPE);
-
- /***************************************************/
- mapset = G_find_cell2(dem, "");
- if (mapset == NULL) {
- G_fatal_error(_("cell file [%s] not found"), dem);
- }
- data_type_dem = G_raster_map_type(dem, mapset);
- if ((infd_dem = G_open_cell_old(dem, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), dem);
- if (G_get_cellhd(dem, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), dem);
- inrast_dem = G_allocate_raster_buf(data_type_dem);
-
-
- /***************************************************/
- mapset = G_find_cell2(lat, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), lat);
- }
- data_type_lat = G_raster_map_type(lat, mapset);
- if ((infd_lat = G_open_cell_old(lat, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), lat);
- if (G_get_cellhd(lat, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), lat);
- inrast_lat = G_allocate_raster_buf(data_type_lat);
-
-
- /***************************************************/
- mapset = G_find_cell2(doy, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), doy);
- }
- data_type_doy = G_raster_map_type(doy, mapset);
- if ((infd_doy = G_open_cell_old(doy, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), doy);
- if (G_get_cellhd(doy, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), doy);
- inrast_doy = G_allocate_raster_buf(data_type_doy);
-
-
- /***************************************************/
- mapset = G_find_cell2(tsw, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), tsw);
- }
- data_type_tsw = G_raster_map_type(tsw, mapset);
- if ((infd_tsw = G_open_cell_old(tsw, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), tsw);
- if (G_get_cellhd(tsw, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), tsw);
- inrast_tsw = G_allocate_raster_buf(data_type_tsw);
-
-
- /***************************************************/
- if (flag2->answer) {
- mapset = G_find_cell2(slope, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), slope);
- }
- data_type_slope = G_raster_map_type(slope, mapset);
- if ((infd_slope = G_open_cell_old(slope, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), slope);
- if (G_get_cellhd(slope, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), slope);
- inrast_slope = G_allocate_raster_buf(data_type_slope);
- }
-
-
- /***************************************************/
- if (flag2->answer) {
- mapset = G_find_cell2(aspect, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), aspect);
- }
- data_type_aspect = G_raster_map_type(aspect, mapset);
- if ((infd_aspect = G_open_cell_old(aspect, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), aspect);
- if (G_get_cellhd(aspect, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), aspect);
- inrast_aspect = G_allocate_raster_buf(data_type_aspect);
- }
-
-
- /***************************************************/
- if (flag3->answer) {
- mapset = G_find_cell2(tair, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), tair);
- }
- data_type_tair = G_raster_map_type(tair, mapset);
- if ((infd_tair = G_open_cell_old(tair, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), tair);
- if (G_get_cellhd(tair, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), tair);
- inrast_tair = G_allocate_raster_buf(data_type_tair);
- }
-
-
- /***************************************************/
- if (flag3->answer) {
- mapset = G_find_cell2(e0, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), e0);
- }
- data_type_e0 = G_raster_map_type(e0, mapset);
- if ((infd_e0 = G_open_cell_old(e0, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), e0);
- if (G_get_cellhd(e0, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), e0);
- inrast_e0 = G_allocate_raster_buf(data_type_e0);
- }
-
-
- /***************************************************/
- mapset = G_find_cell2(ndvi, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), ndvi);
- }
- data_type_ndvi = G_raster_map_type(ndvi, mapset);
- if ((infd_ndvi = G_open_cell_old(ndvi, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), ndvi);
- if (G_get_cellhd(ndvi, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), ndvi);
- inrast_ndvi = G_allocate_raster_buf(data_type_ndvi);
-
-
- /***************************************************/
- G_debug(3, "number of rows %d", cellhd.rows);
- nrows = G_window_rows();
- ncols = G_window_cols();
- outrast1 = G_allocate_raster_buf(data_type_output);
- if (result2) {
- outrast2 = G_allocate_raster_buf(data_type_output);
- }
-
- /* Create New raster files */
- if ((outfd1 = G_open_raster_new(result1, data_type_output)) < 0)
- G_fatal_error(_("Could not open <%s>"), result1);
- if (result2) {
- if ((outfd2 = G_open_raster_new(result2, data_type_output)) < 0)
- G_fatal_error(_("Could not open <%s>"), result2);
- }
-
- /* Process tempk min / max pixels for SENAY Evapfr */
- for (row = 0; row < nrows; row++)
- {
+ /* Process tempk min / max pixels for SENAY Evapfr */
+ for (row = 0; row < nrows; row++)
+ {
DCELL d;
- DCELL d_albedo;
- DCELL d_tempk;
- DCELL d_dem;
- DCELL d_t0dem;
- DCELL d_ndvi;
+ DCELL d_albedo,d_tempk;
+ DCELL d_dem,d_t0dem,d_ndvi;
+
G_percent(row, nrows, 2);
- if (G_get_raster_row
- (infd_albedo, inrast_albedo, row, data_type_albedo) < 0)
- G_fatal_error(_("Could not read from <%s>"), albedo);
- if (G_get_raster_row(infd_tempk, inrast_tempk, row, data_type_tempk)
- < 0)
- G_fatal_error(_("Could not read from <%s>"), tempk);
- if (G_get_raster_row(infd_dem, inrast_dem, row, data_type_dem) < 0)
- G_fatal_error(_("Could not read from <%s>"), dem);
- if (G_get_raster_row(infd_ndvi, inrast_ndvi, row, data_type_ndvi) <
- 0)
- G_fatal_error(_("Could not read from <%s>"), ndvi);
+
+ Rast_get_row(infd_albedo, inrast_albedo, row, DCELL_TYPE);
+ Rast_get_row(infd_tempk, inrast_tempk, row, DCELL_TYPE);
+ Rast_get_row(infd_dem, inrast_dem, row, DCELL_TYPE);
+ Rast_get_row(infd_ndvi, inrast_ndvi, row, DCELL_TYPE);
- /*process the data */
- for (col = 0; col < ncols; col++)
- {
- switch (data_type_albedo) {
- case CELL_TYPE:
- d_albedo = (double)((CELL *) inrast_albedo)[col];
- break;
- case FCELL_TYPE:
- d_albedo = (double)((FCELL *) inrast_albedo)[col];
- break;
- case DCELL_TYPE:
- d_albedo = (double)((DCELL *) inrast_albedo)[col];
- break;
- }
- switch (data_type_tempk) {
- case CELL_TYPE:
- d_tempk = (double)((CELL *) inrast_tempk)[col];
- break;
- case FCELL_TYPE:
- d_tempk = (double)((FCELL *) inrast_tempk)[col];
- break;
- case DCELL_TYPE:
- d_tempk = (double)((DCELL *) inrast_tempk)[col];
- break;
- }
- switch (data_type_dem) {
- case CELL_TYPE:
- d_dem = (double)((CELL *) inrast_dem)[col];
- break;
- case FCELL_TYPE:
- d_dem = (double)((FCELL *) inrast_dem)[col];
- break;
- case DCELL_TYPE:
- d_dem = (double)((DCELL *) inrast_dem)[col];
- break;
- }
- switch (data_type_ndvi) {
- case CELL_TYPE:
- d_ndvi = (double)((CELL *) inrast_ndvi)[col];
- break;
- case FCELL_TYPE:
- d_ndvi = (double)((FCELL *) inrast_ndvi)[col];
- break;
- case DCELL_TYPE:
- d_ndvi = (double)((DCELL *) inrast_ndvi)[col];
- break;
- }
- if (G_is_d_null_value(&d_albedo) || G_is_d_null_value(&d_tempk)
- || G_is_d_null_value(&d_dem) ||
- G_is_d_null_value(&d_ndvi)) {
-
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ d_albedo = (double)((DCELL *) inrast_albedo)[col];
+ d_tempk = (double)((DCELL *) inrast_tempk)[col];
+ d_dem = (double)((DCELL *) inrast_dem)[col];
+ d_ndvi = (double)((DCELL *) inrast_ndvi)[col];
+ if (Rast_is_d_null_value(&d_albedo) || Rast_is_d_null_value(&d_tempk)
+ || Rast_is_d_null_value(&d_dem) ||
+ Rast_is_d_null_value(&d_ndvi)) {
/* do nothing */
}
else {
d_t0dem = d_tempk + 0.00649 * d_dem;
- if (d_t0dem < 0.0 || d_albedo < 0.001) {
-
+ if (d_t0dem < 0.0 || d_albedo < 0.001){
/* do nothing */
}
else {
@@ -490,257 +310,179 @@
}
}
}
- }
}
- G_message("tempk_min=%f\ntempk_max=%f\n", tempk_min, tempk_max);
-
- /* Process pixels */
- for (row = 0; row < nrows; row++)
- {
- DCELL d;
- DCELL d_albedo;
- DCELL d_tempk;
- DCELL d_lat;
- DCELL d_doy;
- DCELL d_tsw;
- DCELL d_solar;
- DCELL d_rnetd;
- DCELL d_slope;
- DCELL d_aspect;
- DCELL d_tair;
- DCELL d_e0;
- DCELL d_ndvi;
- DCELL d_etpotd;
- DCELL d_evapfr;
- G_percent(row, nrows, 2);
+ G_message("tempk_min=%f\ntempk_max=%f\n", tempk_min, tempk_max);
+
+ /**********************************************************/
+ /* If ET Potential/Reference is an input compute ETa Senay*/
+ if(flag1->answer){
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d;
+ DCELL d_etpotd,d_evapfr;
+
+ G_percent(row, nrows, 2);
- /* read input maps */
- if (G_get_raster_row
- (infd_albedo, inrast_albedo, row, data_type_albedo) < 0)
- G_fatal_error(_("Could not read from <%s>"), albedo);
- if (G_get_raster_row(infd_tempk, inrast_tempk, row, data_type_tempk)
- < 0)
- G_fatal_error(_("Could not read from <%s>"), tempk);
- if (G_get_raster_row(infd_lat, inrast_lat, row, data_type_lat) < 0)
- G_fatal_error(_("Could not read from <%s>"), lat);
- if (G_get_raster_row(infd_doy, inrast_doy, row, data_type_doy) < 0)
- G_fatal_error(_("Could not read from <%s>"), doy);
- if (G_get_raster_row(infd_tsw, inrast_tsw, row, data_type_tsw) < 0)
- G_fatal_error(_("Could not read from <%s>"), tsw);
- if (flag2->answer) {
- if (G_get_raster_row
- (infd_slope, inrast_slope, row, data_type_slope) < 0)
- G_fatal_error(_("Could not read from <%s>"), slope);
- if (G_get_raster_row
- (infd_aspect, inrast_aspect, row, data_type_aspect) < 0)
- G_fatal_error(_("Could not read from <%s>"), aspect);
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ d_evapfr = evapfr_senay(tempk_max, tempk_min, d_tempk);
+ /*If water then no water stress */
+ if (d_albedo <= 0.1 && d_ndvi <= 0.0)
+ d_evapfr = 1.0;
+ /*some points are colder than selected low */
+ if (d_evapfr > 1.0)
+ d_evapfr = 1.0;
+ if (result2)
+ outrast2[col] = d_evapfr;
+ d = d_etpotd * d_evapfr;
+ outrast1[col] = d;
+ }
+ Rast_put_row(outfd1, outrast1, DCELL_TYPE);
+ if (result2)
+ Rast_put_row(outfd2, outrast2, DCELL_TYPE);
+ }
}
- if (flag3->answer) {
- if (G_get_raster_row(infd_tair, inrast_tair, row, data_type_tair)
- < 0)
- G_fatal_error(_("Could not read from <%s>"), tair);
- if (G_get_raster_row(infd_e0, inrast_e0, row, data_type_e0) < 0)
- G_fatal_error(_("Could not read from <%s>"), e0);
- }
+ /*********************************************/
+ /* If ET Potential/Reference is NOT an input */
+ /* compute ETPOTD and ETa Senay*/
+ else{
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d;
+ DCELL d_albedo,d_tempk,d_lat,d_doy,d_tsw,d_solar;
+ DCELL d_rnetd,d_slope,d_aspect,d_tair;
+ DCELL d_e0,d_ndvi,d_etpotd,d_evapfr;
+
+ G_percent(row, nrows, 2);
- /*process the data */
- for (col = 0; col < ncols; col++)
- {
- switch (data_type_albedo) {
- case CELL_TYPE:
- d_albedo = (double)((CELL *) inrast_albedo)[col];
- break;
- case FCELL_TYPE:
- d_albedo = (double)((FCELL *) inrast_albedo)[col];
- break;
- case DCELL_TYPE:
- d_albedo = (double)((DCELL *) inrast_albedo)[col];
- break;
- }
- switch (data_type_tempk) {
- case CELL_TYPE:
- d_tempk = (double)((CELL *) inrast_tempk)[col];
- break;
- case FCELL_TYPE:
- d_tempk = (double)((FCELL *) inrast_tempk)[col];
- break;
- case DCELL_TYPE:
- d_tempk = (double)((DCELL *) inrast_tempk)[col];
- break;
- }
- switch (data_type_lat) {
- case CELL_TYPE:
- d_lat = (double)((CELL *) inrast_lat)[col];
- break;
- case FCELL_TYPE:
- d_lat = (double)((FCELL *) inrast_lat)[col];
- break;
- case DCELL_TYPE:
- d_lat = (double)((DCELL *) inrast_lat)[col];
- break;
- }
- switch (data_type_doy) {
- case CELL_TYPE:
- d_doy = (double)((CELL *) inrast_doy)[col];
- break;
- case FCELL_TYPE:
- d_doy = (double)((FCELL *) inrast_doy)[col];
- break;
- case DCELL_TYPE:
- d_doy = (double)((DCELL *) inrast_doy)[col];
- break;
- }
- switch (data_type_tsw) {
- case CELL_TYPE:
- d_tsw = (double)((CELL *) inrast_tsw)[col];
- break;
- case FCELL_TYPE:
- d_tsw = (double)((FCELL *) inrast_tsw)[col];
- break;
- case DCELL_TYPE:
- d_tsw = (double)((DCELL *) inrast_tsw)[col];
- break;
- }
- if (flag2->answer) {
- switch (data_type_slope) {
- case CELL_TYPE:
- d_slope = (double)((CELL *) inrast_slope)[col];
- break;
- case FCELL_TYPE:
- d_slope = (double)((FCELL *) inrast_slope)[col];
- break;
- case DCELL_TYPE:
+ /* read input maps */
+ Rast_get_row(infd_albedo, inrast_albedo, row, DCELL_TYPE);
+ Rast_get_row(infd_tempk, inrast_tempk, row, DCELL_TYPE);
+ Rast_get_row(infd_lat, inrast_lat, row, DCELL_TYPE);
+ Rast_get_row(infd_doy, inrast_doy, row, DCELL_TYPE);
+ Rast_get_row(infd_tsw, inrast_tsw, row, DCELL_TYPE);
+ if (flag2->answer) {
+ Rast_get_row(infd_slope, inrast_slope, row, DCELL_TYPE);
+ Rast_get_row(infd_aspect, inrast_aspect, row, DCELL_TYPE);
+ }
+ if (flag3->answer) {
+ Rast_get_row(infd_tair, inrast_tair, row, DCELL_TYPE);
+ Rast_get_row(infd_e0, inrast_e0, row, DCELL_TYPE);
+ }
+
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ d_albedo = (double)((DCELL *) inrast_albedo)[col];
+ d_tempk = (double)((DCELL *) inrast_tempk)[col];
+ d_lat = (double)((DCELL *) inrast_lat)[col];
+ d_doy = (double)((DCELL *) inrast_doy)[col];
+ d_tsw = (double)((DCELL *) inrast_tsw)[col];
+ if (flag2->answer) {
d_slope = (double)((DCELL *) inrast_slope)[col];
- break;
- }
- switch (data_type_aspect) {
- case CELL_TYPE:
- d_aspect = (double)((CELL *) inrast_aspect)[col];
- break;
- case FCELL_TYPE:
- d_aspect = (double)((FCELL *) inrast_aspect)[col];
- break;
- case DCELL_TYPE:
d_aspect = (double)((DCELL *) inrast_aspect)[col];
- break;
- }
- }
- if (flag3->answer) {
- switch (data_type_tair) {
- case CELL_TYPE:
- d_tair = (double)((CELL *) inrast_tair)[col];
- break;
- case FCELL_TYPE:
- d_tair = (double)((FCELL *) inrast_tair)[col];
- break;
- case DCELL_TYPE:
+ }
+ if (flag3->answer) {
d_tair = (double)((DCELL *) inrast_tair)[col];
- break;
- }
- switch (data_type_e0) {
- case CELL_TYPE:
- d_e0 = (double)((CELL *) inrast_e0)[col];
- break;
- case FCELL_TYPE:
- d_e0 = (double)((FCELL *) inrast_e0)[col];
- break;
- case DCELL_TYPE:
d_e0 = (double)((DCELL *) inrast_e0)[col];
- break;
- }
- }
- if (G_is_d_null_value(&d_albedo) || G_is_d_null_value(&d_tempk)
- || d_tempk < 10.0 || G_is_d_null_value(&d_lat) ||
- G_is_d_null_value(&d_doy) || G_is_d_null_value(&d_tsw) ||
- G_is_d_null_value(&d_ndvi) || ((flag2->answer) &&
- G_is_d_null_value(&d_slope))
- || ((flag2->answer) && G_is_d_null_value(&d_aspect)) ||
- ((flag3->answer) && G_is_d_null_value(&d_tair)) ||
- ((flag3->answer) && G_is_d_null_value(&d_e0))) {
- G_set_d_null_value(&outrast1[col], 1);
- if (result2)
- G_set_d_null_value(&outrast2[col], 1);
- }
- else {
- if (flag2->answer) {
- d_solar =
- solar_day_3d(d_lat, d_doy, d_tsw, d_slope, d_aspect);
- }
- else {
- d_solar = solar_day(d_lat, d_doy, d_tsw);
- }
- d_evapfr = evapfr_senay(tempk_max, tempk_min, d_tempk);
-
- /*If water then no water stress */
- if (d_albedo <= 0.1 && d_ndvi <= 0.0) {
- d_evapfr = 1.0;
- }
-
- /*some points are colder than selected low */
- if (d_evapfr > 1.0) {
- d_evapfr = 1.0;
- }
- if (result2) {
- outrast2[col] = d_evapfr;
- }
- if (flag3->answer) {
- d_rnetd =
- r_net_day_bandara98(d_albedo, d_solar, e_atm, d_e0,
- d_tair);
- }
- else {
- d_rnetd = r_net_day(d_albedo, d_solar, d_tsw);
- }
- d_etpotd = et_pot_day(d_rnetd, d_tempk, roh_w);
- d = d_etpotd * d_evapfr;
- outrast1[col] = d;
- }
- }
- if (G_put_raster_row(outfd1, outrast1, data_type_output) < 0)
- G_fatal_error(_("Cannot write to output raster file"));
- if (result2) {
- if (G_put_raster_row(outfd2, outrast2, data_type_output) < 0)
- G_fatal_error(_("Cannot write to output raster file"));
- }
- }
+ }
+ if (Rast_is_d_null_value(&d_albedo) || Rast_is_d_null_value(&d_tempk) ||
+ d_tempk < 10.0 || Rast_is_d_null_value(&d_lat) ||
+ Rast_is_d_null_value(&d_doy) || Rast_is_d_null_value(&d_tsw) ||
+ Rast_is_d_null_value(&d_ndvi) ||
+ ((flag2->answer) && Rast_is_d_null_value(&d_slope)) ||
+ ((flag2->answer) && Rast_is_d_null_value(&d_aspect)) ||
+ ((flag3->answer) && Rast_is_d_null_value(&d_tair)) ||
+ ((flag3->answer) && Rast_is_d_null_value(&d_e0))) {
+ Rast_set_d_null_value(&outrast1[col], 1);
+ if (result2)
+ Rast_set_d_null_value(&outrast2[col], 1);
+ }
+ else {
+ if (flag2->answer) {
+ d_solar =
+ solar_day_3d(d_lat, d_doy, d_tsw, d_slope, d_aspect);
+ }
+ else {
+ d_solar = solar_day(d_lat, d_doy, d_tsw);
+ }
+ d_evapfr = evapfr_senay(tempk_max, tempk_min, d_tempk);
+ /*If water then no water stress */
+ if (d_albedo <= 0.1 && d_ndvi <= 0.0) {
+ d_evapfr = 1.0;
+ }
+ /*some points are colder than selected low */
+ if (d_evapfr > 1.0) {
+ d_evapfr = 1.0;
+ }
+ if (result2) {
+ outrast2[col] = d_evapfr;
+ }
+ if (flag3->answer) {
+ d_rnetd =
+ r_net_day_bandara98(d_albedo,
+ d_solar, e_atm, d_e0,d_tair);
+ }
+ else {
+ d_rnetd = r_net_day(d_albedo, d_solar, d_tsw);
+ }
+ d_etpotd = et_pot_day(d_rnetd, d_tempk, roh_w);
+ d = d_etpotd * d_evapfr;
+ outrast1[col] = d;
+ }
+ }
+ Rast_put_row(outfd1, outrast1, DCELL_TYPE);
+ if (result2)
+ Rast_put_row(outfd2, outrast2, DCELL_TYPE);
+ }
+ }
+ }
G_free(inrast_albedo);
G_free(inrast_tempk);
G_free(inrast_dem);
G_free(inrast_lat);
G_free(inrast_doy);
G_free(inrast_tsw);
- G_close_cell(infd_albedo);
- G_close_cell(infd_tempk);
- G_close_cell(infd_dem);
- G_close_cell(infd_lat);
- G_close_cell(infd_doy);
- G_close_cell(infd_tsw);
+ Rast_close(infd_albedo);
+ Rast_close(infd_tempk);
+ Rast_close(infd_dem);
+ Rast_close(infd_lat);
+ Rast_close(infd_doy);
+ Rast_close(infd_tsw);
if (flag3->answer) {
G_free(inrast_tair);
- G_close_cell(infd_tair);
+ Rast_close(infd_tair);
G_free(inrast_e0);
- G_close_cell(infd_e0);
+ Rast_close(infd_e0);
}
G_free(outrast1);
- G_close_cell(outfd1);
+ Rast_close(outfd1);
if (result2) {
G_free(outrast2);
- G_close_cell(outfd2);
+ Rast_close(outfd2);
- /* Color table for evaporative fraction */
- G_init_colors(&colors);
- G_add_color_rule(0, 0, 0, 0, 1, 255, 255, 255, &colors);
- G_short_history(result2, "raster", &history);
- G_command_history(&history);
- G_write_history(result2, &history);
+ /* Color table for evaporative fraction */
+ Rast_init_colors(&colors);
+ val1=0;
+ val2=1;
+ Rast_add_c_color_rule(&val1, 0, 0, 0, &val2, 255, 255, 255, &colors);
+ Rast_short_history(result2, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result2, &history);
}
/* Color table for evapotranspiration */
- G_init_colors(&colors);
- G_add_color_rule(0, 0, 0, 0, 10, 255, 255, 255, &colors);
- G_short_history(result1, "raster", &history);
- G_command_history(&history);
- G_write_history(result1, &history);
+ Rast_init_colors(&colors);
+ val1=0;
+ val2=20;
+ Rast_add_c_color_rule(&val1, 0, 0, 0, &val2, 255, 255, 255, &colors);
+ Rast_short_history(result1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result1, &history);
+
exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list