[GRASS-SVN] r51066 - in grass-addons/grass7/imagery: .
i.eb.h_SEBAL95
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Mar 15 05:48:13 EDT 2012
Author: ychemin
Date: 2012-03-15 02:48:13 -0700 (Thu, 15 Mar 2012)
New Revision: 51066
Added:
grass-addons/grass7/imagery/i.eb.h_SEBAL95/
grass-addons/grass7/imagery/i.eb.h_SEBAL95/Makefile
grass-addons/grass7/imagery/i.eb.h_SEBAL95/U_0.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair_desert.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/functions.h
grass-addons/grass7/imagery/i.eb.h_SEBAL95/h1.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/h_0.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/i.eb.h_SEBAL95.html
grass-addons/grass7/imagery/i.eb.h_SEBAL95/main.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_h.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_m.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah1.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah_0.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air_0.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/sensi_h.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/u_star.c
grass-addons/grass7/imagery/i.eb.h_SEBAL95/zom_0.c
grass-addons/grass7/imagery/i.evapo.potrad/
Log:
more maintenance on gipe, upgrade some modules to G7 addons
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/Makefile 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.eb.h_SEBAL95
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/U_0.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/U_0.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/U_0.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double U_0(double zom_0, double u2m)
+{
+ double u_0;
+
+ u_0 =
+ u2m * 0.41 * log10(200 / (0.15 / 7)) / (log10(2 / (0.15 / 7)) *
+ log10(200 / zom_0));
+
+ return (u_0);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,19 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+/* Pixel-based input required are: tempk water & desert
+ * additionally, dtair in Desert point should be given
+ */
+
+
+double dt_air(double t0_dem, double tempk_water, double tempk_desert,
+ double dtair_desert)
+{
+ double a, b, result;
+
+ a = (dtair_desert - 0.0) / (tempk_desert - tempk_water);
+ b = dtair_desert - a * tempk_desert;
+ result = t0_dem * a + b;
+ return result;
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair_desert.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair_desert.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/dtair_desert.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include"functions.h"
+ double dt_air_desert(double h_desert, double roh_air_desert,
+ double rah_desert)
+{
+ double result;
+
+ result = (h_desert * rah_desert) / (roh_air_desert * 1004);
+ return result;
+}
+
+
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/functions.h
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/functions.h (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/functions.h 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,82 @@
+/* These are the headers for the RS functions used in SEBAL */
+/* 2004 */
+
+/* Initial functions */
+double DN2Rad_landsat7(double Lmin, double LMax, double QCalMax,
+ double QCalmin, int DN);
+double Rad2Ref_landsat7(double radiance, double doy, double sun_elevation,
+ double k_exo);
+double tempk_landsat7(double chan6);
+
+double bb_alb_aster(double greenchan, double redchan, double nirchan,
+ double swirchan1, double swirchan2, double swirchan3,
+ double swirchan4, double swirchan5, double swirchan6);
+double bb_alb_landsat(double bluechan, double greenchan, double redchan,
+ double nirchan, double chan5, double chan7);
+double bb_alb_noaa(double redchan, double nirchan);
+
+double bb_alb_modis(double redchan, double nirchan, double chan3,
+ double chan4, double chan5, double chan6, double chan7);
+double nd_vi(double redchan, double nirchan);
+
+double emissivity_generic(double ndvi);
+
+double emissivity_modis(double e31, double e32);
+
+double t0_dem(double dem, double tempk);
+
+double t_air_approx(double temp_k);
+
+/* Instantaneous rnet and g0 */
+double r_net(double bbalb, double ndvi, double tempk, double tair, double e0,
+ double tsw, double doy, double utc, double sunzangle);
+double g_0(double bbalb, double ndvi, double tempk, double rnet, double utc);
+
+/* Diurnal r_net and etpot */
+double solar_day(double lat, double doy, double tsw);
+
+double r_net_day(double bbalb, double solar, double tsw);
+
+double et_pot_day(double bbalb, double solar, double tempk, double tsw);
+
+/* Sensible heat flux functions */
+double sensi_h(int iteration, double tempk_water, double tempk_desert,
+ double t0_dem, double tempk, double ndvi, double ndvi_max,
+ double dem, double rnet_desert, double g0_desert,
+ double t0_dem_desert, double u2m, double dem_desert);
+double roh_air_0(double tempk);
+
+double zom_0(double ndvi, double ndvi_max);
+
+double U_0(double zom_0, double u2m);
+
+double rah_0(double zom_0, double u_0);
+
+double h_0(double roh_air, double rah, double dtair);
+
+double dt_air_approx(double temp_k);
+
+double dt_air_0(double t0_dem, double tempk_water, double tempk_desert);
+
+double dt_air_desert(double h_desert, double roh_air_desert,
+ double rah_desert);
+double dt_air(double t0_dem, double tempk_water, double tempk_desert,
+ double dtair_desert);
+double rohair(double dem, double tempk, double dtair);
+
+double h1(double roh_air, double rah, double dtair);
+
+double u_star(double t0_dem, double h, double ustar, double roh_air,
+ double zom, double u2m);
+double psi_h(double t0_dem, double h, double U_0, double roh_air);
+
+double rah1(double psih, double u_star);
+
+/* Final outputs */
+double evap_fr(double r_net, double g0, double h);
+
+double et_a(double r_net_day, double evap_f, double tempk);
+
+double biomass(double ndvi, double solar_day, double evap_fr,
+ double light_use_ef);
+double soilmoisture(double ndvi, double evap_fr);
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/h1.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/h1.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/h1.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,12 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double h1(double roh_air, double rah, double dtair)
+{
+ double result;
+
+ result = roh_air * 1004 * dtair / rah;
+
+ return result;
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/h_0.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/h_0.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/h_0.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,12 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double h_0(double roh_air, double rah, double dtair)
+{
+ double result;
+
+ result = roh_air * 1004 * (dtair) / rah;
+
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/i.eb.h_SEBAL95.html
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/i.eb.h_SEBAL95.html (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/i.eb.h_SEBAL95.html 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,141 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<h2>NAME</h2> <I>i.eb.h_SEBAL95 </I>- computation of <i>sensible heat flux</i> [W/m2] after Bastiaanssen, 1995 in [1].
+
+<p><I>(GRASS Raster Program)</I>
+
+<h2>SYNOPSIS</h2>
+<b>i.eb.h_SEBAL95</b>
+<br>
+
+<b>i.eb.h_SEBAL95</b> help
+<br>
+
+<b>i.eb.h_SEBAL95</b> <b>[ -qzn ]</b>
+
+<b>DEM</b>=name
+<b>T</b>=name
+<b>RH</b>=name
+<b>WS</b>=name
+<b>NSR</b>=name
+<b>Vh</b>=name
+<b>ETP</b>=name
+
+<h2>DESCRIPTION</h2>
+
+<p><em>i.eb.h_SEBAL95</em> given the vegetation height (hc), humidity (RU),
+wind speed at two meters height (WS), temperature (T), digital terrain model (DEM),
+and net radiation (NSR) raster input maps,
+calculates the sensible heat flux map (h0).
+
+<p>Optionally the user can activate a flag (-z)
+that allows him setting to zero all of the negative evapotranspiration cells;
+in fact these negative values motivated by the condensation of the air water
+vapour content, are sometime undesired because they can produce computational
+problems. The usage of the flag -n detect that the module is run in night hours
+and the appropriate soil heat flux is calculated.
+
+<p>The algorithm implements well known approaches: the hourly Penman-Monteith method as presented in Allen et al. (1998) for land surfaces and the Penman method (Penman, 1948) for water surfaces.<br>
+
+<p>Land and water surfaces are idenfyied by Vh:<br>
+- where Vh less than 0 vegetation is present and evapotranspiration is calculated;<br>
+- where Vh=0 bare ground is present and evapotranspiration is calculated;<br>
+- where Vh more than 0 water surface is present and evaporation is calculated;<br>
+
+<p>For more details on the algorithms see [1].
+
+
+<h2>OPTIONS</h2>
+
+The program will run non-interactively if the user specifies program
+arguments and flag settings on the command line using the following
+form:
+
+<p><b>i.eb.h_SEBAL95</b> <b>[ -qzd ]</b>
+<b>DEM</b>=name
+<b>T</b>=name
+<b>RH</b>=name
+<b>WS</b>=name
+<b>NSR</b>=name
+<b>Vh</b>=name
+<b>ETP</b>=name
+
+
+
+<p>Alternatively, the user can simply type <em>i.eb.h_SEBAL95</em> on the
+command line and the program will ask for parameter values and flag
+settings interactively, using the standard GRASS parser interface.
+
+
+<h3>Parameters:</h3>
+<dl>
+ <dt><b>DEM</b>=<I>name</I>
+ <dd>Input elevation raster [m a.s.l.]. Required.</dd>
+
+ <dt><b>T</b>=<I>name</I>
+ <dd>Input temperature raster [�C]. Required.</dd>
+
+ <dt><b>RH</b> =<I>name</I>
+ <dd>Input relative humidity raster [%]. Required.</dd>
+
+ <dt><b>WS</b> =<I>name</I>
+ <dd>Input wind speed at two meters raster [m/s]. Required.</dd>
+
+ <dt><b>NSR</b> =<I>name</I>
+ <dd>Input net solar radiation raster [MJ/(m2*h)]. Required.</dd>
+
+ <dt><b>Vh</b> =<I>name</I>
+ <dd>Input vegetation heigth raster [m]. Required.</dd>
+
+ <dt><b>ETP</b> =<I>name</I>
+ <dd>Output evapotranspiration raster [mm/h]. Required.</dd>
+
+</dl>
+
+
+<h2>NOTES</h2>
+
+<p>Net solar radiation map in MJ/(m2*h) can be computed from the combination of the r.sun ,
+run in mode 1, and the r.mapcalc commands.
+
+<p>The sum of the three radiation components outputted by r.sun (beam, diffuse, and reflected)
+multiplied by the Wh to Mj conversion factor (0.0036) and optionally by a
+clear sky factor [0-1] allows the generation of a map to be used as
+an NSR input for the <em>i.evapo.pm</em> command.
+<br>
+example:
+<br>r.sun -s elevin=dem aspin=aspect slopein=slope lin=2 albedo=alb_Mar incidout=out beam_rad=beam diff_rad=diffuse refl_rad=reflected day=73 time=13:00 dist=100;
+<br>r.mapcalc 'NSR=0.0036*(beam+diffuse+reflected)';
+
+
+<h2>SEE ALSO</h2>
+<ul>
+ <li><a href=i.eb.h_iter.html>i.eb.h_iter</a>,
+ <a href=i.eb.h0.html>i.eb.h0</a>
+ <a href=i.evapo.pm.html>i.evapo.pm</a>
+</ul>
+
+
+
+<h2>AUTHORS</h2>
+ <p>
+ <i>
+ <br>Yann Chemin
+ <br>International Rice Research Institute, Los Banos, The Philippines.
+ <br>International Water management Institute, Colombo, Sri Lanka.
+ </i>
+ <p>Contact: <a href="mailto:y.chemin at cgiar.org"> Yann chemin</a>
+
+
+<h2>REFERENCES</h2>
+
+ <p>[1] Bastiaanssen, W.G.M., 1995.
+ Estimation of Land surface paramters by remote sensing under clear-sky conditions. PhD thesis, Wageningen University, Wageningen, The Netherlands.
+
+ <p>[2] Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998.
+ Crop Evapotranspiration: Guidelines for computing crop water requirements.
+ Irrigation and Drainage Paper 56, Food and Agriculture Organization of the
+ United Nations, Rome, pp. 300
+
+ <p>[3] Penman, H. L. 1948. Natural evaporation from open water,
+ bare soil and grass. Proc. Roy. Soc. London, A193, pp. 120-146.
+
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/main.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/main.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,920 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.eb.h_SEBAL95
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculates sensible heat flux by SEBAL iteration
+ * Delta T will be reassessed in the iterations !
+ * This has been seen in Bastiaanssen (1995).
+ *
+ * COPYRIGHT: (C) 2002-2009 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.
+ *
+ * CHANGELOG:
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+double sensi_h(int iteration, double tempk_water, double tempk_desert,
+ double t0_dem, double tempk, double ndvi, double ndvi_max,
+ double dem, double rnet_desert, double g0_desert,
+ double t0_dem_desert, double u2m, double dem_desert);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd;
+ char *mapset; /* mapset name */
+ /* buffer for in out raster */
+ void *inrast_T, *inrast_ndvi, *inrast_dem, *inrast_u2m;
+ void *inrast_Rn, *inrast_g0, *inrast_albedo;
+ DCELL *outrast;
+ int nrows, ncols;
+ int row, col;
+ int row_wet, col_wet;
+ int row_dry, col_dry;
+ double m_row_wet, m_col_wet;
+ double m_row_dry, m_col_dry;
+ int infd_T, infd_ndvi, infd_dem, infd_u2m;
+ int infd_Rn, infd_g0, infd_albedo;
+ int outfd;
+ char *mapset_T, *mapset_ndvi, *mapset_dem, *mapset_u2m;
+ char *mapset_Rn, *mapset_g0, *mapset_albedo;
+ char *T, *ndvi, *dem, *u2m, *Rn, *g0, *albedo;
+ char *h0;
+ struct History history;
+ struct GModule *module;
+ struct Option *input_T, *input_ndvi, *input_dem, *input_u2m;
+ struct Option *input_Rn, *input_g0, *input_albedo;
+ struct Option *output;
+ struct Option *input_row_wet, *input_col_wet;
+ struct Option *input_row_dry, *input_col_dry;
+ struct Option *input_iter;
+ struct Flag *flag1, *flag2, *flag3, *day, *zero;
+ /*******************************/
+ RASTER_MAP_TYPE data_type_T;
+ RASTER_MAP_TYPE data_type_ndvi;
+ RASTER_MAP_TYPE data_type_dem;
+ RASTER_MAP_TYPE data_type_u2m;
+ RASTER_MAP_TYPE data_type_Rn;
+ RASTER_MAP_TYPE data_type_g0;
+ RASTER_MAP_TYPE data_type_albedo;
+ RASTER_MAP_TYPE data_type_output = DCELL_TYPE;
+ /*******************************/
+ int iteration = 10; /*SEBAL95 loop number */
+ /********************************/
+ /* Stats for dry/wet pixels */
+ double t0dem_min = 400.0, t0dem_max = 200.0;
+ double tempk_min = 400.0, tempk_max = 200.0;
+ /********************************/
+ double xp, yp;
+ double xmin, ymin;
+ double xmax, ymax;
+ double stepx, stepy;
+ double latitude, longitude;
+ /********************************/
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->description = _("Sensible Heat Flux iteration SEBAL 95");
+
+ /* Define different options */
+ input_T = G_define_standard_option(G_OPT_R_INPUT);
+ input_T->key = "tempk";
+ input_T->description =
+ _("Name of Surface Skin Temperature input map [K]");
+
+ input_dem = G_define_standard_option(G_OPT_R_INPUT);
+ input_dem->key = "dem";
+ input_dem->description = _("Name of dem input map [m a.s.l.]");
+
+ input_u2m = G_define_standard_option(G_OPT_R_INPUT);
+ input_u2m->key = "u2m";
+ input_u2m->description =
+ _("Name of Wind speed at 2m height input map [m/s]");
+
+ input_ndvi = G_define_standard_option(G_OPT_R_INPUT);
+ input_ndvi->key = "ndvi";
+ input_ndvi->description = _("Name of NDVI input map [-]");
+
+ input_albedo = G_define_standard_option(G_OPT_R_INPUT);
+ input_albedo->key = "albedo";
+ input_albedo->description = _("Name of Albedo input map [-]");
+
+ input_Rn = G_define_standard_option(G_OPT_R_INPUT);
+ input_Rn->key = "rnet";
+ input_Rn->description =
+ _("Name of instantaneous Net Solar Radiation input map [W/m2]");
+
+ input_g0 = G_define_standard_option(G_OPT_R_INPUT);
+ input_g0->key = "g0";
+ input_g0->description =
+ _("Name of instantaneous Soil Heat Flux input map [W/m2]");
+
+ input_iter = G_define_option();
+ input_iter->key = "iteration";
+ input_iter->type = TYPE_INTEGER;
+ input_iter->required = NO;
+ input_iter->gisprompt = "old,value";
+ input_iter->description =
+ _("Value of the number of SEBAL95 loops (default is 10)");
+ input_iter->guisection = _("Optional");
+
+ input_row_wet = G_define_option();
+ input_row_wet->key = "row_wet";
+ input_row_wet->type = TYPE_INTEGER;
+ input_row_wet->required = NO;
+ input_row_wet->gisprompt = "old,value";
+ input_row_wet->description = _("Row value of the wet pixel");
+ input_row_wet->guisection = _("Optional");
+
+ input_col_wet = G_define_option();
+ input_col_wet->key = "col_wet";
+ input_col_wet->type = TYPE_INTEGER;
+ input_col_wet->required = NO;
+ input_col_wet->gisprompt = "old,value";
+ input_col_wet->description = _("Column value of the wet pixel");
+ input_col_wet->guisection = _("Optional");
+
+ input_row_dry = G_define_option();
+ input_row_dry->key = "row_dry";
+ input_row_dry->type = TYPE_INTEGER;
+ input_row_dry->required = NO;
+ input_row_dry->gisprompt = "old,value";
+ input_row_dry->description = _("Row value of the dry pixel");
+ input_row_dry->guisection = _("Optional");
+
+ input_col_dry = G_define_option();
+ input_col_dry->key = "col_dry";
+ input_col_dry->type = TYPE_INTEGER;
+ input_col_dry->required = NO;
+ input_col_dry->gisprompt = "old,value";
+ input_col_dry->description = _("Column value of the dry pixel");
+ input_col_dry->guisection = _("Optional");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->description = _("Name of output sensible heat flux layer [W/m2]");
+
+ /* Define the different flags */
+ flag1 = G_define_flag();
+ flag1->key = 't';
+ flag1->description = _("Temperature histogram check (careful!)");
+
+ flag2 = G_define_flag();
+ flag2->key = 'a';
+ flag2->description = _("Automatic wet/dry pixel (careful!)");
+
+ flag3 = G_define_flag();
+ flag3->key = 'c';
+ flag3->description =
+ _("Coordinates of manual dry/wet pixels are in image projection and not row/col");
+
+ zero = G_define_flag();
+ zero->key = 'z';
+ zero->description = _("set negative evapo to zero");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get entered parameters */
+ T = input_T->answer;
+ dem = input_dem->answer;
+ u2m = input_u2m->answer;
+ ndvi = input_ndvi->answer;
+ Rn = input_Rn->answer;
+ g0 = input_g0->answer;
+ albedo = input_albedo->answer;
+
+ h0 = output->answer;
+
+ if (input_iter->answer) {
+ iteration = atoi(input_iter->answer);
+ }
+ if (input_row_wet->answer && input_row_dry &&
+ input_col_wet->answer && input_col_dry) {
+ m_row_wet = atof(input_row_wet->answer);
+ m_col_wet = atof(input_col_wet->answer);
+ m_row_dry = atof(input_row_dry->answer);
+ m_col_dry = atof(input_col_dry->answer);
+ if (flag3->answer) {
+ G_message("Manual wet/dry pixels in image coordinates");
+ G_message("Wet Pixel=> x:%f y:%f", m_col_wet, m_row_wet);
+ G_message("Dry Pixel=> x:%f y:%f", m_col_dry, m_row_dry);
+ }
+ else {
+ G_message("Wet Pixel=> row:%.0f col:%.0f", m_row_wet, m_col_wet);
+ G_message("Dry Pixel=> row:%.0f col:%.0f", m_row_dry, m_col_dry);
+ }
+ }
+ /* determine the input map type (CELL/FCELL/DCELL) */
+ data_type_T = Rast_map_type(T, "");
+ data_type_dem = Rast_map_type(dem, "");
+ data_type_u2m = Rast_map_type(u2m, "");
+ data_type_ndvi = Rast_map_type(ndvi, "");
+ data_type_Rn = Rast_map_type(Rn, "");
+ data_type_g0 = Rast_map_type(g0, "");
+ data_type_albedo = Rast_map_type(albedo, "");
+
+ infd_T = Rast_open_old(T, "");
+ infd_dem = Rast_open_old(dem, "");
+ infd_u2m = Rast_open_old(u2m, "");
+ infd_ndvi = Rast_open_old(ndvi, "");
+ infd_Rn = Rast_open_old(Rn, "");
+ infd_g0 = Rast_open_old(g0, "");
+ infd_albedo = Rast_open_old(albedo, "");
+
+ Rast_get_cellhd(T, "", &cellhd);
+ Rast_get_cellhd(dem, "", &cellhd);
+ Rast_get_cellhd(u2m, "", &cellhd);
+ Rast_get_cellhd(ndvi, "", &cellhd);
+ Rast_get_cellhd(Rn, "", &cellhd);
+ Rast_get_cellhd(g0, "", &cellhd);
+ Rast_get_cellhd(albedo, "", &cellhd);
+
+ /* Allocate input buffer */
+ inrast_T = Rast_allocate_buf(data_type_T);
+ inrast_dem = Rast_allocate_buf(data_type_dem);
+ inrast_u2m = Rast_allocate_buf(data_type_u2m);
+ inrast_ndvi = Rast_allocate_buf(data_type_ndvi);
+ inrast_Rn = Rast_allocate_buf(data_type_Rn);
+ inrast_g0 = Rast_allocate_buf(data_type_g0);
+ inrast_albedo = Rast_allocate_buf(data_type_albedo);
+ /***************************************************/
+ /* Setup pixel location variables */
+ /***************************************************/
+ stepx = cellhd.ew_res;
+ stepy = cellhd.ns_res;
+
+ xmin = cellhd.west;
+ xmax = cellhd.east;
+ ymin = cellhd.south;
+ ymax = cellhd.north;
+
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ /***************************************************/
+ /* Allocate output buffer */
+ /***************************************************/
+ outrast = Rast_allocate_d_buf();
+ outfd = Rast_open_new(h0, DCELL_TYPE);
+ DCELL d_ndvi; /* Input raster */
+ DCELL d_ndvi_max = 0.0; /* Generated here */
+
+ /* THREAD 1 */
+ /* NDVI Max */
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ Rast_get_d_row(infd_ndvi, inrast_ndvi, row);
+ for (col = 0; col < ncols; col++) {
+ 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 (Rast_is_d_null_value(&d_ndvi)) {
+ /* do nothing */
+ }
+ else if ((d_ndvi) > d_ndvi_max && (d_ndvi) < 0.999) {
+ d_ndvi_max = d_ndvi;
+ }
+ }
+ }
+ G_message("ndvi_max=%f\n", d_ndvi_max);
+ /* Pick up wet and dry pixel values */
+ DCELL d_Rn; /* Input raster */
+ DCELL d_g0; /* Input raster */
+ DCELL d_tempk_wet;
+ DCELL d_tempk_dry;
+ DCELL d_Rn_dry;
+ DCELL d_g0_dry;
+ DCELL d_t0dem_dry;
+ DCELL d_dem_dry;
+ DCELL d_dT_dry;
+ DCELL d_ndvi_dry;
+ DCELL d_dT;
+
+ /*START Temperature minimum search */
+ /* THREAD 1 */
+ /*This is correcting for un-Earthly temperatures */
+ /*It finds when histogram is actually starting to pull up... */
+ int peak1, peak2, peak3;
+ int i_peak1, i_peak2, i_peak3;
+ int bottom1a, bottom1b;
+ int bottom2a, bottom2b;
+ int bottom3a, bottom3b;
+ int i_bottom1a, i_bottom1b;
+ int i_bottom2a, i_bottom2b;
+ int i_bottom3a, i_bottom3b;
+ if (flag1->answer) {
+ int i = 0;
+ int histogram[400];
+ for (i = 0; i < 400; i++) {
+ histogram[i] = 0;
+ }
+ DCELL d_T;
+ /****************************/
+ /* Process pixels histogram */
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+ /* read input map */
+ Rast_get_row(infd_T, inrast_T, row, data_type_T);
+ /*process the data */
+ for (col = 0; col < ncols; col++) {
+ switch (data_type_T) {
+ case CELL_TYPE:
+ d_T = (double)((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_T = (double)((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_T = (double)((DCELL *) inrast_T)[col];
+ break;
+ }
+ if (Rast_is_d_null_value(&d_T)) {
+ /*Do nothing */
+ }
+ else {
+ int temp;
+
+ temp = (int)d_T;
+ if (temp > 0) {
+ histogram[temp] = histogram[temp] + 1.0;
+ }
+ }
+ }
+ }
+ int sum = 0;
+
+ double average = 0.0;
+
+ for (i = 0; i < 400; i++) {
+ sum += histogram[i];
+ }
+ average = (double)sum / 400.0;
+ G_message
+ ("Histogram of Temperature map (if it has rogue values to clean)");
+ peak1 = 0;
+ peak2 = 0;
+ peak3 = 0;
+ i_peak1 = 0;
+ i_peak2 = 0;
+ i_peak3 = 0;
+ bottom1a = 100000;
+ bottom1b = 100000;
+ bottom2a = 100000;
+ bottom2b = 100000;
+ bottom3a = 100000;
+ bottom3b = 100000;
+ i_bottom1a = 1000;
+ i_bottom1b = 1000;
+ i_bottom2a = 1000;
+ i_bottom2b = 1000;
+ i_bottom3a = 1000;
+ i_bottom3b = 1000;
+ for (i = 0; i < 400; i++) {
+ /* Search for highest peak of dataset (2) */
+ /* Highest Peak */
+ if (histogram[i] > peak2) {
+ peak2 = histogram[i];
+ i_peak2 = i;
+ }
+ }
+ int stop = 0;
+
+ for (i = i_peak2; i > 5; i--) {
+ if (((histogram[i] + histogram[i - 1] + histogram[i - 2] +
+ histogram[i - 3] + histogram[i - 4]) / 5) < histogram[i] &&
+ stop == 0) {
+ bottom2a = histogram[i];
+ i_bottom2a = i;
+ }
+ else if (((histogram[i] + histogram[i - 1] + histogram[i - 2] +
+ histogram[i - 3] + histogram[i - 4]) / 5) >
+ histogram[i] && stop == 0) {
+ /*Search for peaks of datasets (1) */
+ peak1 = histogram[i];
+ i_peak1 = i;
+ stop = 1;
+ }
+ }
+ stop = 0;
+ for (i = i_peak2; i < 395; i++) {
+ if (((histogram[i] + histogram[i + 1] + histogram[i + 2] +
+ histogram[i + 3] + histogram[i + 4]) / 5) < histogram[i] &&
+ stop == 0) {
+ bottom2b = histogram[i];
+ i_bottom2b = i;
+ }
+ else if (((histogram[i] + histogram[i + 1] + histogram[i + 2] +
+ histogram[i + 3] + histogram[i + 4]) / 5) >
+ histogram[i] && stop == 0) {
+ /*Search for peaks of datasets (3) */
+ peak3 = histogram[i];
+ i_peak3 = i;
+ stop = 1;
+ }
+ }
+ /* First histogram lower bound */
+ for (i = 250; i < i_peak1; i++) {
+ if (histogram[i] < bottom1a) {
+ bottom1a = histogram[i];
+ i_bottom1a = i;
+ }
+ }
+ /* First histogram higher bound */
+ for (i = i_peak2; i > i_peak1; i--) {
+ if (histogram[i] <= bottom1b) {
+ bottom1b = histogram[i];
+ i_bottom1b = i;
+ }
+ }
+ /* Third histogram lower bound */
+ for (i = i_peak2; i < i_peak3; i++) {
+ if (histogram[i] < bottom3a) {
+ bottom3a = histogram[i];
+ i_bottom3a = i;
+ }
+ }
+ /* Third histogram higher bound */
+ for (i = 399; i > i_peak3; i--) {
+ if (histogram[i] < bottom3b) {
+ bottom3b = histogram[i];
+ i_bottom3b = i;
+ }
+ }
+ G_message("bottom1a: [%i]=>%i\n", i_bottom1a, bottom1a);
+ G_message("peak1: [%i]=>%i\n", i_peak1, peak1);
+ G_message("bottom1b: [%i]=>%i\n", i_bottom1b, bottom1b);
+ G_message("bottom2a: [%i]=>%i\n", i_bottom2a, bottom2a);
+ G_message("peak2: [%i]=>%i\n", i_peak2, peak2);
+ G_message("bottom2b: [%i]=>%i\n", i_bottom2b, bottom2b);
+ G_message("bottom3a: [%i]=>%i\n", i_bottom3a, bottom3a);
+ G_message("peak3: [%i]=>%i\n", i_peak3, peak3);
+ G_message("bottom3b: [%i]=>%i\n", i_bottom3b, bottom3b);
+ } /*END OF FLAG1 */
+
+ /* End of processing histogram */
+ /*******************/
+
+ /*Process wet pixel values */
+ /* FLAG2 */
+ if (flag2->answer) {
+ /* THREAD 3 */
+ /* Process tempk min / max pixels */
+ /* Internal use only */
+ DCELL d_Rn_wet;
+ DCELL d_g0_wet;
+ /*********************/
+ for (row = 0; row < nrows; row++) {
+ DCELL d_albedo;
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_t0dem;
+ DCELL d_h0 = 100.0; /*for flag 1 */
+ DCELL d_h0_max = 100.0; /*for flag 1 */
+ G_percent(row, nrows, 2);
+ Rast_get_row(infd_albedo, inrast_albedo, row, data_type_albedo);
+ Rast_get_row(infd_T, inrast_T, row, data_type_T);
+ Rast_get_row(infd_dem, inrast_dem, row, data_type_dem);
+ Rast_get_row(infd_Rn, inrast_Rn, row, data_type_Rn);
+ Rast_get_row(infd_g0, inrast_g0, row, data_type_g0);
+ Rast_get_d_row(infd_ndvi, inrast_ndvi, row);
+ /*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_T) {
+ case CELL_TYPE:
+ d_tempk = (double)((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double)((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double)((DCELL *) inrast_T)[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_Rn) {
+ case CELL_TYPE:
+ d_Rn = (double)((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double)((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double)((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch (data_type_g0) {
+ case CELL_TYPE:
+ d_g0 = (double)((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double)((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double)((DCELL *) inrast_g0)[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 (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_Rn) ||
+ Rast_is_d_null_value(&d_g0) || Rast_is_d_null_value(&d_ndvi)) {
+ /* do nothing */
+ }
+ else {
+ d_t0dem = d_tempk + 0.00627 * d_dem;
+ if (d_t0dem <= 250.0 || d_tempk <= 250.0) {
+ /* do nothing */
+ }
+ else {
+ d_h0 = d_Rn - d_g0;
+ if (d_t0dem < t0dem_min && d_albedo < 0.15) {
+ t0dem_min = d_t0dem;
+ tempk_min = d_tempk;
+ d_tempk_wet = d_tempk;
+ d_Rn_wet = d_Rn;
+ d_g0_wet = d_g0;
+ col_wet = col;
+ row_wet = row;
+ }
+ if (flag1->answer &&
+ d_tempk >= (double)i_peak1 - 5.0 &&
+ d_tempk < (double)i_peak1 + 1.0) {
+ tempk_min = d_tempk;
+ d_tempk_wet = d_tempk;
+ d_Rn_wet = d_Rn;
+ d_g0_wet = d_g0;
+ col_wet = col;
+ row_wet = row;
+ }
+ if (d_t0dem > t0dem_max) {
+ t0dem_max = d_t0dem;
+ d_t0dem_dry = d_t0dem;
+ tempk_max = d_tempk;
+ d_tempk_dry = d_tempk;
+ d_Rn_dry = d_Rn;
+ d_g0_dry = d_g0;
+ d_dem_dry = d_dem;
+ col_dry = col;
+ row_dry = row;
+ d_ndvi_dry = d_ndvi;
+ }
+ if (flag1->answer &&
+ d_tempk >= (double)i_peak3 - 0.0 &&
+ d_tempk < (double)i_peak3 + 7.0 &&
+ d_h0 > 100.0 && d_h0 > d_h0_max &&
+ d_g0 > 10.0 && d_Rn > 100.0 && d_albedo > 0.3) {
+ tempk_max = d_tempk;
+ d_tempk_dry = d_tempk;
+ d_Rn_dry = d_Rn;
+ d_g0_dry = d_g0;
+ d_h0_max = d_h0;
+ d_dem_dry = d_dem;
+ col_dry = col;
+ row_dry = row;
+ d_ndvi_dry = d_ndvi;
+ }
+ }
+ }
+ }
+ }
+ G_message("tempk_min=%f\ntempk_max=%f\n", tempk_min, tempk_max);
+ G_message("row_wet=%d\tcol_wet=%d\n", row_wet, col_wet);
+ G_message("row_dry=%d\tcol_dry=%d\n", row_dry, col_dry);
+ G_message("tempk_wet=%f\n", d_tempk_wet);
+ G_message("g0_wet=%f\n", d_g0_wet);
+ G_message("Rn_wet=%f\n", d_Rn_wet);
+ G_message("LE_wet=%f\n", d_Rn_wet - d_g0_wet);
+ G_message("tempk_dry=%f\n", d_tempk_dry);
+ G_message("dem_dry=%f\n", d_dem_dry);
+ G_message("t0dem_dry=%f\n", d_t0dem_dry);
+ G_message("rnet_dry=%f\n", d_Rn_dry);
+ G_message("g0_dry=%f\n", d_g0_dry);
+ G_message("h0_dry=%f\n", d_Rn_dry - d_g0_dry);
+ } /* END OF FLAG2 */
+
+ /* MANUAL WET/DRY PIXELS */
+ if (input_row_wet->answer && input_row_dry->answer &&
+ input_col_wet->answer && input_col_dry->answer) {
+ /*DRY PIXEL */
+ if (flag3->answer) {
+ /*Calculate coordinates of row/col from projected ones */
+ row = (int)((ymax - m_row_dry) / (double)stepy);
+ col = (int)((m_col_dry - xmin) / (double)stepx);
+ G_message("Dry Pixel | row:%i col:%i", row, col);
+ m_row_dry = row;
+ m_col_dry = col;
+ }
+ row = (int)m_row_dry;
+ col = (int)m_col_dry;
+ G_message("Dry Pixel | row:%i col:%i", row, col);
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_t0dem;
+ Rast_get_row(infd_T, inrast_T, row, data_type_T);
+ Rast_get_row(infd_dem, inrast_dem, row, data_type_dem);
+ Rast_get_row(infd_Rn, inrast_Rn, row, data_type_Rn);
+ Rast_get_row(infd_g0, inrast_g0, row, data_type_g0);
+ switch (data_type_T) {
+ case CELL_TYPE:
+ d_tempk = (double)((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double)((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double)((DCELL *) inrast_T)[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_Rn) {
+ case CELL_TYPE:
+ d_Rn = (double)((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double)((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double)((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch (data_type_g0) {
+ case CELL_TYPE:
+ d_g0 = (double)((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double)((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double)((DCELL *) inrast_g0)[col];
+ break;
+ }
+ d_t0dem = d_tempk + 0.00627 * d_dem;
+ d_t0dem_dry = d_t0dem;
+ d_tempk_dry = d_tempk;
+ d_Rn_dry = d_Rn;
+ d_g0_dry = d_g0;
+ d_dem_dry = d_dem;
+ /*WET PIXEL */
+ if (flag3->answer) {
+ /*Calculate coordinates of row/col from projected ones */
+ row = (int)((ymax - m_row_wet) / (double)stepy);
+ col = (int)((m_col_wet - xmin) / (double)stepx);
+ G_message("Wet Pixel | row:%i col:%i", row, col);
+ m_row_wet = row;
+ m_col_wet = col;
+ }
+ row = m_row_wet;
+ col = m_col_wet;
+ G_message("Wet Pixel | row:%i col:%i", row, col);
+ Rast_get_row(infd_T, inrast_T, row, data_type_T);
+ switch (data_type_T) {
+ case CELL_TYPE:
+ d_tempk = (double)((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double)((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double)((DCELL *) inrast_T)[col];
+ break;
+ }
+ d_tempk_wet = d_tempk;
+ G_message("Manual Pixels\n");
+ G_message("***************************\n");
+ G_message("row_wet=%d\tcol_wet=%d\n", (int)m_row_wet, (int)m_col_wet);
+ G_message("row_dry=%d\tcol_dry=%d\n", (int)m_row_dry, (int)m_col_dry);
+ G_message("tempk_wet=%f\n", d_tempk_wet);
+ G_message("tempk_dry=%f\n", d_tempk_dry);
+ G_message("dem_dry=%f\n", d_dem_dry);
+ G_message("t0dem_dry=%f\n", d_t0dem_dry);
+ G_message("rnet_dry=%f\n", d_Rn_dry);
+ G_message("g0_dry=%f\n", d_g0_dry);
+ G_message("h0_dry=%f\n", d_Rn_dry - d_g0_dry);
+ }
+ /* END OF MANUAL WET/DRY PIXELS */
+
+ for (row = 0; row < nrows; row++) {
+ DCELL d_albedo;
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_u2m;
+ DCELL d_t0dem;
+ DCELL d; /* Output pixel */
+ G_percent(row, nrows, 2);
+ /* read a line input maps into buffers */
+ Rast_get_row(infd_albedo, inrast_albedo, row, data_type_albedo);
+ Rast_get_row(infd_T, inrast_T, row, data_type_T);
+ Rast_get_row(infd_dem, inrast_dem, row, data_type_dem);
+ Rast_get_row(infd_u2m, inrast_u2m, row, data_type_u2m);
+ Rast_get_row(infd_ndvi, inrast_ndvi, row, data_type_ndvi);
+ Rast_get_row(infd_Rn, inrast_Rn, row, data_type_Rn);
+ Rast_get_row(infd_g0, inrast_g0, row, data_type_g0);
+ /* read every cell in the line buffers */
+ 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_T) {
+ case CELL_TYPE:
+ d_tempk = (double)((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double)((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double)((DCELL *) inrast_T)[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_u2m) {
+ case CELL_TYPE:
+ d_u2m = (double)((CELL *) inrast_u2m)[col];
+ break;
+ case FCELL_TYPE:
+ d_u2m = (double)((FCELL *) inrast_u2m)[col];
+ break;
+ case DCELL_TYPE:
+ d_u2m = (double)((DCELL *) inrast_u2m)[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;
+ }
+ switch (data_type_Rn) {
+ case CELL_TYPE:
+ d_Rn = (double)((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double)((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double)((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch (data_type_g0) {
+ case CELL_TYPE:
+ d_g0 = (double)((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double)((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double)((DCELL *) inrast_g0)[col];
+ break;
+ }
+ if (Rast_is_d_null_value(&d_tempk) ||
+ Rast_is_d_null_value(&d_dem) ||
+ Rast_is_d_null_value(&d_u2m) ||
+ Rast_is_d_null_value(&d_ndvi) ||
+ Rast_is_d_null_value(&d_Rn) ||
+ Rast_is_d_null_value(&d_g0) ||
+ d_g0 < 0.0 || d_Rn < 0.0 ||
+ d_dem <= -100.0 || d_dem > 9000.0 || d_tempk < 200.0) {
+ Rast_set_d_null_value(&outrast[col], 1);
+ }
+ else {
+ /* Albedo < 0 */
+ if (d_albedo < 0.001) {
+ d_albedo = 0.001;
+ }
+ /* Calculate T0dem */
+ d_t0dem = (double)d_tempk + 0.00627 * (double)d_dem;
+ /* G_message("**InLoop d_t0dem=%5.3f",d_t0dem);
+ G_message(" d_dem=%5.3f",d_dem);
+ G_message(" d_tempk=%5.3f",d_tempk);
+ G_message(" d_albedo=%5.3f",d_albedo);
+ G_message(" d_Rn=%5.3f",d_Rn);
+ G_message(" d_g0=%5.3f",d_g0);
+ G_message(" d_ndvi=%5.3f",d_ndvi);
+ G_message(" d_u_hu=%5.3f",d_u_hu);
+ *//* Calculate sensible heat flux */
+ d = sensi_h(iteration,
+ d_tempk_wet,
+ d_tempk_dry,
+ d_t0dem,
+ d_tempk,
+ d_ndvi,
+ d_ndvi_max,
+ d_dem,
+ d_Rn_dry,
+ d_g0_dry, d_t0dem_dry, d_u2m, d_dem_dry);
+ /* G_message(" d_h0=%5.3f",d); */
+ if (zero->answer && d < 0.0) {
+ d = 0.0;
+ }
+ outrast[col] = d;
+ }
+ }
+ Rast_put_d_row(outfd, outrast);
+ }
+ G_free(inrast_T);
+ Rast_close(infd_T);
+ G_free(inrast_dem);
+ Rast_close(infd_dem);
+ G_free(inrast_u2m);
+ Rast_close(infd_u2m);
+ G_free(inrast_ndvi);
+ Rast_close(infd_ndvi);
+ G_free(inrast_Rn);
+ Rast_close(infd_Rn);
+ G_free(inrast_g0);
+ Rast_close(infd_g0);
+ G_free(inrast_albedo);
+ Rast_close(infd_albedo);
+
+ G_free(outrast);
+ Rast_close(outfd);
+ /* add command line incantation to history file */
+ Rast_short_history(h0, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(h0, &history);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_h.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_h.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_h.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,33 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double psi_h(double t0_dem, double h, double U_0, double roh_air)
+{
+ double result;
+
+ double n5_temp, n11_mem, n12_mem;
+
+ if (h != 0.0) {
+ n5_temp =
+ (-1004 * roh_air * pow(U_0, 3) * t0_dem) / (0.41 * 9.81 * h);
+ }
+ else {
+ n5_temp = -1000.0;
+ }
+
+ if (n5_temp < 0.0) {
+ n12_mem = pow((1 - 16 * (2 / n5_temp)), 0.25);
+ n11_mem = (2 * log10((1 + pow(n12_mem, 2)) / 2));
+ }
+ else {
+ n12_mem = 1.0;
+ n11_mem = -5 * 2 / n5_temp;
+ }
+
+ result = n11_mem;
+
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_m.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_m.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/psi_m.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,31 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double psi_m(double t0_dem, double h, double ustar, double roh_air, double hu)
+{
+ double result;
+
+ double n5, n10, n12;
+
+ if (h != 0.0) {
+ n5 = (-1004 * roh_air * pow(ustar, 3) * t0_dem) / (0.41 * 9.81 * h);
+ }
+ else {
+ n5 = -1000.0;
+ }
+ if (n5 < 0.0) {
+ n12 = pow((1 - 16 * (hu / n5)), 0.25);
+ n10 =
+ (2 * log((1 + n12) / 2)) + log((1 + pow(n12, 2)) / 2) -
+ 2 * atan(n12) + 0.5 * PI;
+ }
+ else {
+ n12 = 1.0;
+ n10 = -5 * hu / n5;
+ }
+ result = n10;
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah1.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah1.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah1.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double rah1(double psih, double ustar)
+{
+ double result;
+
+ result = (log10(2 / 0.01) - psih) / (ustar * 0.41);
+
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah_0.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah_0.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/rah_0.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,12 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double rah_0(double zom_0, double u_0)
+{
+ double result;
+
+ result = log10(2 / 0.01) / (u_0 * 0.41);
+
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double rohair(double dem, double tempk, double dtair)
+{
+ double a, b, result;
+
+ a = tempk - dtair;
+ b = ((a - 0.00627 * dem) / a);
+ result = 349.467 * pow(b, 5.26) / a;
+
+ return result;
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air_0.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air_0.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/roh_air_0.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,13 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double roh_air_0(double tempk)
+{
+ double A, B, result;
+
+ A = (18 * (6.11 * exp(17.27 * 34.8 / (34.8 + 237.3))) / 100.0);
+ B = (18 * (6.11 * exp(17.27 * 34.8 / (34.8 + 237.3))) / 100.0);
+ result = (1000.0 - A) / (tempk * 2.87) + B / (tempk * 4.61);
+ return result;
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/sensi_h.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/sensi_h.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/sensi_h.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,131 @@
+/* This is the main loop used in SEBAL */
+/* April 2004 */
+
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include "functions.h"
+
+/* Arrays Declarations */
+#define ITER_MAX 10
+
+double sensi_h(int iteration, double tempk_water, double tempk_desert,
+ double t0_dem, double tempk, double ndvi, double ndvi_max,
+ double dem, double rnet_desert, double g0_desert,
+ double t0_dem_desert, double u2m, double dem_desert)
+{
+ /* Arrays Declarations */
+ double dtair[ITER_MAX], roh_air[ITER_MAX], rah[ITER_MAX];
+
+ double h[ITER_MAX];
+
+ double ustar[ITER_MAX], zom[ITER_MAX];
+
+ /* Declarations */
+ int i, j, ic, debug = 0;
+
+ double u_0, zom0;
+
+ double h_desert, rah_desert, roh_air_desert;
+
+ double dtair_desert;
+
+ double psih_desert, ustar_desert, ustar_desertold, zom_desert;
+
+ double psih;
+
+ double result;
+
+ /* Fat-free junk food */
+ if (iteration > ITER_MAX) {
+ iteration = ITER_MAX;
+ }
+
+ if (debug == 1) {
+ printf("*****************************\n");
+ printf("t0_dem = %5.3f\n", t0_dem);
+ printf("ndvi = %5.3f ndvimax = %5.3f\n", ndvi, ndvi_max);
+ printf("*****************************\n");
+ }
+ /* dtair[0] = dt_air_0(t0_dem, tempk_water, tempk_desert); */
+ dtair[0] = 5.0;
+ /* printf("*****************************dtair = %5.3f\n",dtair[0]); */
+ roh_air[0] = roh_air_0(tempk);
+ /* printf("*****************************rohair=%5.3f\n",roh_air[0]); */
+ roh_air_desert = roh_air_0(tempk_desert);
+ /* printf("**rohairdesert = %5.3f\n",roh_air_desert); */
+ zom0 = zom_0(ndvi, ndvi_max);
+ /* printf("*****************************zom = %5.3f\n",zom0); */
+ u_0 = U_0(zom0, u2m);
+ /* printf("*****************************u0\n"); */
+ rah[0] = rah_0(zom0, u_0);
+ /* printf("*****************************rah = %5.3f\n",rah[0]); */
+ h[0] = h_0(roh_air[0], rah[0], dtair[0]);
+ /* printf("*****************************h\n"); */
+ if (debug == 1) {
+ printf("dtair[0] = %5.3f K\n", dtair[0]);
+ printf("roh_air[0] = %5.3f kg/m3\n", roh_air[0]);
+ printf("roh_air_desert0 = %5.3f kg/m3\n", roh_air_desert);
+ printf("zom_0 = %5.3f\n", zom0);
+ printf("u_0 = %5.3f\n", u_0);
+ printf("rah[0] = %5.3f s/m\n", rah[0]);
+ printf("h[0] = %5.3f W/m2\n", h[0]);
+ }
+
+/*----------------------------------------------------------------*/
+ /*Main iteration loop of SEBAL */
+ zom[0] = zom0;
+ for (ic = 1; ic < iteration + 1; ic++) {
+ if (debug == 1) {
+ printf("\n ******** ITERATION %i *********\n", ic);
+ }
+ /* Where is roh_air[i]? */
+ psih = psi_h(t0_dem, h[ic - 1], u_0, roh_air[ic - 1]);
+ ustar[ic] =
+ u_star(t0_dem, h[ic - 1], u_0, roh_air[ic - 1], zom[0], u2m);
+ rah[ic] = rah1(psih, ustar[ic]);
+ /* get desert point values from maps */
+ if (ic == 1) {
+ h_desert = rnet_desert - g0_desert;
+ zom_desert = 0.002;
+ psih_desert = psi_h(t0_dem_desert, h_desert, u_0, roh_air_desert);
+ ustar_desert =
+ u_star(t0_dem_desert, h_desert, u_0, roh_air_desert,
+ zom_desert, u2m);
+ }
+ else {
+ roh_air_desert = rohair(dem_desert, tempk_desert, dtair_desert);
+ h_desert = h1(roh_air_desert, rah_desert, dtair_desert);
+ ustar_desertold = ustar_desert;
+ psih_desert =
+ psi_h(t0_dem_desert, h_desert, ustar_desertold,
+ roh_air_desert);
+ ustar_desert =
+ u_star(t0_dem_desert, h_desert, ustar_desertold,
+ roh_air_desert, zom_desert, u2m);
+ }
+ rah_desert = rah1(psih_desert, ustar_desert);
+ dtair_desert = dt_air_desert(h_desert, roh_air_desert, rah_desert);
+ /* This should find the new dtair from inversed h equation... */
+ dtair[ic] = dt_air(t0_dem, tempk_water, tempk_desert, dtair_desert);
+ /* This produces h[ic] and roh_air[ic+1] */
+ roh_air[ic] = rohair(dem, tempk, dtair[ic]);
+ h[ic] = h1(roh_air[ic], rah[ic], dtair[ic]);
+ /* Output values of the iteration parameters */
+ if (debug == 1) {
+ printf("psih[%i] = %5.3f\n", ic, psih);
+ printf("ustar[%i] = %5.3f\n", ic, ustar[ic]);
+ printf("rah[%i] = %5.3f s/m\n", ic, rah[ic]);
+ printf("h_desert = %5.3f\n", h_desert);
+ printf("rohair_desert = %5.3f\n", roh_air_desert);
+ printf
+ ("psih_desert = %5.3f\tustar_desert = %5.3f\trah_desert = %5.3f\n",
+ psih_desert, ustar_desert, rah_desert);
+ printf("dtair_desert = %8.5f\n", dtair_desert);
+ printf("dtair[%i] = %5.3f K\n", ic, dtair[ic]);
+ printf("roh_air[%i] = %5.3f kg/m3\n", ic, roh_air[ic]);
+ printf("h[%i] = %5.3f W/m2\n", ic, h[ic]);
+ }
+ }
+ return h[iteration];
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/u_star.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/u_star.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/u_star.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,50 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double u_star(double t0_dem, double h, double ustar, double roh_air,
+ double zom, double u2m)
+{
+ double result;
+
+ double n5_temp; /* Monin-Obukov Length */
+
+ double n10_mem; /* psi m */
+
+ double n31_mem; /* x for U200 (that is bh...) */
+
+ double hv = 0.15; /* crop height (m) */
+
+ double bh = 200; /* blending height (m) */
+
+ if (h != 0.0) {
+ n5_temp =
+ (-1004 * roh_air * pow(ustar, 3) * t0_dem) / (0.41 * 9.81 * h);
+ }
+ else {
+ n5_temp = -1000.0;
+ }
+
+ if (n5_temp < 0.0) {
+
+ n31_mem = pow((1 - 16 * (200 / n5_temp)), 0.25);
+ n10_mem =
+ (2 * log10((1 + n31_mem) / 2) + log10((1 + pow(n31_mem, 2)) / 2) -
+ 2 * atan(n31_mem) + 0.5 * PI);
+
+ }
+ else {
+
+ /* n31_mem = 1.0; */
+ n10_mem = -5 * 2 / n5_temp;
+
+ }
+
+ result =
+ ((u2m * 0.41 / log10(2 / (hv / 7))) / 0.41 *
+ log10(bh / (hv / 7) * 0.41)) / (log10(bh / zom) - n10_mem);
+
+ return (result);
+}
Added: grass-addons/grass7/imagery/i.eb.h_SEBAL95/zom_0.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.h_SEBAL95/zom_0.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.h_SEBAL95/zom_0.c 2012-03-15 09:48:13 UTC (rev 51066)
@@ -0,0 +1,20 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double zom_0(double ndvi, double ndvi_max)
+{
+ double a, b, zom;
+
+ double hv_ndvimax = 1.5; /* crop vegetation height (m) */
+
+ double hv_desert = 0.002; /* desert base vegetation height (m) */
+
+ a = (log(hv_desert) -
+ ((log(hv_ndvimax / 7) - log(hv_desert)) / (ndvi_max - 0.02) * 0.02));
+ b = ((log(hv_ndvimax / 7) - log(hv_desert)) / (ndvi_max - 0.02)) * ndvi;
+ zom = exp(a + b);
+ /* Greece works and SEBAL01 */
+ /*zom = exp(3.3219*ndvi-3.9939); */
+ return (zom);
+}
More information about the grass-commit
mailing list