[GRASS-SVN] r56873 - in grass-addons/grass7/imagery: . i.eb.z0m
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 22 01:10:37 PDT 2013
Author: ychemin
Date: 2013-06-22 01:10:37 -0700 (Sat, 22 Jun 2013)
New Revision: 56873
Added:
grass-addons/grass7/imagery/i.eb.z0m/
grass-addons/grass7/imagery/i.eb.z0m/Makefile
grass-addons/grass7/imagery/i.eb.z0m/i.eb.z0m.html
grass-addons/grass7/imagery/i.eb.z0m/main.c
grass-addons/grass7/imagery/i.eb.z0m/z0m.c
Modified:
grass-addons/grass7/imagery/Makefile
Log:
upgraded to G7, merged i.eb.z0m and i.eb.z0m0
Modified: grass-addons/grass7/imagery/Makefile
===================================================================
--- grass-addons/grass7/imagery/Makefile 2013-06-22 04:38:21 UTC (rev 56872)
+++ grass-addons/grass7/imagery/Makefile 2013-06-22 08:10:37 UTC (rev 56873)
@@ -4,6 +4,7 @@
SUBDIRS = \
i.eb.h_sebal95 \
+ i.eb.z0m \
i.edge \
i.evapo.potrad \
i.evapo.senay \
Added: grass-addons/grass7/imagery/i.eb.z0m/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.eb.z0m/Makefile (rev 0)
+++ grass-addons/grass7/imagery/i.eb.z0m/Makefile 2013-06-22 08:10:37 UTC (rev 56873)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.eb.z0m
+
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/imagery/i.eb.z0m/i.eb.z0m.html
===================================================================
--- grass-addons/grass7/imagery/i.eb.z0m/i.eb.z0m.html (rev 0)
+++ grass-addons/grass7/imagery/i.eb.z0m/i.eb.z0m.html 2013-06-22 08:10:37 UTC (rev 56873)
@@ -0,0 +1,28 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.eb.z0m</em> calculates the momentum roughness length (z0m) and optionally the surface roughness for heat transport (z0h) as per SEBAL requirements from bastiaanssen (1995).
+Default: calculating from a NDVI with an deterministic equation, as seen in Bastiaanssen (1995).
+Flag -p : calculating from a SAVI with an empirical equation, as seen in Pawan (2004).
+This is a typical input to sensible heat flux computations of any energy balance modeling.
+<h2>NOTES</h2>
+
+The NDVI map input and the ndvi_max operation set, is only to get a linear relationship from NDVI to vegetation height. The latter being related to z0m by a factor 7.
+
+If you happen to have a vegetation height (hv) map, then z0m=hv/7 and z0h=0.1*z0m.
+<h2>TODO</h2>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.eb.h0.html">i.eb.h0</a><br>
+</em>
+
+
+<h2>AUTHORS</h2>
+Yann Chemin
+<br>
+
+
+<p>
+<i>Last changed: 2013-06-29</i>
Added: grass-addons/grass7/imagery/i.eb.z0m/main.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.z0m/main.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.z0m/main.c 2013-06-22 08:10:37 UTC (rev 56873)
@@ -0,0 +1,236 @@
+
+/****************************************************************************
+ *
+ * MODULE: i.eb.z0m
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculates the momentum roughness length
+ * as seen in Bastiaanssen (1995), Pawan (2004)
+ *
+ * COPYRIGHT: (C) 2002-2013 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 <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+double z0m_pawan(double sa_vi);
+double z0m_bastiaanssen(double ndvi, double ndvi_max, double hv_ndvimax,
+ double hv_desert);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /*region+header info */
+ char *mapset; /*mapset name */
+ int nrows, ncols;
+ int row, col;
+ int heat = 0; /*Flag for surf. roughness for heat transport output */
+ struct GModule *module;
+ struct Option *input1, *input2, *input3, *input4;
+ struct Option *output1, *output2;
+ struct Flag *flag1, *flag2;
+ struct History history; /*metadata */
+ /************************************/
+ char *name; /*input raster name */
+ char *result1, *result2; /*output raster name */
+ /*File Descriptors */
+ int infd_ndvi;
+ int outfd1, outfd2;
+ char *ndvi;
+ char *z0m, *z0h;
+ double coef_z0h, hv_max, hmr;
+ int i = 0, j = 0;
+ void *inrast_ndvi;
+ DCELL * outrast1, *outrast2;
+ RASTER_MAP_TYPE data_type_ndvi;
+ /************************************/
+ G_gisinit(argv[0]);
+ module = G_define_module();
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("roughness length"));
+ G_add_keyword(_("energy balance"));
+ G_add_keyword(_("SEBAL"));
+ module->description =
+ _("Momentum roughness length (z0m) and surface roughness for heat transport (z0h) as seen in Bastiaanssen (2004)");
+
+ /* Define the different options */
+ input1 = G_define_standard_option(G_OPT_R_INPUT);
+ input1->key = _("ndvi");
+ input1->description = _("Name of the NDVI map [-1.0;1.0], SAVI if -p flag");
+
+ input2 = G_define_option();
+ input2->key = _("coef");
+ input2->type = TYPE_DOUBLE;
+ input2->required = NO;
+ input2->gisprompt = _("parameter,value");
+ input2->description =
+ _("Value of the conversion factor from z0m and z0h (Bastiaanssen (2005) used 0.1), for -f flag");
+ input2->answer = _("0.1");
+
+ input3 = G_define_option();
+ input3->key = _("hv_max");
+ input3->type = TYPE_DOUBLE;
+ input3->required = NO;
+ input3->gisprompt = _("parameter,value");
+ input3->description =
+ _("Value of the vegetation height at max(NDVI) i.e. standard C3 crop could be 1.5m, not needed for -p flag");
+ input3->answer = _("1.5");
+
+ input4 = G_define_option();
+ input4->key = _("hmr");
+ input4->type = TYPE_DOUBLE;
+ input4->required = NO;
+ input4->gisprompt = _("parameter,value");
+ input4->description =
+ _("Value of the micro-relief height (h.m-r.) on flat bare ground, most references point to 2cm, not need for -p flag");
+ input4->answer = _("0.02");
+
+ output1 = G_define_standard_option(G_OPT_R_OUTPUT);
+ output1->description = _("Name of the output z0m layer");
+
+ output2 = G_define_standard_option(G_OPT_R_OUTPUT);
+ output2->key = _("z0h");
+ output2->required = NO;
+ output2->description = _("Name of the output z0h layer");
+
+ flag2 = G_define_flag();
+ flag2->key = 'p';
+ flag2->description = _("Use SAVI input only with equation of Pawan (2004)");
+
+ /********************/
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ ndvi = input1->answer;
+ if (input2->answer)
+ coef_z0h = atof(input2->answer);
+ if (input3->answer)
+ hv_max = atof(input3->answer);
+ if (input4->answer)
+ hmr = atof(input4->answer);
+ result1 = output1->answer;
+ if (output2->answer)
+ result2 = output2->answer;
+ /***************************************************/
+ data_type_ndvi = Rast_map_type(ndvi, "");
+ infd_ndvi = Rast_open_old(ndvi, "");
+ Rast_get_cellhd(ndvi, "", &cellhd);
+ inrast_ndvi = Rast_allocate_buf(data_type_ndvi);
+ /***************************************************/
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ /***************************************************/
+ outrast1 = Rast_allocate_d_buf();
+ Rast_open_new(result1, DCELL_TYPE);
+ if (input2->answer && output2->answer){
+ outrast2 = Rast_allocate_d_buf();
+ Rast_open_new(result2, DCELL_TYPE);
+ }
+ DCELL d_ndvi; /* Input raster */
+ DCELL d_ndvi_max = 0.0; /* Generated here */
+
+ /* NDVI Max */
+ if(flag2->answer){
+ /* skip it */
+ } else {
+ for (row = 0; row < nrows; row++)
+ {
+ G_percent(row, nrows, 2);
+ Rast_get_row(infd_ndvi,inrast_ndvi,row,data_type_ndvi);
+ 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.98) {
+ d_ndvi_max = d_ndvi;
+ }
+ }
+ }
+ G_message("ndvi_max=%f\n", d_ndvi_max);
+ }
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d;
+ DCELL d_z0h;
+ G_percent(row, nrows, 2);
+ /* read input maps */
+ Rast_get_row(infd_ndvi,inrast_ndvi,row,data_type_ndvi);
+ /*process the data */
+ 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)) {
+ Rast_set_d_null_value(&outrast1[col], 1);
+ if (input2->answer && output2->answer) {
+ Rast_set_d_null_value(&outrast2[col], 1);
+ }
+ }
+ else {
+ /****************************/
+ /* calculate z0m */
+ if(flag2->answer){
+ d = z0m_pawan(d_ndvi);
+ } else {
+ d = z0m_bastiaanssen(d_ndvi, d_ndvi_max, hv_max, hmr);
+ }
+ outrast1[col] = d;
+ if (input2->answer && output2->answer) {
+ d_z0h = d * coef_z0h;
+ outrast2[col] = d_z0h;
+ }
+ }
+ }
+ Rast_put_d_row(outfd1, outrast1);
+ if (input2->answer && output2->answer) {
+ Rast_put_d_row(outfd2, outrast2);
+ }
+ }
+ G_free(inrast_ndvi);
+ Rast_close(infd_ndvi);
+ G_free(outrast1);
+ Rast_close(outfd1);
+ if (input2->answer && output2->answer) {
+ G_free(outrast2);
+ Rast_close(outfd2);
+ }
+ Rast_short_history(result1, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result1, &history);
+ if (input2->answer && output2->answer) {
+ Rast_short_history(result2, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result2, &history);
+ }
+ exit(EXIT_SUCCESS);
+}
+
+
Added: grass-addons/grass7/imagery/i.eb.z0m/z0m.c
===================================================================
--- grass-addons/grass7/imagery/i.eb.z0m/z0m.c (rev 0)
+++ grass-addons/grass7/imagery/i.eb.z0m/z0m.c 2013-06-22 08:10:37 UTC (rev 56873)
@@ -0,0 +1,21 @@
+#include<stdio.h>
+#include<math.h>
+
+double z0m_bastiaanssen(double ndvi, double ndvi_max, double hv_ndvimax,
+ double hv_desert)
+{
+ double a, b, z0m;
+ /* hv_ndvimax = crop vegetation height (m) */
+ /* 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;
+ z0m = exp(a + b);
+ return (z0m);
+}
+
+/* Momentum roughness length (z0m) as seen in Pawan (2004) */
+double z0m_pawan(double sa_vi)
+{
+ return(exp(-5.809 + 5.62 * sa_vi));
+}
More information about the grass-commit
mailing list