added i.sunhours and i.emissivity

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@33939 15284696-431f-4ddb-bdfa-cd5b030d7da7

gui/tcltk/gis.m/gmmenu.tcl

@@ -550,10 +550,12 @@ set descmenu [subst  {
 	{cascad {[G_msg "Basic RS processing"]} {} "" $tmenu {			
 				{command {[G_msg "MODIS Quality flags"]} {} "i.modis.qc: Extracts Modis Quality flags as raster values." {} -command {execute i.modis.qc }}
 				{command {[G_msg "Latitude/longitude maps"]} {} "i.latlong: Latitude or Longitude maps." {} -command {execute i.latlong }}
+				{command {[G_msg "Potential Sunshine Hours"]} {} "i.sunhours: Potential Sunshine hours." {} -command {execute i.sunhours }}
 				{separator}
 				{command {[G_msg "Albedo"]} {} "i.albedo: Calculates Albedo from Modis, Aster, Landsat or AVHRR" {} -command {execute i.albedo }}
-				{separator}
 				{command {[G_msg "Vegetation Indices"]} {} "i.vi: Calculates 14 types of Vegetation Indices" {} -command {execute i.vi }}
+				{separator}
+				{command {[G_msg "Emissivity"]} {} "i.emissivity: Calculates emissivity from NDVI" {} -command {execute i.emissivity }}
 	}}
 	{command {[G_msg "Brovey sharpening"]} {} "i.fusion.brovey: Brovey transformation and pan sharpening" {} -command {execute i.fusion.brovey }}
 	{cascad {[G_msg "Classify image"]} {} "" $tmenu {			

gui/wxpython/gui_modules/menudata.py

@@ -1635,6 +1635,11 @@ class Data:
                                  _("Creates Latitude or Longitude map."),
                                  "self.OnMenuCmd",
                                  "i.latlong"),
+
+                                (_("Potential Sunshine Hours"),
+                                 _("Creates Potential Sunshine Hours map."),
+                                 "self.OnMenuCmd",
+                                 "i.sunhours"),
 	                        ("","","", ""),
 
                                 (_("Albedo"),
@@ -1646,6 +1651,12 @@ class Data:
                                  _("14 types of vegetation Indices."),
                                  "self.OnMenuCmd",
                                  "i.vi"),
+	                        ("","","", ""),
+
+                                (_("Emissivity"),
+                                 _("Emissivity from NDVI."),
+                                 "self.OnMenuCmd",
+                                 "i.emissivity"),
 				)
                          ),
                         ("","","", ""),

imagery/Makefile

@@ -7,6 +7,7 @@ SUBDIRS = \
 	i.atcorr \
 	i.cca \
 	i.cluster \
+	i.emissivity \
 	i.find \
 	i.gensig \
 	i.gensigset \
@@ -18,6 +19,7 @@ SUBDIRS = \
 	i.rectify \
 	i.rgb.his \
 	i.smap \
+	i.sunhours \
 	i.target \
 	i.pca \
 	i.vi 

imagery/i.emissivity/Makefile

@@ -0,0 +1,14 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.emissivity
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

imagery/i.emissivity/emissivity_generic.c

@@ -0,0 +1,21 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+
+    /* Emissivity Generic mode (Reads directly from NDVI)
+     * Estimation in the 8-14 micrometers range for sparse canopy
+     */ 
+double emissivity_generic(double ndvi) 
+{
+    double result;
+
+    if (ndvi < 0.16) 
+	result = 1.0;
+    else if (ndvi > 0.74) 
+	result = 0.9;
+    else 
+	result = 1.009 + 0.047 * log(ndvi);
+    return result;
+}
+
+

imagery/i.emissivity/i.emissivity.html

@@ -0,0 +1,25 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.emissivity</EM> calculates the emissivity in the longwave radiation spectrum, according to the semi-empirical equation related to NDVI by Caselles and Colles (1997), valid in the NDVI range of 0.16 to 0.74.
+
+Estimation in the 8-14 micrometers range for sparse canopy
+
+<H2>NOTES</H2>
+
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.eb.netrad.html">i.eb.netrad</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+Yann Chemin, GRASS Development Team<BR>
+
+
+<p>
+<i>Last changed: $Date: 2007/03/29 22:19:55 $</i>

imagery/i.emissivity/main.c

@@ -0,0 +1,114 @@
+/****************************************************************************
+ *
+ * MODULE:       i.emissivity
+ * AUTHOR(S):    Yann Chemin - yann.chemin@gmail.com
+ * PURPOSE:      Calculates the emissivity from NDVI (empirical)
+ *                as seen in Caselles and Colles (1997). 
+ *
+ * COPYRIGHT:    (C) 2002-2008 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/glocale.h>
+
+double emissivity_generic(double ndvi);
+
+int main(int argc, char *argv[]) 
+{
+    int nrows, ncols;
+    int row, col;
+    struct GModule *module;
+    struct Option *input, *output;
+    struct History history;	/*metadata */
+    
+    /************************************/ 
+    char *result1;    /*output raster name */
+    int infd, outfd;    /*File Descriptors */ 
+    char *ndvi;
+    char *emissivity;
+    void *inr;
+    DCELL * outr;
+
+    /************************************/ 
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    module->keywords = _("emissivity, land flux, energy balance");
+    module->description =
+	_("Emissivity from NDVI, generic method for spares land.");
+    
+    /* Define the different options */ 
+    input = G_define_standard_option(G_OPT_R_INPUT);
+    input->description = _("Name of the NDVI map [-]");
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+    output->description = _("Name of the output emissivity layer");
+    
+    /********************/ 
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    ndvi = input->answer;
+    result1 = output->answer;
+    
+    /***************************************************/ 
+    if ((infd = G_open_cell_old(ndvi, "")) < 0)
+	G_fatal_error(_("Cannot open cell file [%s]"), ndvi);
+    inr = G_allocate_d_raster_buf();
+    
+    /***************************************************/ 
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+    outr = G_allocate_d_raster_buf();
+    
+    /* Create New raster files */ 
+    if ((outfd = G_open_raster_new(result1, DCELL_TYPE)) < 0)
+	G_fatal_error(_("Could not open <%s>"), result1);
+    
+    /* Process pixels */ 
+    for (row = 0; row < nrows; row++)
+    {
+	DCELL d;
+	DCELL d_ndvi;
+	G_percent(row, nrows, 2);
+	
+        /* read input maps */ 
+        if (G_get_raster_row(infd,inr,row,DCELL_TYPE)< 0)
+	    G_fatal_error(_("Could not read from <%s>"), ndvi);
+	
+        /*process the data */ 
+        for (col = 0; col < ncols; col++)
+        {
+            d_ndvi = ((DCELL *) inr)[col];
+	    if (G_is_d_null_value(&d_ndvi)) 
+		G_set_d_null_value(&outr[col], 1);
+	    else {
+                /****************************/ 
+                /* calculate emissivity     */ 
+                d = emissivity_generic(d_ndvi);
+		outr[col] = d;
+	    }
+        }
+	if (G_put_raster_row(outfd, outr, DCELL_TYPE) < 0)
+	    G_fatal_error(_("Cannot write to output raster file"));
+    }
+    G_free(inr);
+    G_close_cell(infd);
+    G_free(outr);
+    G_close_cell(outfd);
+
+    G_short_history(result1, "raster", &history);
+    G_command_history(&history);
+    G_write_history(result1, &history);
+
+    exit(EXIT_SUCCESS);
+}
+
+

imagery/i.latlong/main.c

@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
     input1->description = _("Name of the input map");
 
     output1 = G_define_standard_option(G_OPT_R_OUTPUT);
-    output1->description = _("Name of the output longitude layer");
+    output1->description = _("Name of the output latitude or longitude layer");
     
     flag1 = G_define_flag();
     flag1->key = 'l';

imagery/i.sunhours/Makefile

@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.sunhours
+
+LIBES =  $(GPROJLIB) $(GISLIB)
+DEPENDENCIES = $(GPROJDEP) $(GISDEP)
+EXTRA_INC = $(PROJINC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_LARGEFILES),)
+	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
+endif
+
+default: cmd

imagery/i.sunhours/i.sunhours.html

@@ -0,0 +1,24 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.sunhours</EM> creates a sunshine hours map from any map, it considers a perfect clear day. This method follows Iqbal (1983) as found in the AHAS manual (Parodi, 2000).
+<H2>NOTES</H2>
+Iqbal, M., 1983. An Introduction to Solar Radiation. Iqbal, M., Editorial: Academic Press. Toronto, Canada.
+Parodi, G., 2000. AVHRR Hydrological Analysis System. Algorithms and Theory, Version 1.0. WRES - ITC, The Netherlands.
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.evapo.TSA.html">i.evapo.TSA</A><br>
+<A HREF="i.latitude.html">i.latitude</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, GRASS Development Team<BR>
+
+
+<p>
+<i>Last changed: $Date: 2007/02/19 10:19:59 $</i>

imagery/i.sunhours/main.c

@@ -0,0 +1,134 @@
+/****************************************************************************
+ *
+ * MODULE:       i.sunhours
+ * AUTHOR(S):    Yann Chemin - yann.chemin@gmail.com
+ * PURPOSE:      Calculates the sunshine hours (also called daytime period)
+ * 		 under a perfect clear sky condition.
+ * 		 Called generally "N" in meteorology. 
+ *
+ * COPYRIGHT:    (C) 2002-2008 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 <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+    
+#define PI 3.1415927
+
+int main(int argc, char *argv[]) 
+{
+    int nrows, ncols;
+    int row, col;
+    struct GModule *module;
+    struct Option *input1, *input2, *output1;
+    struct History history;	/*metadata */
+
+    /************************************/ 
+    char *name;			/*input raster name */
+    char *result1;		/*output raster name */
+    int infd_lat, infd_doy;    /*File Descriptors */ 
+    int outfd1;
+    char *lat, *doy;
+    void *inrast_lat, *inrast_doy;
+    DCELL * outrast1;
+
+    /************************************/ 
+    G_gisinit(argv[0]);
+    module = G_define_module();
+    module->keywords = _("sunshine, hours, daytime");
+    module->description = _("creates a sunshine hours map");
+    
+    /* Define the different options */ 
+    input1 = G_define_standard_option(G_OPT_R_INPUT);
+    input1->key = _("doy");
+    input1->description = _("Name of the doy input map");
+
+    input2 = G_define_standard_option(G_OPT_R_INPUT);
+    input2->key = _("lat");
+    input2->description = _("Name of the latitude input map");
+
+    output1 = G_define_standard_option(G_OPT_R_OUTPUT);
+    output1->description = _("Name of the output sunshine hours map");
+
+    /********************/ 
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    doy = input1->answer;
+    lat = input2->answer;
+    result1 = output1->answer;
+
+    /***************************************************/ 
+    if ((infd_doy = G_open_cell_old(doy, "")) < 0)
+	G_fatal_error(_("Cannot open cell file [%s]"), doy);
+    inrast_doy = G_allocate_d_raster_buf();
+
+    /***************************************************/ 
+    if ((infd_lat = G_open_cell_old(lat, "")) < 0)
+	G_fatal_error(_("Cannot open cell file [%s]"), lat);
+    inrast_lat = G_allocate_d_raster_buf();
+
+    /***************************************************/ 
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+
+    outrast1 = G_allocate_d_raster_buf();
+    if ((outfd1 = G_open_raster_new(result1, DCELL_TYPE)) < 0)
+	G_fatal_error(_("Could not open <%s>"), result1);
+
+    for (row = 0; row < nrows; row++)
+    {
+	DCELL d;
+	DCELL d_da;
+	DCELL d_delta;
+	DCELL d_Ws;
+	DCELL d_N;
+	DCELL d_lat;
+	DCELL d_doy;
+	G_percent(row, nrows, 2);
+
+	if (G_get_raster_row(infd_doy, inrast_doy, row, DCELL_TYPE) < 0)
+	    G_fatal_error(_("Could not read from <%s>"), doy);
+	if (G_get_raster_row(infd_lat, inrast_lat, row, DCELL_TYPE) < 0)
+	    G_fatal_error(_("Could not read from <%s>"), lat);
+
+	for (col = 0; col < ncols; col++)
+        {
+            d_doy = ((DCELL *) inrast_doy)[col];
+            d_lat = ((DCELL *) inrast_lat)[col];
+
+	    d_da = 2 * PI * (d_doy - 1) / 365.0;
+	    d_delta = 
+		0.006918 - 0.399912 * cos(d_da) + 0.070257 * sin(d_da) -
+		0.006758 * cos(2 * d_da) + 0.000907 * sin(2 * d_da) -
+		0.002697 * cos(3 * d_da) + 0.00148 * sin(3 * d_da);
+	    d_Ws = acos(-tan(d_lat * PI / 180) * tan(d_delta));
+	    d_N = (360.0 / (15.0 * PI)) * d_Ws;
+	    ((DCELL *) outrast1)[col] = d_N;
+        }
+	if (G_put_raster_row(outfd1, outrast1, DCELL_TYPE) < 0)
+	    G_fatal_error(_("Cannot write to output raster file"));
+    }
+    G_free(inrast_lat);
+    G_free(inrast_doy);
+    G_close_cell(infd_lat);
+    G_close_cell(infd_doy);
+    G_free(outrast1);
+    G_close_cell(outfd1);
+
+    G_short_history(result1, "raster", &history);
+    G_command_history(&history);
+    G_write_history(result1, &history);
+
+    exit(EXIT_SUCCESS);
+}
+
+