[GRASS-SVN] r44969 - grass/trunk/imagery/i.eb.netrad
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 11 19:45:35 EST 2011
Author: ychemin
Date: 2011-01-11 16:45:35 -0800 (Tue, 11 Jan 2011)
New Revision: 44969
Added:
grass/trunk/imagery/i.eb.netrad/i.eb.netrad.html
Removed:
grass/trunk/imagery/i.eb.netrad/description.html
Modified:
grass/trunk/imagery/i.eb.netrad/Makefile
grass/trunk/imagery/i.eb.netrad/main.c
grass/trunk/imagery/i.eb.netrad/r_net.c
Log:
Ported i.eb.netrad to grass 7
Modified: grass/trunk/imagery/i.eb.netrad/Makefile
===================================================================
--- grass/trunk/imagery/i.eb.netrad/Makefile 2011-01-12 00:42:06 UTC (rev 44968)
+++ grass/trunk/imagery/i.eb.netrad/Makefile 2011-01-12 00:45:35 UTC (rev 44969)
@@ -2,8 +2,8 @@
PGM = i.eb.netrad
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Deleted: grass/trunk/imagery/i.eb.netrad/description.html
===================================================================
--- grass/trunk/imagery/i.eb.netrad/description.html 2011-01-12 00:42:06 UTC (rev 44968)
+++ grass/trunk/imagery/i.eb.netrad/description.html 2011-01-12 00:45:35 UTC (rev 44969)
@@ -1,25 +0,0 @@
-<H2>DESCRIPTION</H2>
-
-<EM>i.eb.netrad</EM> calculates the net radiation at the time of satellite overpass, the way it is in the SEBAL model of bastiaanssen (1995).
-
-It takes input of Albedo, NDVI, Surface Skin temperature, time of satellite overpass, surface emissivity, difference of temperature from surface skin and about 2 m height (dT), instantaneous satellite overpass single-way atmospheric transmissivity (tsw), Day of Year (DOY), and sun zenith angle.
-
-<H2>NOTES</H2>
-In the old methods, dT was taken as flat images (dT=5.0), if you don't have a dT map from ground data, you would want to try something in this line, this is to calculate atmospherical energy balance. In the same way, a standard tsw is used in those equations. Refer to r_net.c for that and for other non-used equations, but stored in there for further research convenience.
-<H2>TODO</H2>
-
-<H2>SEE ALSO</H2>
-
-<em>
-<A HREF="i.eb.g0.html">i.eb.g0</A><br>
-<A HREF="i.albedo.html">i.albedo</A><br>
-</em>
-
-
-<H2>AUTHORS</H2>
-
-Yann Chemin, International Rice Research Institute, The Philippines<BR>
-
-
-<p>
-<i>Last changed: $Date$</i>
Copied: grass/trunk/imagery/i.eb.netrad/i.eb.netrad.html (from rev 44968, grass/trunk/imagery/i.eb.netrad/description.html)
===================================================================
--- grass/trunk/imagery/i.eb.netrad/i.eb.netrad.html (rev 0)
+++ grass/trunk/imagery/i.eb.netrad/i.eb.netrad.html 2011-01-12 00:45:35 UTC (rev 44969)
@@ -0,0 +1,25 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.eb.netrad</EM> calculates the net radiation at the time of satellite overpass, the way it is in the SEBAL model of bastiaanssen (1995).
+
+It takes input of Albedo, NDVI, Surface Skin temperature, time of satellite overpass, surface emissivity, difference of temperature from surface skin and about 2 m height (dT), instantaneous satellite overpass single-way atmospheric transmissivity (tsw), Day of Year (DOY), and sun zenith angle.
+
+<H2>NOTES</H2>
+In the old methods, dT was taken as flat images (dT=5.0), if you don't have a dT map from ground data, you would want to try something in this line, this is to calculate atmospherical energy balance. In the same way, a standard tsw is used in those equations. Refer to r_net.c for that and for other non-used equations, but stored in there for further research convenience.
+<H2>TODO</H2>
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.eb.g0.html">i.eb.g0</A><br>
+<A HREF="i.albedo.html">i.albedo</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, International Rice Research Institute, The Philippines<BR>
+
+
+<p>
+<i>Last changed: $Date$</i>
Modified: grass/trunk/imagery/i.eb.netrad/main.c
===================================================================
--- grass/trunk/imagery/i.eb.netrad/main.c 2011-01-12 00:42:06 UTC (rev 44968)
+++ grass/trunk/imagery/i.eb.netrad/main.c 2011-01-12 00:45:35 UTC (rev 44969)
@@ -7,7 +7,7 @@
* as seen in Bastiaanssen (1995) using time of
* satellite overpass.
*
- * COPYRIGHT: (C) 2006-2008 by the GRASS Development Team
+ * 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
@@ -20,6 +20,7 @@
#include <stdlib.h>
#include <string.h>
#include <grass/gis.h>
+#include <grass/raster.h>
#include <grass/glocale.h>
double r_net(double bbalb, double ndvi, double tempk, double dtair,
double e0, double tsw, double doy, double utc,
@@ -27,123 +28,96 @@
int main(int argc, char *argv[])
{
struct Cell_head cellhd; /*region+header info */
-
char *mapset; /*mapset name */
-
int nrows, ncols;
-
int row, col;
-
struct GModule *module;
-
struct Option *input1, *input2, *input3, *input4, *input5;
-
struct Option *input6, *input7, *input8, *input9, *output1;
-
struct Flag *flag1;
-
struct History history; /*metadata */
-
struct Colors colors; /*Color rules */
-
-
-
- /************************************/
- /* FMEO Declarations**************** */
char *name; /*input raster name */
-
char *result; /*output raster name */
-
-
- /*File Descriptors */
int infd_albedo, infd_ndvi, infd_tempk, infd_time, infd_dtair;
-
int infd_emissivity, infd_tsw, infd_doy, infd_sunzangle;
-
int outfd;
-
char *albedo, *ndvi, *tempk, *time, *dtair, *emissivity;
-
char *tsw, *doy, *sunzangle;
-
int i = 0, j = 0;
-
void *inrast_albedo, *inrast_ndvi, *inrast_tempk, *inrast_rnet;
-
void *inrast_time, *inrast_dtair, *inrast_emissivity, *inrast_tsw;
-
void *inrast_doy, *inrast_sunzangle;
-
DCELL * outrast;
- RASTER_MAP_TYPE data_type_output = DCELL_TYPE;
- RASTER_MAP_TYPE data_type_albedo;
- RASTER_MAP_TYPE data_type_ndvi;
- RASTER_MAP_TYPE data_type_tempk;
- RASTER_MAP_TYPE data_type_time;
- RASTER_MAP_TYPE data_type_dtair;
- RASTER_MAP_TYPE data_type_emissivity;
- RASTER_MAP_TYPE data_type_tsw;
- RASTER_MAP_TYPE data_type_doy;
- RASTER_MAP_TYPE data_type_sunzangle;
-
-
- /************************************/
- G_gisinit(argv[0]);
+ CELL val1,val2; /*For color range*/
+ /************************************/
+ G_gisinit(argv[0]);
module = G_define_module();
- module->keywords = _("net radiation, energy balance, SEBAL");
+ G_add_keyword(_("net radiation, energy balance, SEBAL"));
+ G_add_keyword(_("energy balance"));
+ G_add_keyword(_("SEBAL"));
module->description =
_("net radiation approximation (Bastiaanssen, 1995)");
/* Define the different options */
- input1 = G_define_standard_option(G_OPT_R_INPUT);
+ input1 = G_define_standard_option(G_OPT_R_INPUT);
input1->key = _("albedo");
input1->description = _("Name of the Albedo map [0.0;1.0]");
input1->answer = _("albedo");
+
input2 = G_define_standard_option(G_OPT_R_INPUT);
input2->key = _("ndvi");
input2->description = _("Name of the ndvi map [-1.0;+1.0]");
input2->answer = _("ndvi");
+
input3 = G_define_standard_option(G_OPT_R_INPUT);
input3->key = _("tempk");
input3->description =
_("Name of the Surface temperature map [degree Kelvin]");
input3->answer = _("tempk");
+
input4 = G_define_standard_option(G_OPT_R_INPUT);
input4->key = _("time");
input4->description =
_("Name of the map of local UTC time of satellite overpass [hh.hhh]");
input4->answer = _("time");
+
input5 = G_define_standard_option(G_OPT_R_INPUT);
input5->key = _("dtair");
input5->description =
_("Name of the difference of temperature from surface skin to about 2 m height [K]");
input5->answer = _("dtair");
+
input6 = G_define_standard_option(G_OPT_R_INPUT);
input6->key = _("emissivity");
input6->description = _("Name of the emissivity map [-]");
input6->answer = _("emissivity");
+
input7 = G_define_standard_option(G_OPT_R_INPUT);
input7->key = _("tsw");
input7->description =
_("Name of the single-way atmospheric transmissivitymap [-]");
input7->answer = _("tsw");
+
input8 = G_define_standard_option(G_OPT_R_INPUT);
input8->key = _("doy");
input8->description = _("Name of the Day Of Year (DOY) map [-]");
input8->answer = _("doy");
+
input9 = G_define_standard_option(G_OPT_R_INPUT);
input9->key = _("sunzangle");
input9->description = _("Name of the sun zenith angle map [degrees]");
input9->answer = _("sunzangle");
+
output1 = G_define_standard_option(G_OPT_R_INPUT);
output1->key = _("rnet");
output1->description = _("Name of the output rnet layer");
output1->answer = _("rnet");
-
- /********************/
- if (G_parser(argc, argv))
+ /********************/
+ if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+
albedo = input1->answer;
ndvi = input2->answer;
tempk = input3->answer;
@@ -155,137 +129,43 @@
sunzangle = input9->answer;
result = output1->answer;
-
- /***************************************************/
- mapset = G_find_cell2(albedo, "");
- if (mapset == NULL) {
- G_fatal_error(_("cell file [%s] not found"), albedo);
- }
- data_type_albedo = G_raster_map_type(albedo, mapset);
- if ((infd_albedo = G_open_cell_old(albedo, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), albedo);
- if (G_get_cellhd(albedo, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s])"), albedo);
- inrast_albedo = G_allocate_raster_buf(data_type_albedo);
+ /* Open access to image files and allocate row access memory */
+ infd_albedo = Rast_open_old(albedo, "");
+ inrast_albedo = Rast_allocate_d_buf();
+ infd_ndvi = Rast_open_old(ndvi, "");
+ inrast_ndvi = Rast_allocate_d_buf();
- /***************************************************/
- mapset = G_find_cell2(ndvi, "");
- if (mapset == NULL) {
- G_fatal_error(_("cell file [%s] not found"), ndvi);
- }
- data_type_ndvi = G_raster_map_type(ndvi, mapset);
- if ((infd_ndvi = G_open_cell_old(ndvi, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), ndvi);
- if (G_get_cellhd(ndvi, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), ndvi);
- inrast_ndvi = G_allocate_raster_buf(data_type_ndvi);
-
+ infd_tempk = Rast_open_old(tempk, "");
+ inrast_tempk = Rast_allocate_d_buf();
- /***************************************************/
- mapset = G_find_cell2(tempk, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), tempk);
- }
- data_type_tempk = G_raster_map_type(tempk, mapset);
- if ((infd_tempk = G_open_cell_old(tempk, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), tempk);
- if (G_get_cellhd(tempk, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), tempk);
- inrast_tempk = G_allocate_raster_buf(data_type_tempk);
-
+ infd_dtair = Rast_open_old(dtair, "");
+ inrast_dtair = Rast_allocate_d_buf();
- /***************************************************/
- mapset = G_find_cell2(dtair, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), dtair);
- }
- data_type_dtair = G_raster_map_type(dtair, mapset);
- if ((infd_dtair = G_open_cell_old(dtair, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), dtair);
- if (G_get_cellhd(dtair, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), dtair);
- inrast_dtair = G_allocate_raster_buf(data_type_dtair);
-
+ infd_time = Rast_open_old(time, "");
+ inrast_time = Rast_allocate_d_buf();
- /***************************************************/
- mapset = G_find_cell2(time, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), time);
- }
- data_type_time = G_raster_map_type(time, mapset);
- if ((infd_time = G_open_cell_old(time, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), time);
- if (G_get_cellhd(time, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), time);
- inrast_time = G_allocate_raster_buf(data_type_time);
+ infd_emissivity = Rast_open_old(emissivity, "");
+ inrast_emissivity = Rast_allocate_d_buf();
-
- /***************************************************/
- mapset = G_find_cell2(emissivity, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), emissivity);
- }
- data_type_emissivity = G_raster_map_type(emissivity, mapset);
- if ((infd_emissivity = G_open_cell_old(emissivity, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), emissivity);
- if (G_get_cellhd(emissivity, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), emissivity);
- inrast_emissivity = G_allocate_raster_buf(data_type_emissivity);
+ infd_tsw = Rast_open_old(tsw, "");
+ inrast_tsw = Rast_allocate_d_buf();
-
- /***************************************************/
- mapset = G_find_cell2(tsw, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), tsw);
- }
- data_type_tsw = G_raster_map_type(tsw, mapset);
- if ((infd_tsw = G_open_cell_old(tsw, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), tsw);
- if (G_get_cellhd(tsw, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), tsw);
- inrast_tsw = G_allocate_raster_buf(data_type_tsw);
+ infd_doy = Rast_open_old(doy, "");
+ inrast_doy = Rast_allocate_d_buf();
-
- /***************************************************/
- mapset = G_find_cell2(doy, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), doy);
- }
- data_type_doy = G_raster_map_type(doy, mapset);
- if ((infd_doy = G_open_cell_old(doy, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), doy);
- if (G_get_cellhd(doy, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), doy);
- inrast_doy = G_allocate_raster_buf(data_type_doy);
+ infd_sunzangle = Rast_open_old(sunzangle, "");
+ inrast_sunzangle = Rast_allocate_d_buf();
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
- /***************************************************/
- mapset = G_find_cell2(sunzangle, "");
- if (mapset == NULL) {
- G_fatal_error(_("Cell file [%s] not found"), sunzangle);
- }
- data_type_sunzangle = G_raster_map_type(sunzangle, mapset);
- if ((infd_sunzangle = G_open_cell_old(sunzangle, mapset)) < 0)
- G_fatal_error(_("Cannot open cell file [%s]"), sunzangle);
- if (G_get_cellhd(sunzangle, mapset, &cellhd) < 0)
- G_fatal_error(_("Cannot read file header of [%s]"), sunzangle);
- inrast_sunzangle = G_allocate_raster_buf(data_type_sunzangle);
+ outfd = Rast_open_new(result, DCELL_TYPE);
+ outrast = Rast_allocate_d_buf();
-
- /***************************************************/
- G_debug(3, "number of rows %d", cellhd.rows);
- nrows = G_window_rows();
- ncols = G_window_cols();
- outrast = G_allocate_raster_buf(data_type_output);
-
- /* Create New raster files */
- if ((outfd = G_open_raster_new(result, data_type_output)) < 0)
- G_fatal_error(_("Could not open <%s>"), result);
-
- /* Process pixels */
- for (row = 0; row < nrows; row++)
- {
+ /* Process pixels */
+ for (row = 0; row < nrows; row++)
+ {
DCELL d;
DCELL d_albedo;
DCELL d_ndvi;
@@ -296,160 +176,54 @@
DCELL d_tsw;
DCELL d_doy;
DCELL d_sunzangle;
+
+ /* Display row process percentage */
G_percent(row, nrows, 2);
+
+ /* Load rows for each input image */
+ Rast_get_d_row(infd_albedo, inrast_albedo, row);
+ Rast_get_d_row(infd_ndvi, inrast_ndvi, row);
+ Rast_get_d_row(infd_tempk, inrast_tempk, row);
+ Rast_get_d_row(infd_dtair, inrast_dtair, row);
+ Rast_get_d_row(infd_time, inrast_time, row);
+ Rast_get_d_row(infd_emissivity, inrast_emissivity, row);
+ Rast_get_d_row(infd_tsw, inrast_tsw, row);
+ Rast_get_d_row(infd_doy, inrast_doy, row);
+ Rast_get_d_row(infd_sunzangle, inrast_sunzangle, row);
- /* read input maps */
- if (G_get_raster_row
- (infd_albedo, inrast_albedo, row, data_type_albedo) < 0)
- G_fatal_error(_("Could not read from <%s>"), albedo);
- if (G_get_raster_row(infd_ndvi, inrast_ndvi, row, data_type_ndvi) <
- 0)
- G_fatal_error(_("Could not read from <%s>"), ndvi);
- if (G_get_raster_row(infd_tempk, inrast_tempk, row, data_type_tempk)
- < 0)
- G_fatal_error(_("Could not read from <%s>"), tempk);
- if (G_get_raster_row(infd_dtair, inrast_dtair, row, data_type_dtair)
- < 0)
- G_fatal_error(_("Could not read from <%s>"), dtair);
- if (G_get_raster_row(infd_time, inrast_time, row, data_type_time) <
- 0)
- G_fatal_error(_("Could not read from <%s>"), time);
- if (G_get_raster_row
- (infd_emissivity, inrast_emissivity, row,
- data_type_emissivity) < 0)
- G_fatal_error(_("Could not read from <%s>"), emissivity);
- if (G_get_raster_row(infd_tsw, inrast_tsw, row, data_type_tsw) < 0)
- G_fatal_error(_("Could not read from <%s>"), tsw);
- if (G_get_raster_row(infd_doy, inrast_doy, row, data_type_doy) < 0)
- G_fatal_error(_("Could not read from <%s>"), doy);
- if (G_get_raster_row
- (infd_sunzangle, inrast_sunzangle, row, data_type_sunzangle) < 0)
- G_fatal_error(_("Could not read from <%s>"), sunzangle);
-
- /*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;
+ /*process the data */
+ for (col = 0; col < ncols; col++)
+ {
+ d_albedo = (double)((DCELL *) inrast_albedo)[col];
+ d_ndvi = (double)((DCELL *) inrast_ndvi)[col];
+ d_tempk = (double)((DCELL *) inrast_tempk)[col];
+ d_dtair = (double)((DCELL *) inrast_dtair)[col];
+ d_time = (double)((DCELL *) inrast_time)[col];
+ d_emissivity = (double)((DCELL *) inrast_emissivity)[col];
+ d_tsw = (double)((DCELL *) inrast_tsw)[col];
+ d_doy = (double)((DCELL *) inrast_doy)[col];
+ d_sunzangle = (double)((DCELL *) inrast_sunzangle)[col];
+ /* process NULL Values */
+ if (Rast_is_d_null_value(&d_albedo) ||
+ Rast_is_d_null_value(&d_ndvi) ||
+ Rast_is_d_null_value(&d_tempk) ||
+ Rast_is_d_null_value(&d_dtair) ||
+ Rast_is_d_null_value(&d_time) ||
+ Rast_is_d_null_value(&d_emissivity) ||
+ Rast_is_d_null_value(&d_tsw) ||
+ Rast_is_d_null_value(&d_doy) ||
+ Rast_is_d_null_value(&d_sunzangle)) {
+ Rast_set_d_null_value(&outrast[col], 1);
}
- 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 = ((DCELL *) inrast_ndvi)[col];
- break;
- }
- switch (data_type_tempk) {
- case CELL_TYPE:
- d_tempk = (double)((CELL *) inrast_tempk)[col];
- break;
- case FCELL_TYPE:
- d_tempk = (double)((FCELL *) inrast_tempk)[col];
- break;
- case DCELL_TYPE:
- d_tempk = (double)((DCELL *) inrast_tempk)[col];
- break;
- }
- switch (data_type_dtair) {
- case CELL_TYPE:
- d_dtair = (double)((CELL *) inrast_dtair)[col];
- break;
- case FCELL_TYPE:
- d_dtair = (double)((FCELL *) inrast_dtair)[col];
- break;
- case DCELL_TYPE:
- d_dtair = (double)((DCELL *) inrast_dtair)[col];
- break;
- }
- switch (data_type_time) {
- case CELL_TYPE:
- d_time = (double)((CELL *) inrast_time)[col];
- break;
- case FCELL_TYPE:
- d_time = (double)((FCELL *) inrast_time)[col];
- break;
- case DCELL_TYPE:
- d_time = (double)((DCELL *) inrast_time)[col];
- break;
- }
- switch (data_type_emissivity) {
- case CELL_TYPE:
- d_emissivity = (double)((CELL *) inrast_emissivity)[col];
- break;
- case FCELL_TYPE:
- d_emissivity = (double)((FCELL *) inrast_emissivity)[col];
- break;
- case DCELL_TYPE:
- d_emissivity = (double)((DCELL *) inrast_emissivity)[col];
- break;
- }
- switch (data_type_tsw) {
- case CELL_TYPE:
- d_tsw = (double)((CELL *) inrast_tsw)[col];
- break;
- case FCELL_TYPE:
- d_tsw = (double)((FCELL *) inrast_tsw)[col];
- break;
- case DCELL_TYPE:
- d_tsw = (double)((DCELL *) inrast_tsw)[col];
- break;
- }
- switch (data_type_doy) {
- case CELL_TYPE:
- d_doy = (double)((CELL *) inrast_doy)[col];
- break;
- case FCELL_TYPE:
- d_doy = (double)((FCELL *) inrast_doy)[col];
- break;
- case DCELL_TYPE:
- d_doy = (double)((DCELL *) inrast_doy)[col];
- break;
- }
- switch (data_type_sunzangle) {
- case CELL_TYPE:
- d_sunzangle = (double)((CELL *) inrast_sunzangle)[col];
- break;
- case FCELL_TYPE:
- d_sunzangle = (double)((FCELL *) inrast_sunzangle)[col];
- break;
- case DCELL_TYPE:
- d_sunzangle = (double)((DCELL *) inrast_sunzangle)[col];
- break;
- }
- if (G_is_d_null_value(&d_albedo) || G_is_d_null_value(&d_ndvi)
- || G_is_d_null_value(&d_tempk) ||
- G_is_d_null_value(&d_dtair) || G_is_d_null_value(&d_time)
- || G_is_d_null_value(&d_emissivity) ||
- G_is_d_null_value(&d_tsw) || G_is_d_null_value(&d_doy) ||
- G_is_d_null_value(&d_sunzangle)) {
- G_set_d_null_value(&outrast[col], 1);
- }
else {
-
-
- /************************************/
- /* calculate the net radiation */
- d =
- r_net(d_albedo, d_ndvi, d_tempk, d_dtair, d_emissivity,
- d_tsw, d_doy, d_time, d_sunzangle);
- outrast[col] = d;
+ /************************************/
+ /* calculate the net radiation */
+ d = r_net(d_albedo, d_ndvi, d_tempk, d_dtair, d_emissivity, d_tsw, d_doy, d_time, d_sunzangle);
+ outrast[col] = d;
}
- }
- if (G_put_raster_row(outfd, outrast, data_type_output) < 0)
- G_fatal_error(_("Cannot write to output raster file"));
}
+ Rast_put_d_row(outfd, outrast);
+ }
G_free(inrast_albedo);
G_free(inrast_ndvi);
G_free(inrast_tempk);
@@ -459,26 +233,28 @@
G_free(inrast_tsw);
G_free(inrast_doy);
G_free(inrast_sunzangle);
- G_close_cell(infd_albedo);
- G_close_cell(infd_ndvi);
- G_close_cell(infd_tempk);
- G_close_cell(infd_dtair);
- G_close_cell(infd_time);
- G_close_cell(infd_emissivity);
- G_close_cell(infd_tsw);
- G_close_cell(infd_doy);
- G_close_cell(infd_sunzangle);
+ Rast_close(infd_albedo);
+ Rast_close(infd_ndvi);
+ Rast_close(infd_tempk);
+ Rast_close(infd_dtair);
+ Rast_close(infd_time);
+ Rast_close(infd_emissivity);
+ Rast_close(infd_tsw);
+ Rast_close(infd_doy);
+ Rast_close(infd_sunzangle);
G_free(outrast);
- G_close_cell(outfd);
+ Rast_close(outfd);
- /* Colors in grey shade */
- G_init_colors(&colors);
- G_add_color_rule(0, 0, 0, 0, 900, 255, 255, 255, &colors);
+ /* Colors in grey shade */
+ Rast_init_colors(&colors);
+ val1=0;
+ val2=900;
+ Rast_add_c_color_rule(&val1, 0, 0, 0, &val2, 255, 255, 255, &colors);
- /* Metadata */
- G_short_history(result, "raster", &history);
- G_command_history(&history);
- G_write_history(result, &history);
+ /* Metadata */
+ Rast_short_history(result, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(result, &history);
exit(EXIT_SUCCESS);
}
Modified: grass/trunk/imagery/i.eb.netrad/r_net.c
===================================================================
--- grass/trunk/imagery/i.eb.netrad/r_net.c 2011-01-12 00:42:06 UTC (rev 44968)
+++ grass/trunk/imagery/i.eb.netrad/r_net.c 2011-01-12 00:45:35 UTC (rev 44969)
@@ -8,39 +8,37 @@
double sunzangle)
{
- /* Tsw = atmospheric transmissivity single-way (~0.7 -) */
- /* DOY = Day of Year */
- /* utc = UTC time of sat overpass */
- /* sunzangle = sun zenith angle at sat. overpass */
- /* tair = air temperature (approximative, or met.station) */
+ /* Tsw = atmospheric transmissivity single-way (~0.7 -) */
+ /* DOY = Day of Year */
+ /* utc = UTC time of sat overpass */
+ /* sunzangle = sun zenith angle at sat. overpass */
+ /* tair = air temperature (approximative, or met.station) */
double Kin = 0.0, Lin = 0.0, Lout = 0.0, Lcorr = 0.0, result = 0.0;
-
double temp = 0.0, ds = 0.0, e_atm = 0.0, delta = 0.0;
-
double tsw_for_e_atm = 0.7; /*Special tsw, consider it a coefficient */
-
- /* Atmospheric emissivity (Bastiaanssen, 1995) */
- e_atm = 1.08 * pow(-log(tsw), 0.265);
+ /* Atmospheric emissivity (Bastiaanssen, 1995) */
+ e_atm = 1.08 * pow(-log(tsw), 0.265);
- /* Atmospheric emissivity (Pawan, 2004) */
- /* e_atm = 0.85 * pow(-log(tsw),0.09); */
+ /* Atmospheric emissivity (Pawan, 2004) */
+ /* e_atm = 0.85 * pow(-log(tsw),0.09); */
- /* ds = 1.0 + 0.01672 * sin(2*PI*(doy-93.5)/365); */
- ds = 1.0 / pow((1 + 0.033 * cos(2 * PI * doy / 365)), 2);
+ /* ds = 1.0 + 0.01672 * sin(2*PI*(doy-93.5)/365); */
+ ds = 1.0 / pow((1 + 0.033 * cos(2 * PI * doy / 365)), 2);
delta = 0.4093 * sin((2 * PI * doy / 365) - 1.39);
- /* Kin is the shortwave incoming radiation */
- Kin = 1358.0 * (cos(sunzangle * PI / 180) * tsw / (ds * ds));
+ /* Kin is the shortwave incoming radiation */
+ Kin = 1358.0 * (cos(sunzangle * PI / 180) * tsw / (ds * ds));
- /* Lin is incoming longwave radiation */
- Lin = (e_atm) * 5.67 * pow(10, -8) * pow((tempk - dtair), 4);
+ /* Lin is incoming longwave radiation */
+ Lin = (e_atm) * 5.67 * pow(10, -8) * pow((tempk - dtair), 4);
- /* Lout is surface grey body emission in Longwave spectrum */
- Lout = e0 * 5.67 * pow(10, -8) * pow(tempk, 4);
+ /* Lout is surface grey body emission in Longwave spectrum */
+ Lout = e0 * 5.67 * pow(10, -8) * pow(tempk, 4);
- /* Lcorr is outgoing longwave radiation "reflected" by the emissivity */
- Lcorr = (1.0 - e0) * Lin;
+ /* Lcorr is outgoing longwave radiation "reflected" by the emissivity */
+ Lcorr = (1.0 - e0) * Lin;
+
result = (1.0 - bbalb) * Kin + Lin - Lout - Lcorr;
return result;
}
More information about the grass-commit
mailing list