[GRASS-SVN] r44972 - in grass/trunk/imagery: . i.evapo.PM
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 11 21:55:04 EST 2011
Author: ychemin
Date: 2011-01-11 18:55:03 -0800 (Tue, 11 Jan 2011)
New Revision: 44972
Added:
grass/trunk/imagery/i.evapo.PM/
grass/trunk/imagery/i.evapo.PM/i.evapo.PM.html
Removed:
grass/trunk/imagery/i.evapo.PM/description.html
Modified:
grass/trunk/imagery/i.evapo.PM/Makefile
grass/trunk/imagery/i.evapo.PM/functions.c
grass/trunk/imagery/i.evapo.PM/local_proto.h
grass/trunk/imagery/i.evapo.PM/main.c
Log:
Added i.evapo.PM
Modified: grass/trunk/imagery/i.evapo.PM/Makefile
===================================================================
--- grass-addons/HydroFOSS/r.evapo.PM/Makefile 2011-01-11 18:37:06 UTC (rev 44966)
+++ grass/trunk/imagery/i.evapo.PM/Makefile 2011-01-12 02:55:03 UTC (rev 44972)
@@ -1,10 +0,0 @@
-MODULE_TOPDIR = ../../..
-
-PGM = r.evapo.PM
-
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
-
-include $(MODULE_TOPDIR)/include/Make/Module.make
-
-default: cmd
Deleted: grass/trunk/imagery/i.evapo.PM/description.html
===================================================================
--- grass-addons/HydroFOSS/r.evapo.PM/description.html 2011-01-11 18:37:06 UTC (rev 44966)
+++ grass/trunk/imagery/i.evapo.PM/description.html 2011-01-12 02:55:03 UTC (rev 44972)
@@ -1,167 +0,0 @@
-<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
-<html>
-<head>
-<title>r.evapo.PM</title>
-<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
-<link rel="stylesheet" href="grassdocs.css" type="text/css">
-</head>
-<body bgcolor="white">
-
-<img src="grass.smlogo.gif" alt="GRASS logo"><hr align=center size=6 noshade>
-
-<H2>NAME</H2> <B><I>r.evapo.PM </I></B>- computation of <i>potential evapotranspiration</i> [mm/h] for hourly time step.
-
-<P><I>(GRASS Raster Program)</I>
-
-<H2>SYNOPSIS</H2>
-<B>r.evapo.PM</B>
-<BR>
-<B>r.evapo.PM</B> help</br>
-<BR>
-
-<B>r.evapo.PM</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>r.evapo.PM</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 potential evapotranspiration map (EPo).
-
-<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>0 vegetation is present and evapotranspiration is calculated;<br>
-- where Vh=0 bare ground is present and evapotranspiration is calculated;<br>
-- where Vh<0 water surface is present and evaporation is calculated;<br>
-
-<p>For more details on the algorithms see [1,2,3].
-
-
-<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>r.evapo.PM</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>r.evapo.PM</EM> on the
-command line and the program will ask for parameter values and flag
-settings interactively, using the standard GRASS parser interface.
-
-
-<H3>Flags:</H3>
-<dl>
- <dt><B>-q</B>
- <dd>Run quietly (do not display status messages). By default
- <EM>r.evapo.PM</EM> is run verbosely.
- <dt><B>-z</B>
- <dd>Set negative calculated evapotranspiration values to zero.
- <dt><B>-n</B>
- <dd>Calculate soil heat flux for night time. By default
- <EM>r.evapo.PM</EM> calculate it for day time.
-</dl>
-
-
-<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->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>r.evapo.PM</EM> command.
-<dt>example:
-<br><dd>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><dd>r.mapcalc 'NSR=0.0036*(beam+diffuse+reflected)';
-
-
-<H2>SEE ALSO</H2>
-<ul>
- <li>The <a href="http://istgis.ist.supsi.ch:8001/geomatica/">HydroFOSS</a>
-project at IST-SUPSI (Institute of Earth Sciences - University school of applied science for the Southern Switzerland)
-
- <li><a href=r.sun.html>r.sun</a>,
- <a href=r.mapcalc.html>r.mapcal</a>
-</ul>
-
-
-
-<H2>AUTHORS</H2>
-
- <p>Original version of program: The <a href="http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1">HydroFOSS</a> project, 2006, IST-SUPSI. (http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1)
- <i>
- <br>Massimiliano Cannata, Scuola Universitaria Professionale della Svizzera Italiana - Istituto Scienze della Terra
- <br>Maria A. Brovelli, Politecnico di Milano - Polo regionale di Como
- </i>
-
- <p>Contact: <a href="mailto:massimiliano.cannata at supsi.ch"> Massimiliano Cannata</a>
-
-
-<H2>REFERENCES</H2>
-
- <p>[1] Cannata M., 2006. <A HREF="http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1">
- GIS embedded approach for Free & Open Source Hydrological Modelling</A>. PhD thesis, Department of Geodesy and Geomatics, Polytechnic of Milan, Italy.
-
- <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.
-
-<p><i>Last changed: $Date$</i>
-</body>
-</html>
Modified: grass/trunk/imagery/i.evapo.PM/functions.c
===================================================================
--- grass-addons/HydroFOSS/r.evapo.PM/functions.c 2011-01-11 18:37:06 UTC (rev 44966)
+++ grass/trunk/imagery/i.evapo.PM/functions.c 2011-01-12 02:55:03 UTC (rev 44972)
@@ -5,7 +5,7 @@
#include <math.h>
#include "local_proto.h"
-extern CELL f_c(CELL);
+extern CELL f_c(CELL);
extern FCELL f_f(FCELL);
extern DCELL f_d(DCELL);
@@ -13,148 +13,153 @@
//#define FALSE 1
/* constant definition */
-//#define k_sb 4.903 //[MJ/m2*h] Stefan Bolzman constant
-#define cp 1.013 //[kJ/kg*°C] specific heat of moist air
-#define epsilon 0.622 //[-] ratio of molecular weigth of water to dry air
-#define Po 101.3 //[kPa] atmospheric pressure at sea level
-#define Tko 293.16 //[K] reference temperature at sea level
-#define eta 0.0065 //[K/m] constant lapse rate
-#define Ao 0 //[m] altitude at sea level
-#define g 9.81 //[m/s] gravitational accelleration
-#define R 287 //[J/kg*K] specific gas constant
-#define Zw 2 //[m] height of wind measurements
-#define Zh 2 //[m] height of humidity measurements
-#define k 0.41 //[-] Von Karman constant
+//#define k_sb 4.903 //[MJ/m2*h] Stefan Bolzman constant
+#define cp 1.013 //[kJ/kg*°C] specific heat of moist air
+#define epsilon 0.622 //[-] ratio of molecular weigth of water to dry air
+#define Po 101.3 //[kPa] atmospheric pressure at sea level
+#define Tko 293.16 //[K] reference temperature at sea level
+#define eta 0.0065 //[K/m] constant lapse rate
+#define Ao 0 //[m] altitude at sea level
+#define g 9.81 //[m/s] gravitational accelleration
+#define R 287 //[J/kg*K] specific gas constant
+#define Zw 2 //[m] height of wind measurements
+#define Zh 2 //[m] height of humidity measurements
+#define k 0.41 //[-] Von Karman constant
-DCELL calc_ETp(DCELL T, DCELL Z, DCELL u2, DCELL Rn, int night, DCELL Rh, DCELL hc){
-
- DCELL ea, delta, gamma, gstar, lambda;
- DCELL P, ra, d, Zom, Zoh, G, ETrad, u10, rs, ed, Tkv, rho, ETaero, ETp;
-
- /* calculus: mean saturation vapoure pressure [KPa] */
- ea = 0.61078*exp((17.27*T)/(T+237.3));
-
- /* calculus: slope of vapoure pressure curve [KPa/°C] */
- delta = (4098*ea)/pow((237.3+T),2);
-
- /* calculus: latent heat vapourisation [MJ/kg] */
- lambda = 2.501 - (0.002361*T);
-
- /* calculus: atmospheric pressure [KPa] */
- P = Po * pow(((Tko-eta*(Z-Ao))/Tko),(g/(eta*R)));
-
- /* calculus: psichiometric constant [kPa/°C] */
- gamma = ((cp*P)/(epsilon*lambda))*0.001;
-
- /* calculus: aerodynamic resistance [s/m] */
- if (hc<2) {
- d = (2/3)*hc;
- Zom = 0.123*hc;
- Zoh = 0.1*Zom;
- ra = ( log((Zw-d)/Zom) * log((Zh-d)/Zoh) ) / (k*k*u2);
- }
- else {
- u10 = u2*(log((67.8*10)-5.42))/4.87;
- ra = 94 / u10;
- }
-
- /* calculus: surface resistance [s/m] */
- rs = 100/(0.5*24*hc);
-
- /*calculus: modified psichiometric constant [kPa/°C] */
- gstar = gamma*(1+(rs/ra));
+DCELL calc_ETp(DCELL T, DCELL Z, DCELL u2, DCELL Rn, int night, DCELL Rh,
+ DCELL hc)
+{
- /*calculus: net radiation [MJ/m2*d] */
- /*Rn derived from r.sun */
-
- /*calculus: soil heat flux [MJ/m2*d] */
- if (night==FALSE)
- G=0.1*Rn;
- else
- G=0.5*Rn;
-
- /* calculus: radiation term [mm/h] */
- /* ETrad = (delta/(delta+gstar))*((Rn-G)/(lambda*1000000)); */
- ETrad = (delta/(delta+gstar))*((Rn-G)/lambda); /* torna da analisi dimensionale */
-
- /* calculus: actual vapoure pressure [kPa] */
- ed = Rh*ea/100;
-
- /* calculus: virtual temperature [°C] */
- Tkv = (T+273.15)/(1-(0.378*ed/P));
-
- /* calculus: atmospheric density [Kg/m^3] */
- rho = P/(Tkv*R/100);
-
- /* calculus: aerodynamic term [mm/h] */
- /* ETaero = (0.001/lambda)*(1/(delta+gstar))*(rho*cp/ra)*(ea-ed); */
- ETaero = (3.6/lambda)*(1/(delta+gstar))*(rho*cp/ra)*(ea-ed); /* torna da analisi dimensionale */
-
- /* calculus: potential evapotranspiration [mm/h] */
- ETp = ETrad + ETaero;
+ DCELL ea, delta, gamma, gstar, lambda;
+ DCELL P, ra, d, Zom, Zoh, G, ETrad, u10, rs, ed, Tkv, rho, ETaero, ETp;
- return ETp;
-
+ /* calculus: mean saturation vapoure pressure [KPa] */
+ ea = 0.61078 * exp((17.27 * T) / (T + 237.3));
+
+ /* calculus: slope of vapoure pressure curve [KPa/°C] */
+ delta = (4098 * ea) / pow((237.3 + T), 2);
+
+ /* calculus: latent heat vapourisation [MJ/kg] */
+ lambda = 2.501 - (0.002361 * T);
+
+ /* calculus: atmospheric pressure [KPa] */
+ P = Po * pow(((Tko - eta * (Z - Ao)) / Tko), (g / (eta * R)));
+
+ /* calculus: psichiometric constant [kPa/°C] */
+ gamma = ((cp * P) / (epsilon * lambda)) * 0.001;
+
+ /* calculus: aerodynamic resistance [s/m] */
+ if (hc < 2) {
+ d = (2 / 3) * hc;
+ Zom = 0.123 * hc;
+ Zoh = 0.1 * Zom;
+ ra = (log((Zw - d) / Zom) * log((Zh - d) / Zoh)) / (k * k * u2);
+ }
+ else {
+ u10 = u2 * (log((67.8 * 10) - 5.42)) / 4.87;
+ ra = 94 / u10;
+ }
+
+ /* calculus: surface resistance [s/m] */
+ rs = 100 / (0.5 * 24 * hc);
+
+ /*calculus: modified psichiometric constant [kPa/°C] */
+ gstar = gamma * (1 + (rs / ra));
+
+ /*calculus: net radiation [MJ/m2*d] */
+ /*Rn derived from r.sun */
+
+ /*calculus: soil heat flux [MJ/m2*d] */
+ if (night == FALSE)
+ G = 0.1 * Rn;
+ else
+ G = 0.5 * Rn;
+
+ /* calculus: radiation term [mm/h] */
+ /* ETrad = (delta/(delta+gstar))*((Rn-G)/(lambda*1000000)); */
+ ETrad = (delta / (delta + gstar)) * ((Rn - G) / lambda); /* torna da analisi dimensionale */
+
+ /* calculus: actual vapoure pressure [kPa] */
+ ed = Rh * ea / 100;
+
+ /* calculus: virtual temperature [°C] */
+ Tkv = (T + 273.15) / (1 - (0.378 * ed / P));
+
+ /* calculus: atmospheric density [Kg/m^3] */
+ rho = P / (Tkv * R / 100);
+
+ /* calculus: aerodynamic term [mm/h] */
+ /* ETaero = (0.001/lambda)*(1/(delta+gstar))*(rho*cp/ra)*(ea-ed); */
+ ETaero = (3.6 / lambda) * (1 / (delta + gstar)) * (rho * cp / ra) * (ea - ed); /* torna da analisi dimensionale */
+
+ /* calculus: potential evapotranspiration [mm/h] */
+ ETp = ETrad + ETaero;
+
+ return ETp;
+
}
-DCELL calc_openwaterETp(DCELL T, DCELL Z, DCELL u2, DCELL Rn, int day, DCELL Rh, DCELL hc){
- DCELL ea, delta, gamma, lambda;
- DCELL P, ed, ETaero, ETp;
-
- /* calculus: mean saturation vapoure pressure [KPa] */
- ea = 0.61078*exp((17.27*T)/(T+237.3));
-
- /* calculus: slope of vapoure pressure curve [KPa/°C] */
- delta = (4098*ea)/pow((237.3+T),2);
-
- /* calculus: latent heat vapourisation [MJ/kg] */
- lambda = 2.501 - (0.002361*T);
-
- /* calculus: atmospheric pressure [KPa] */
- P = Po * pow(((Tko-eta*(Z-Ao))/Tko),(g/(eta*R)));
-
- /* calculus: di psichiometric constant [kPa/°C] */
- gamma = ((cp*P)/(epsilon*lambda))*0.001;
-
- /*calculus: net radiation [MJ/m2*h] */
- /*Rn derived from r.sun
-
- /*calculus: actual vapoure pressure [kPa] */
- ed = Rh*ea/100;
-
- /*calculus: aerodynamic term [mm/h] */
- /*ETaero = 0.35*(0.5+(0.621375*u2/100))*7.500638*(ea-ed); */
- /*to convert mm/d to mm/h it results: */
- ETaero = (0.35/24)*(0.5+(0.621375*u2/100))*7.500638*(ea-ed);
-
- /*calculus: potential evapotranspiration [mm/h] */
- ETp = (((Rn*delta)/lambda)+(gamma*ETaero))/(delta+gamma);
+DCELL calc_openwaterETp(DCELL T, DCELL Z, DCELL u2, DCELL Rn, int day,
+ DCELL Rh, DCELL hc)
+{
+ DCELL ea, delta, gamma, lambda;
+ DCELL P, ed, ETaero, ETp;
- return ETp;
-
+ /* calculus: mean saturation vapoure pressure [KPa] */
+ ea = 0.61078 * exp((17.27 * T) / (T + 237.3));
+
+ /* calculus: slope of vapoure pressure curve [KPa/°C] */
+ delta = (4098 * ea) / pow((237.3 + T), 2);
+
+ /* calculus: latent heat vapourisation [MJ/kg] */
+ lambda = 2.501 - (0.002361 * T);
+
+ /* calculus: atmospheric pressure [KPa] */
+ P = Po * pow(((Tko - eta * (Z - Ao)) / Tko), (g / (eta * R)));
+
+ /* calculus: di psichiometric constant [kPa/°C] */
+ gamma = ((cp * P) / (epsilon * lambda)) * 0.001;
+
+ /*calculus: net radiation [MJ/m2*h] */
+ /*Rn derived from r.sun
+
+ /*calculus: actual vapoure pressure [kPa] */
+ ed = Rh * ea / 100;
+
+ /*calculus: aerodynamic term [mm/h] */
+ /*ETaero = 0.35*(0.5+(0.621375*u2/100))*7.500638*(ea-ed); */
+ /*to convert mm/d to mm/h it results: */
+ ETaero =
+ (0.35 / 24) * (0.5 + (0.621375 * u2 / 100)) * 7.500638 * (ea - ed);
+
+ /*calculus: potential evapotranspiration [mm/h] */
+ ETp = (((Rn * delta) / lambda) + (gamma * ETaero)) / (delta + gamma);
+
+ return ETp;
+
}
/*
-DCELL calc_ETp_WaSiM(){
-
+ DCELL calc_ETp_WaSiM(){
- //calculus of saturation vapour pressure at the temperature T: es[hPa]
- //with T[°C]
- es = 6.1078*exp((17.27*T)/(237.3+T);
- //tangent of the saturated vapour pressure curve: D[hPa/K] with T[°C]
- D = (25029/pow((273.3+T),2))*exp((17.27*T)/(237.3+T);
-
- //air pressure at level hM
- P = 1013*exp(-hM/(7991+29.33*Tv));
-
- //calculus of lambda [KJ/Kg] with T[°C]
- lambda = 2500.8 - 2.372*T;
-
- //calculus psichiometric constant gamma [hPa/°C] with cp=1.005[KJ/(Kg*K)]
- gamma = ((1.005*P)/(0.622*lambda))*0.001;
-
- //
-*/
+ //calculus of saturation vapour pressure at the temperature T: es[hPa]
+ //with T[°C]
+ es = 6.1078*exp((17.27*T)/(237.3+T);
+
+ //tangent of the saturated vapour pressure curve: D[hPa/K] with T[°C]
+ D = (25029/pow((273.3+T),2))*exp((17.27*T)/(237.3+T);
+
+ //air pressure at level hM
+ P = 1013*exp(-hM/(7991+29.33*Tv));
+
+ //calculus of lambda [KJ/Kg] with T[°C]
+ lambda = 2500.8 - 2.372*T;
+
+ //calculus psichiometric constant gamma [hPa/°C] with cp=1.005[KJ/(Kg*K)]
+ gamma = ((1.005*P)/(0.622*lambda))*0.001;
+
+ //
+ */
Copied: grass/trunk/imagery/i.evapo.PM/i.evapo.PM.html (from rev 44966, grass-addons/HydroFOSS/r.evapo.PM/description.html)
===================================================================
--- grass/trunk/imagery/i.evapo.PM/i.evapo.PM.html (rev 0)
+++ grass/trunk/imagery/i.evapo.PM/i.evapo.PM.html 2011-01-12 02:55:03 UTC (rev 44972)
@@ -0,0 +1,72 @@
+<H2>DESCRIPTION</H2>
+
+<p><EM>i.evapo.PM</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 potential evapotranspiration map (EPo).
+
+<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 gt 0 vegetation is present and evapotranspiration is calculated;<br>
+- where Vh = 0 bare ground is present and evapotranspiration is calculated;<br>
+- where Vh lt 0 water surface is present and evaporation is calculated;<br>
+
+<p>For more details on the algorithms see [1,2,3].
+
+<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>
+<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>The <a href="http://istgis.ist.supsi.ch:8001/geomatica/">HydroFOSS</a>
+project at IST-SUPSI (Institute of Earth Sciences - University school of applied science for the Southern Switzerland)
+
+ <li><a href=r.sun.html>r.sun</a>,
+ <a href=r.mapcalc.html>r.mapcalc</a>
+</ul>
+
+
+
+<H2>AUTHORS</H2>
+
+ <p>Original version of program: The <a href="http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1">HydroFOSS</a> project, 2006, IST-SUPSI. (http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1)
+ <i>
+ <br>Massimiliano Cannata, Scuola Universitaria Professionale della Svizzera Italiana - Istituto Scienze della Terra
+ <br>Maria A. Brovelli, Politecnico di Milano - Polo regionale di Como
+ </i>
+
+ <p>Contact: <a href="mailto:massimiliano.cannata at supsi.ch"> Massimiliano Cannata</a>
+
+
+<H2>REFERENCES</H2>
+
+ <p>[1] Cannata M., 2006. <A HREF="http://istgis.ist.supsi.ch:8001/geomatica/index.php?id=1">
+ GIS embedded approach for Free & Open Source Hydrological Modelling</A>. PhD thesis, Department of Geodesy and Geomatics, Polytechnic of Milan, Italy.
+
+ <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.
+
+<p><i>Last changed: $Date$</i>
Modified: grass/trunk/imagery/i.evapo.PM/local_proto.h
===================================================================
--- grass-addons/HydroFOSS/r.evapo.PM/local_proto.h 2011-01-11 18:37:06 UTC (rev 44966)
+++ grass/trunk/imagery/i.evapo.PM/local_proto.h 2011-01-12 02:55:03 UTC (rev 44972)
@@ -1,4 +1,5 @@
/* function.c */
+#include <grass/raster.h>
DCELL calc_Delta (DCELL T);
DCELL calc_g (DCELL Z);
DCELL calc_Eo (DCELL T);
Modified: grass/trunk/imagery/i.evapo.PM/main.c
===================================================================
--- grass-addons/HydroFOSS/r.evapo.PM/main.c 2011-01-11 18:37:06 UTC (rev 44966)
+++ grass/trunk/imagery/i.evapo.PM/main.c 2011-01-12 02:55:03 UTC (rev 44972)
@@ -1,266 +1,242 @@
/* Created by Anjuta version 1.0.2 */
-/* This file will not be overwritten */
+/****************************************************************************
+ *
+ * MODULE: i.evapo.PM
+ * AUTHOR(S): Massimiliano Cannata - massimiliano.cannata AT supsi.ch
+ * Maria A. Brovelli
+ * PURPOSE: Originally r.evapo.PM from HydroFOSS
+ * Calculates the Penman-Monteith reference evapotranspiration
+ * and Open Water Evaporation.
+ *
+ * COPYRIGHT: (C) 2006-2010 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.
+ *
+ *****************************************************************************/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
-#include <grass/gis.h>
#include <math.h>
-#include "local_proto.h"
+#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
+#include "local_proto.h"
-
int main(int argc, char *argv[])
-{
- struct Cell_head cellhd;
-
- /* buffer for in out raster */
- DCELL *inrast_T,*inrast_RH,*inrast_u2,*inrast_Rn,*inrast_DEM,*inrast_hc,*outrast;
- unsigned char *EPo;
-
- int nrows, ncols;
- int row, col;
- int infd_T,infd_RH,infd_u2,infd_Rn,infd_DEM,infd_hc;
- int outfd;
-
- char *mapset_T,*mapset_RH,*mapset_u2,*mapset_Rn,*mapset_DEM,*mapset_hc;
- char *T, *RH, *u2, *Rn, *DEM, *hc;
- DCELL d_T,d_RH,d_u2,d_Rn,d_Z,d_hc;
- DCELL d_EPo;
-
- int d_night;
-
- struct History history;
- struct GModule *module;
- struct Option *input_DEM, *input_T, *input_RH, *input_u2, *input_Rn, *input_hc, *output;
- struct Flag *flag1, *day, *zero;
-
-
- G_gisinit(argv[0]);
-
- module = G_define_module();
- module->description = _("Potontial Evapotranspiration Calculation with hourly Penman-Monteith");
+{
+ struct Cell_head cellhd;
- /* Define different options */
- input_DEM = G_define_option();
- input_DEM->key = "DEM";
- input_DEM->type = TYPE_STRING;
- input_DEM->required = YES;
- input_DEM->gisprompt = "old,cell,raster";
- input_DEM->description = _("Name of DEM raster map [m a.s.l.]");
-
- input_T = G_define_option();
- input_T->key = "T";
- input_T->type = TYPE_STRING;
- input_T->required = YES;
- input_T->gisprompt = "old,cell,raster";
- input_T->description = _("Name of Temperature raster map [°C]");
-
- input_RH = G_define_option();
- input_RH->key = "RU";
- input_RH->type = TYPE_STRING;
- input_RH->required = YES;
- input_RH->gisprompt = "old,cell,raster";
- input_RH->description = _("Name of Relative Umidity raster map [%]");
-
- input_u2 = G_define_option();
- input_u2->key = "WS";
- input_u2->type = TYPE_STRING;
- input_u2->required = YES;
- input_u2->gisprompt = "old,cell,raster";
- input_u2->description = _("Name of Wind Speed raster map [m/s]");
-
- input_Rn = G_define_option();
- input_Rn->key = "NSR";
- input_Rn->type = TYPE_STRING;
- input_Rn->required = YES;
- input_Rn->gisprompt = "old,cell,raster";
- input_Rn->description = _("Name of Net Solar Radiation raster map [MJ/m2/h]");
-
- input_hc = G_define_option();
- input_hc->key = "Vh";
- input_hc->type = TYPE_STRING;
- input_hc->required = YES;
- input_hc->gisprompt = "old,cell,raster";
- input_hc->description = _("Name of crop height raster map [m]");
-
- output = G_define_option() ;
- output->key = "EPo";
- output->type = TYPE_STRING;
- output->required = YES;
- output->gisprompt = "new,cell,raster" ;
- output->description= _("Name of output Reference Potential Evapotranspiration layer [mm/h]");
-
- /* Define the different flags */
- //flag1 = G_define_flag() ;
- //flag1->key = 'q' ;
- //flag1->description = "Quiet" ;
-
- zero = G_define_flag() ;
- zero->key = 'z' ;
- zero->description = _("set negative evapo to zero");
-
- day = G_define_flag() ;
- day->key = 'n' ;
- day->description = _("night-time");
-
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
-
- /* get entered parameters */
- T=input_T->answer;
- RH=input_RH->answer;
- u2=input_u2->answer;
- Rn=input_Rn->answer;
- EPo=output->answer;
- DEM=input_DEM->answer;
- hc=input_hc->answer;
-
- if (day->answer) {
- d_night = TRUE;
- }
- else {
- d_night=FALSE;
- }
-
- /* find maps in mapset */
- mapset_T = G_find_cell2 (T, "");
- if (mapset_T == NULL)
- G_fatal_error (_("cell file [%s] not found"), T);
- mapset_RH = G_find_cell2 (RH, "");
- if (mapset_RH == NULL)
- G_fatal_error (_("cell file [%s] not found"), RH);
- mapset_u2 = G_find_cell2 (u2, "");
- if (mapset_u2 == NULL)
- G_fatal_error (_("cell file [%s] not found"), u2);
- mapset_Rn = G_find_cell2 (Rn, "");
- if (mapset_Rn == NULL)
- G_fatal_error (_("cell file [%s] not found"), Rn);
- mapset_DEM = G_find_cell2 (DEM, "");
- if (mapset_DEM == NULL)
- G_fatal_error (_("cell file [%s] not found"), DEM);
- mapset_hc = G_find_cell2 (hc, "");
- if (mapset_hc == NULL)
- G_fatal_error (_("cell file [%s] not found"), hc);
+ /* buffer for in out raster */
+ DCELL *inrast_T, *inrast_RH, *inrast_u2;
+ DCELL *inrast_Rn, *inrast_DEM, *inrast_hc, *outrast;
+ unsigned char *EPo;
+ int nrows, ncols;
+ int row, col;
+ int infd_T, infd_RH, infd_u2, infd_Rn, infd_DEM, infd_hc;
+ int outfd;
-
- /* check legal output name */
- if (G_legal_filename (EPo) < 0)
- G_fatal_error (_("[%s] is an illegal name"), EPo);
-
- /* determine the input map type (CELL/FCELL/DCELL) */
- //data_type = G_raster_map_type(T, mapset);
+ char *mapset_T, *mapset_RH, *mapset_u2;
+ char *mapset_Rn, *mapset_DEM, *mapset_hc;
+ char *T, *RH, *u2, *Rn, *DEM, *hc;
+ DCELL d_T, d_RH, d_u2, d_Rn, d_Z, d_hc;
+ DCELL d_EPo;
- if ( (infd_T = G_open_cell_old (T, mapset_T)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"), T);
- if ( (infd_RH = G_open_cell_old (RH, mapset_RH)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"),RH);
- if ( (infd_u2 = G_open_cell_old (u2, mapset_u2)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"),u2);
- if ( (infd_Rn = G_open_cell_old (Rn, mapset_Rn)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"),Rn);
- if ( (infd_DEM = G_open_cell_old (DEM, mapset_DEM)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"),DEM);
- if ( (infd_hc = G_open_cell_old (hc, mapset_hc)) < 0)
- G_fatal_error (_("Cannot open cell file [%s]"),hc);
-
- if (G_get_cellhd (T, mapset_T, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), T);
- if (G_get_cellhd (RH, mapset_RH, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), RH);
- if (G_get_cellhd (u2, mapset_u2, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), u2);
- if (G_get_cellhd (Rn, mapset_Rn, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), Rn);
- if (G_get_cellhd (DEM, mapset_DEM, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), DEM);
- if (G_get_cellhd (hc, mapset_hc, &cellhd) < 0)
- G_fatal_error (_("Cannot read file header of [%s]"), hc);
+ int d_night;
- /* Allocate input buffer */
- inrast_T = G_allocate_d_raster_buf();
- inrast_RH = G_allocate_d_raster_buf();
- inrast_u2 = G_allocate_d_raster_buf();
- inrast_Rn = G_allocate_d_raster_buf();
- inrast_DEM = G_allocate_d_raster_buf();
- inrast_hc = G_allocate_d_raster_buf();
-
- /* Allocate output buffer */
- nrows = G_window_rows();
- ncols = G_window_cols();
- outrast = G_allocate_d_raster_buf();
+ struct History history;
+ struct GModule *module;
+ struct Option *input_DEM, *input_T, *input_RH;
+ struct Option *input_u2, *input_Rn, *input_hc, *output;
+ struct Flag *flag1, *day, *zero;
- if ( (outfd = G_open_raster_new (EPo,DCELL_TYPE)) < 0)
- G_fatal_error (_("Could not open <%s>"),T);
-
- for (row = 0; row < nrows; row++)
- {
-
- /* read a line input maps into buffers*/
- if (G_get_d_raster_row (infd_T, inrast_T, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),T);
- if (G_get_d_raster_row (infd_RH, inrast_RH, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),RH);
- if (G_get_d_raster_row (infd_u2, inrast_u2, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),u2);
- if (G_get_d_raster_row (infd_Rn, inrast_Rn, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),Rn);
- if (G_get_d_raster_row (infd_DEM, inrast_DEM, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),DEM);
- if (G_get_d_raster_row (infd_hc, inrast_hc, row) < 0)
- G_fatal_error (_("Could not read from <%s>"),hc);
-
- /* read every cell in the line buffers */
- for (col=0; col < ncols; col++)
- {
- d_T = ((DCELL *) inrast_T)[col];
- d_RH = ((DCELL *) inrast_RH)[col];
- d_u2 = ((DCELL *) inrast_u2)[col];
- d_Rn = ((DCELL *) inrast_Rn)[col];
- d_Z = ((DCELL *) inrast_DEM)[col];
- d_hc = ((DCELL *) inrast_hc)[col];
-
- //calculate evapotranspiration
- if (d_hc<0){
- //calculate evaporation
- d_EPo=calc_openwaterETp(d_T,d_Z,d_u2,d_Rn,d_night,d_RH,d_hc);
- }
- else {
- //calculate evapotranspiration
- d_EPo=calc_ETp(d_T,d_Z,d_u2,d_Rn,d_night,d_RH,d_hc);
- }
-
- if (zero->answer && d_EPo<0)
- d_EPo=0;
-
- ((DCELL *) outrast)[col] = d_EPo;
- }
-
- if (G_put_d_raster_row (outfd, outrast) < 0)
- G_fatal_error (_("Cannot write to <%s>"),EPo);
-
- }
- G_free(inrast_T);
- G_free(inrast_RH);
- G_free(inrast_u2);
- G_free(inrast_Rn);
- G_free(inrast_DEM);
- G_free(inrast_hc);
- G_free(outrast);
- G_close_cell (infd_T);
- G_close_cell (infd_RH);
- G_close_cell (infd_u2);
- G_close_cell (infd_Rn);
- G_close_cell (infd_DEM);
- G_close_cell (infd_hc);
- G_close_cell (outfd);
-
- /* add command line incantation to history file */
- G_short_history(output, "raster", &history);
- G_command_history(&history);
- G_write_history(output, &history);
+ G_gisinit(argv[0]);
- exit(EXIT_SUCCESS);
+ module = G_define_module();
+ module->description =
+ _("Potontial Evapotranspiration Calculation with hourly Penman-Monteith");
+
+ /* Define different options */
+ input_DEM = G_define_option();
+ input_DEM->key = "DEM";
+ input_DEM->type = TYPE_STRING;
+ input_DEM->required = YES;
+ input_DEM->gisprompt = "old,cell,raster";
+ input_DEM->description = _("Name of DEM raster map [m a.s.l.]");
+
+ input_T = G_define_option();
+ input_T->key = "T";
+ input_T->type = TYPE_STRING;
+ input_T->required = YES;
+ input_T->gisprompt = "old,cell,raster";
+ input_T->description = _("Name of Temperature raster map [°C]");
+
+ input_RH = G_define_option();
+ input_RH->key = "RU";
+ input_RH->type = TYPE_STRING;
+ input_RH->required = YES;
+ input_RH->gisprompt = "old,cell,raster";
+ input_RH->description = _("Name of Relative Umidity raster map [%]");
+
+ input_u2 = G_define_option();
+ input_u2->key = "WS";
+ input_u2->type = TYPE_STRING;
+ input_u2->required = YES;
+ input_u2->gisprompt = "old,cell,raster";
+ input_u2->description = _("Name of Wind Speed raster map [m/s]");
+
+ input_Rn = G_define_option();
+ input_Rn->key = "NSR";
+ input_Rn->type = TYPE_STRING;
+ input_Rn->required = YES;
+ input_Rn->gisprompt = "old,cell,raster";
+ input_Rn->description =
+ _("Name of Net Solar Radiation raster map [MJ/m2/h]");
+
+ input_hc = G_define_option();
+ input_hc->key = "Vh";
+ input_hc->type = TYPE_STRING;
+ input_hc->required = YES;
+ input_hc->gisprompt = "old,cell,raster";
+ input_hc->description = _("Name of crop height raster map [m]");
+
+ output = G_define_option();
+ output->key = "EPo";
+ output->type = TYPE_STRING;
+ output->required = YES;
+ output->gisprompt = "new,cell,raster";
+ output->description =
+ _("Name of output Reference Potential Evapotranspiration layer [mm/h]");
+
+ /* Define the different flags */
+ //flag1 = G_define_flag() ;
+ //flag1->key = 'q' ;
+ //flag1->description = "Quiet" ;
+
+ zero = G_define_flag();
+ zero->key = 'z';
+ zero->description = _("set negative evapo to zero");
+
+ day = G_define_flag();
+ day->key = 'n';
+ day->description = _("night-time");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get entered parameters */
+ T = input_T->answer;
+ RH = input_RH->answer;
+ u2 = input_u2->answer;
+ Rn = input_Rn->answer;
+ EPo = output->answer;
+ DEM = input_DEM->answer;
+ hc = input_hc->answer;
+
+ if (day->answer) {
+ d_night = TRUE;
+ }
+ else {
+ d_night = FALSE;
+ }
+
+ /* check legal output name */
+ if (G_legal_filename(EPo) < 0)
+ G_fatal_error(_("[%s] is an illegal name"), EPo);
+
+ /* determine the input map type (CELL/FCELL/DCELL) */
+ //data_type = G_raster_map_type(T, mapset);
+
+ infd_T = Rast_open_old(T, "");
+ infd_RH = Rast_open_old(RH, "");
+ infd_u2 = Rast_open_old(u2, "");
+ infd_Rn = Rast_open_old(Rn, "");
+ infd_DEM = Rast_open_old(DEM, "");
+ infd_hc = Rast_open_old(hc, "");
+
+ Rast_get_cellhd(T, "", &cellhd);
+ Rast_get_cellhd(RH, "", &cellhd);
+ Rast_get_cellhd(u2, "", &cellhd);
+ Rast_get_cellhd(Rn, "", &cellhd);
+ Rast_get_cellhd(DEM, "", &cellhd);
+ Rast_get_cellhd(hc, "", &cellhd);
+
+ /* Allocate input buffer */
+ inrast_T = Rast_allocate_d_buf();
+ inrast_RH = Rast_allocate_d_buf();
+ inrast_u2 = Rast_allocate_d_buf();
+ inrast_Rn = Rast_allocate_d_buf();
+ inrast_DEM = Rast_allocate_d_buf();
+ inrast_hc = Rast_allocate_d_buf();
+
+ /* Allocate output buffer */
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ outrast = Rast_allocate_d_buf();
+
+ outfd = Rast_open_new(EPo, DCELL_TYPE);
+
+ for (row = 0; row < nrows; row++) {
+
+ /* read a line input maps into buffers */
+ Rast_get_d_row(infd_T, inrast_T, row);
+ Rast_get_d_row(infd_RH, inrast_RH, row);
+ Rast_get_d_row(infd_u2, inrast_u2, row);
+ Rast_get_d_row(infd_Rn, inrast_Rn, row);
+ Rast_get_d_row(infd_DEM, inrast_DEM, row);
+ Rast_get_d_row(infd_hc, inrast_hc, row);
+
+ /* read every cell in the line buffers */
+ for (col = 0; col < ncols; col++) {
+ d_T = ((DCELL *) inrast_T)[col];
+ d_RH = ((DCELL *) inrast_RH)[col];
+ d_u2 = ((DCELL *) inrast_u2)[col];
+ d_Rn = ((DCELL *) inrast_Rn)[col];
+ d_Z = ((DCELL *) inrast_DEM)[col];
+ d_hc = ((DCELL *) inrast_hc)[col];
+
+ //calculate evapotranspiration
+ if (d_hc < 0) {
+ //calculate evaporation
+ d_EPo =
+ calc_openwaterETp(d_T, d_Z, d_u2, d_Rn, d_night, d_RH,
+ d_hc);
+ }
+ else {
+ //calculate evapotranspiration
+ d_EPo = calc_ETp(d_T, d_Z, d_u2, d_Rn, d_night, d_RH, d_hc);
+ }
+
+ if (zero->answer && d_EPo < 0)
+ d_EPo = 0;
+
+ ((DCELL *) outrast)[col] = d_EPo;
+ }
+ Rast_put_d_row(outfd, outrast);
+ }
+ G_free(inrast_T);
+ G_free(inrast_RH);
+ G_free(inrast_u2);
+ G_free(inrast_Rn);
+ G_free(inrast_DEM);
+ G_free(inrast_hc);
+ G_free(outrast);
+ Rast_close(infd_T);
+ Rast_close(infd_RH);
+ Rast_close(infd_u2);
+ Rast_close(infd_Rn);
+ Rast_close(infd_DEM);
+ Rast_close(infd_hc);
+ Rast_close(outfd);
+
+ /* add command line incantation to history file */
+ Rast_short_history(EPo, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(EPo, &history);
+
+ exit(EXIT_SUCCESS);
}
More information about the grass-commit
mailing list