[GRASS-SVN] r35690 - in grass-addons/grass7/imagery: . i.topo.corr

svn_grass at osgeo.org svn_grass at osgeo.org
Sat Jan 31 11:37:37 EST 2009


Author: ejtizado
Date: 2009-01-31 11:37:37 -0500 (Sat, 31 Jan 2009)
New Revision: 35690

Added:
   grass-addons/grass7/imagery/i.topo.corr/
   grass-addons/grass7/imagery/i.topo.corr/Makefile
   grass-addons/grass7/imagery/i.topo.corr/correction.c
   grass-addons/grass7/imagery/i.topo.corr/i.topo.corr.html
   grass-addons/grass7/imagery/i.topo.corr/illumination.c
   grass-addons/grass7/imagery/i.topo.corr/local_proto.h
   grass-addons/grass7/imagery/i.topo.corr/main.c
Log:
Modified to compile with GRASS 7

Added: grass-addons/grass7/imagery/i.topo.corr/Makefile
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/Makefile	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/Makefile	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.topo.corr
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

Added: grass-addons/grass7/imagery/i.topo.corr/correction.c
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/correction.c	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/correction.c	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,146 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+void eval_tcor(int method, Gfile *out, Gfile *cosi, Gfile *band, double zenith)
+{
+	int i, row, col, nrows, ncols;
+	void *pref, *pcos;
+
+	double cos_z, cos_i, ref_i;
+	double n, sx, sxx, sy, sxy, tx, ty;
+	double a, m, cka, ckb, kk;
+
+
+	nrows = G_window_rows();
+	ncols = G_window_cols();
+
+	cos_z = cos(D2R * zenith);
+
+	/* Calculating regression */
+	if (method > NON_LAMBERTIAN)
+	{
+		n = sx = sxx = sy = sxy = 0.;
+		for (row = 0; row < nrows; row++)
+		{
+			G_percent(row, nrows, 2);
+
+			G_get_raster_row(band->fd, band->rast, row, DCELL_TYPE);
+			G_get_raster_row(cosi->fd, cosi->rast, row, cosi->type);
+
+			for (col = 0; col < ncols; col++)
+			{
+				switch( cosi->type )
+				{
+					case FCELL_TYPE:
+						pcos  = (void *)((FCELL *)cosi->rast + col);
+						cos_i = (double)((FCELL *)cosi->rast)[col];
+						break;
+					case DCELL_TYPE:
+						pcos  = (void *)((DCELL *)cosi->rast + col);
+						cos_i = (double)((DCELL *)cosi->rast)[col];
+						break;
+					default:
+						pcos  = NULL;
+						cos_i = 0.;
+						break;
+				}
+				pref = (void *)((DCELL *)band->rast + col);
+
+				if( !G_is_null_value(pref, DCELL_TYPE) &&
+					!G_is_null_value(pcos, cosi->type) )
+				{
+					ref_i = (double)*((DCELL *)pref);
+					switch( method )
+					{
+						case MINNAERT:
+							if (cos_i > 0. && cos_z > 0. && ref_i > 0.)
+							{
+								n++;
+								tx = log(cos_i/cos_z);
+								ty = log(ref_i);
+								sx  += tx; sxx += tx*tx;
+								sy  += ty; sxy += tx*ty;
+							}
+							break;
+						case C_CORRECT:
+							{
+								n++;
+								sx  += cos_i; sxx += cos_i*cos_i;
+								sy  += ref_i; sxy += cos_i*ref_i;
+							}
+							break;
+					}
+				}
+			}
+		}
+		m = (n == 0.) ? 1. : (n * sxy - sx * sy) / (n * sxx - sx * sx);
+		a = (n == 0.) ? 0. : (sy - m * sx) / n;
+	}
+	/* Calculating Constants */
+	switch( method )
+	{
+		case MINNAERT:
+			cka = ckb = 0.;
+			kk = m;
+			G_message("Minnaert constant = %lf", kk);
+			break;
+		case C_CORRECT:
+			cka = ckb = a/m;
+			kk = 1.;
+			G_message("C-factor constant = %lf (a = %.5f, m = %.5f)", cka, a, m);
+			break;
+		case PERCENT:
+			cka = 2. - cos_z; ckb = 1.;
+			kk = 1.;
+			break;
+		default: /* COSINE */
+			cka = ckb = 0.;
+			kk = 1.;
+	}
+	/* Topographic correction */
+	for (row = 0; row < nrows; row++)
+	{
+		G_percent(row, nrows, 2);
+
+		G_get_raster_row(band->fd, band->rast, row, band->type);
+		G_get_raster_row(cosi->fd, cosi->rast, row, cosi->type);
+
+		for (col = 0; col < ncols; col++)
+		{
+			switch( cosi->type )
+			{
+				case FCELL_TYPE:
+					pcos  = (void *)((FCELL *)cosi->rast + col);
+					cos_i = (double)*((FCELL *)pcos);
+					break;
+				case DCELL_TYPE:
+					pcos  = (void *)((DCELL *)cosi->rast + col);
+					cos_i = (double)*((DCELL *)pcos);
+					break;
+				default:
+					pcos  = NULL;
+					cos_i = 0.;
+					break;
+			}
+			pref = (void *)((DCELL *)band->rast + col);
+
+			if( pcos == NULL ||
+				G_is_null_value(pref, DCELL_TYPE) ||
+				G_is_null_value(pcos, cosi->type) )
+			{
+				G_set_null_value((DCELL *) out->rast + col, 1, DCELL_TYPE);
+			}
+			else
+			{
+				ref_i = (double)*((DCELL *)pref);
+				((DCELL *)out->rast)[col] = (DCELL)(ref_i * pow( (cos_z + cka)/(cos_i + ckb), kk ) );
+			}
+		}
+		G_put_raster_row(out->fd, out->rast, out->type);
+	}
+}

Added: grass-addons/grass7/imagery/i.topo.corr/i.topo.corr.html
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/i.topo.corr.html	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/i.topo.corr.html	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,89 @@
+<H2>DESCRIPTION</H2>
+
+<p><EM>i.topo.corr</EM> is to correct topographically reflectance files, e.g. obtained with <EM>i.landsat.toar</EM>, from a sun illumination terrain model. This illumination model represents the cosine of the incident angle, i.e. the  angle between the normal to the ground and the sun rays. This can be obtained with <em>r.sun</em> (parameter incidout) and then calculating their cosine with float precission.</p>
+
+<p>Using flag -i and given a elevation map as basemap (UTM), <em>i.topo.corr</em> permit to obtain a simple illumination model from the formulae:</p>
+<ul>
+	<li> cos_i = cos(s) · cos(z) + sin(s) · sin(z) · sin(a - o) </li>
+</ul>
+<p>where,
+<em>s</em> if the terrain slope angle, <em>z</em> is the solar zenith angle, <em>a</em> the solar azimuth angle, <em>o</em> the terrain aspect angle.</p>
+
+<p>For each band file, the corrected reflectance (ref_c) is calculate from the original reflectance (ref_o) using one of the three methods (one lambertian and two non-lambertian).</p>
+
+<H3>Method: cosine</H3>
+
+<ul>
+	<li> ref_c = ref_o · cos_z / cos_i </li>
+</ul>
+
+
+<H3>Method: minnaert</H3>
+
+<ul>
+	<li>ref_c = ref_o · (cos_z / cos_i) ^ k</li>
+</ul>
+<p>where,
+<em>k</em> is obtained by linear regression of ln(ref_o) = ln(ref_c) - k ln(cos_z/cos_i)</p>
+
+<H3>Method: c-factor</H3>
+
+
+<ul>
+	<li>ref_c = ref_o · (cos_z + c)/ (cos_i + c)</li>
+</ul>
+<p>where,
+<em>c</em> is a/m from ref_o = a + m cos_i.</p>
+
+<H3>Method: percent</H3>
+
+<p>I propose to test a new method. As cos_i varies from -1 to +1, then (cos_i + 1)/2 varied from 0 (surface in the side in opposition to the sun: infinite correction) to 1 (direct exhibition to the sun: no correction) and it can be calculated as</p>
+<ul>
+	<li>ref_c = ref_o · 2 / (cos_i + 1)</li>
+</ul>
+
+
+<H2>NOTES</H2>
+
+<ol>
+	<li>The illumination model (cos_i) with flag -i uses the actual region as limits and the resolution of the elevation map.</li>
+	<li>The topographic correction use the full reflectance file (set to null if is null cos_i) and their resolution.</li>
+	<li>The elevation map to calculate illumination model would have to be UTM.</li>
+</ol>
+
+<H2>EXAMPLES</H2>
+
+<p>First, make a illumination model from the SRTM map, and then make topographic correction of the bands toar.5, toar.4 and toar.3 with output as tcor.toar.5, tcor.toar.4, and tcor.toar.3 using c-factor (= c-correction) method.</p>
+
+<div class="code">
+<pre>
+i.topo.corr -i base=SRTM zenith=33.3631 azimuth=59.8897 out=cosi
+i.topo.corr base=cosi input=toar.5,toar.4,toar.3 out=tcor zenith=33.3631 method=c-factor</pre>
+</div>
+
+<H2>REFERENCES</H2>
+
+<ul>
+	<li>Law K.H. and Nichol J. Topographic Correction For Differential Illumination Effects On Ikonos Satellite Imagery.</li>
+	<li>Meyer P. et al. Radiometric Corrections Of Topographically Induced Effects On Landsat Tm Data In Alpine Terrain</li>
+	<li>Riaño et al. Assessment of Different Topographic Corrections in Landsat-TM Data for Mapping Vegetation Types. IEEE Transactions On Geoscience And Remote Sensing, Vol. 41, No. 5, May 2003</li>
+	<li>Twele A. and Erasmi S. Evaluating topographic correction algorithms for improved land cover discrimination in mountainous areas of Central Sulawesi. Göttinger Geographische Abhandlungen, vol. 113, 2005.</li>
+</ul>
+
+
+
+<H2>SEE ALSO</H2>
+
+<p><em>
+	<A HREF="i.landsat.toar">i.landsat.toar</A>,
+	<A HREF="r.mapcalc.html">r.mapcalc</A>,
+	<A HREF="r.sun">r.sun</A>
+</em></p>
+
+
+<H2>AUTHOR</H2>
+
+<p>E. Jorge Tizado  (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR></p>
+
+<p><i>Last changed: $Date: 2008/02/10 00:00:00 $</i></p>

Added: grass-addons/grass7/imagery/i.topo.corr/illumination.c
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/illumination.c	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/illumination.c	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,98 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+/*
+	Solo UTM mapas, sino mal calculado slope y aspect (ver r.slope.aspect para completarlo)
+	Mapa de alturas en metros.
+ */
+void eval_cosi(Gfile * out, Gfile * dem, double zenith, double azimuth)
+{
+	struct Cell_head window;
+
+	int i, row, col, nrows, ncols;
+	FCELL *cell[3], *temp;
+	FCELL *c1, *c2, *c3, *c4, *c5, *c6, *c7, *c8, *c9;
+	double H, V, dx, dy, key, north, east, south, west, center;
+	double cos_i, cos_z, sin_z, slope, aspect;
+
+	G_get_set_window (&window);
+
+	G_begin_distance_calculations();
+	north  = G_row_to_northing(0.5, &window);
+	center = G_row_to_northing(1.5, &window);
+	south  = G_row_to_northing(2.5, &window);
+	east   = G_col_to_easting(2.5, &window);
+	west   = G_col_to_easting(0.5, &window);
+	V = G_distance(east, north, east, south) * 4;
+	H = G_distance(east, center, west, center) * 4;
+
+	zenith  *= D2R;
+	azimuth *= D2R;
+
+	cos_z = cos(zenith);
+	sin_z = sin(zenith);
+
+	/* Making cos_i raster ... */
+	nrows = G_window_rows();
+	ncols = G_window_cols();
+
+	cell[0] = (FCELL *) G_calloc (ncols + 1, sizeof(FCELL));
+	G_set_f_null_value(cell[0], ncols);
+	cell[1] = (FCELL *) G_calloc (ncols + 1, sizeof(FCELL));
+	G_set_f_null_value(cell[1], ncols);
+	cell[2] = (FCELL *) G_calloc (ncols + 1, sizeof(FCELL));
+	G_set_f_null_value(cell[2], ncols);
+
+	/* First row is null */
+	G_set_null_value((DCELL *)out->rast, G_window_cols(), DCELL_TYPE);
+	G_put_raster_row(out->fd, out->rast, DCELL_TYPE);
+	/* Next rows ... */
+	for( row = 2; row < nrows; row++ )
+	{
+		G_percent(row, nrows, 2);
+		temp = cell[0];
+		cell[0] = cell[1];
+		cell[1] = cell[2];
+		cell[2] = temp;
+		G_get_f_raster_row_nomask (dem->fd, cell[2], row);
+
+		c1 = cell[0]; c2 = c1+1; c3 = c1+2;
+		c4 = cell[1]; c5 = c4+1; c6 = c4+2;
+		c7 = cell[2]; c8 = c7+1; c9 = c7+2;
+
+		for (col = 1; col < ncols-1; col++, c1++,c2++,c3++,c4++,c5++,c6++,c7++,c8++,c9++)
+		{
+			if (G_is_f_null_value(c1) || G_is_f_null_value(c2) ||
+				G_is_f_null_value(c3) || G_is_f_null_value(c4) ||
+				G_is_f_null_value(c5) || G_is_f_null_value(c6) ||
+				G_is_f_null_value(c7) || G_is_f_null_value(c8) ||
+				G_is_f_null_value(c9))
+			{
+				G_set_d_null_value((DCELL *) out->rast + col, 1);
+			}
+			else
+			{
+				dx = ((*c1 + *c4 + *c4 + *c7) - (*c3 + *c6 + *c6 + *c9)) / H;
+				dy = ((*c1 + *c2 + *c2 + *c3) - (*c7 + *c8 + *c8 + *c9)) / V;
+				key = dx*dx + dy*dy;
+				slope  = atan(sqrt(key));
+				aspect = atan2(dx,-dy);
+				if (aspect < 0.0) aspect += 2*PI;
+
+				cos_i = cos_z*cos(slope) + sin_z*sin(slope)*cos(azimuth - aspect);
+
+				((DCELL *)out->rast)[col] = (DCELL)cos_i;
+			}
+		}
+
+		G_put_raster_row(out->fd, out->rast, DCELL_TYPE);
+	}
+	/* Last row is null */
+	G_set_null_value((DCELL *)out->rast, G_window_cols(), DCELL_TYPE);
+	G_put_raster_row(out->fd, out->rast, DCELL_TYPE);
+}

Added: grass-addons/grass7/imagery/i.topo.corr/local_proto.h
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/local_proto.h	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/local_proto.h	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,28 @@
+#ifndef _LOCAL_PROTO_H
+#define _LOCAL_PROTO_H
+
+#define PI   3.1415926535897932384626433832795
+#define R2D 57.295779513082320877
+#define D2R  0.017453292519943295769
+
+typedef struct
+{
+	int     fd;
+	char    name[128];
+	char *	mapset;
+	void *  rast;
+	RASTER_MAP_TYPE type;
+} Gfile;
+
+#define LAMBERTIAN		 	 0
+#define COSINE			 	 1
+#define PERCENT				 2
+#define NON_LAMBERTIAN		10
+#define MINNAERT			11
+#define C_CORRECT			12
+
+
+void eval_cosi(Gfile *, Gfile *, double, double);
+void eval_tcor(int, Gfile *, Gfile *, Gfile *, double);
+
+#endif

Added: grass-addons/grass7/imagery/i.topo.corr/main.c
===================================================================
--- grass-addons/grass7/imagery/i.topo.corr/main.c	                        (rev 0)
+++ grass-addons/grass7/imagery/i.topo.corr/main.c	2009-01-31 16:37:37 UTC (rev 35690)
@@ -0,0 +1,235 @@
+/****************************************************************************
+ *
+ * MODULE:       i.topo.corr
+ *
+ * AUTHOR(S):    E. Jorge Tizado - ej.tizado at unileon.es
+ *
+ * PURPOSE:      Topographic corrections
+ *
+ * COPYRIGHT:    (C) 2002, 2005 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 <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+int full_open_old(Gfile * gf, char * fname)
+{
+	gf->fd = -1;
+
+	snprintf(gf->name, 127, "%s", fname);
+
+	gf->mapset = G_find_cell2(gf->name, "");
+	if( gf->mapset == NULL )
+		G_fatal_error("Cell file <%s> not found", gf->name);
+
+	gf->fd = G_open_cell_old(gf->name, gf->mapset);
+	if( gf->fd < 0 )
+		G_fatal_error(_("Cannot open file [%s] in mapset [%s]"), gf->name, gf->mapset);
+
+	gf->type = G_raster_map_type(gf->name, gf->mapset);
+
+	return gf->fd;
+}
+
+int full_open_new(Gfile * gf, char * fname, RASTER_MAP_TYPE ftype)
+{
+	gf->fd = -1;
+
+	snprintf(gf->name, 127, "%s", fname);
+	if( G_legal_filename(gf->name) < 0 )
+		G_fatal_error("[%s] is an illegal name", gf->name);
+
+	gf->type = ftype;
+	gf->mapset = G_mapset();
+	gf->fd  = G_open_raster_new(gf->name, gf->type);
+	if( gf->fd < 0 )
+		G_fatal_error(_("Cannot create file [%s] in mapset [%s]"), gf->name, gf->mapset);
+
+	return gf->fd;
+}
+
+int main(int argc, char *argv[])
+{
+	struct History history;
+	struct GModule *module;
+	struct Cell_head hd_band, hd_dem, window;
+
+	char bufname[128];
+	int nrows, ncols, row, col;
+
+	int i;
+	struct Option *base, *output, *input, *zeni, *azim, *metho;
+	struct Flag *ilum;
+
+	Gfile dem, out, band;
+	double zenith, azimuth;
+	int method = COSINE;
+
+	/* initialize GIS environment */
+	G_gisinit(argv[0]);
+
+	/* initialize module */
+	module = G_define_module();
+	module->description = _("Topografic correction");
+	module->keywords = _("Topographic correction, Cosine, Minnaert, C-Factor");
+
+	/* It defines the different parameters */
+
+	input = G_define_option();
+	input->key          = _("input");
+	input->type         = TYPE_STRING;
+	input->required     = NO;
+	input->multiple     = YES;
+	input->gisprompt    = _("input,cell,raster");
+	input->description  = _("List of reflectance band files");
+
+	output = G_define_option();
+	output->key         = _("output");
+	output->type        = TYPE_STRING;
+	output->required    = YES;
+	output->gisprompt   = _("output,cell,raster");
+	output->description = _("File name of output");
+
+	base = G_define_option();
+	base->key           = _("basemap");
+	base->type          = TYPE_STRING;
+	base->required      = YES;
+	base->gisprompt     = _("base,cell,raster");
+	base->description   = _("Base map for analysis (elevation or ilumination)");
+
+	zeni = G_define_option();
+	zeni->key          = _("zenith");
+	zeni->type         = TYPE_DOUBLE;
+	zeni->required     = YES;
+	zeni->gisprompt    = _("zenith");
+	zeni->description  = _("Solar zenith in degrees = 90 - [solar elevation])");
+
+	azim = G_define_option();
+	azim->key          = _("azimuth");
+	azim->type         = TYPE_DOUBLE;
+	azim->required     = NO;
+	azim->gisprompt    = _("azimuth");
+	azim->description  = _("Solar azimuth in degrees (only if flag -i)");
+
+	metho = G_define_option();
+	metho->key         = _("method");
+	metho->type        = TYPE_STRING;
+	metho->required    = NO;
+	metho->options     = "cosine,minnaert,c-factor,percent";
+	metho->gisprompt   = _("topographic correction method");
+	metho->description = _("Topographic correction method");
+	metho->answer      = "c-factor";
+
+	ilum = G_define_flag();
+	ilum->key          = 'i';
+	ilum->description  = _("Output sun ilumination terrain model");
+	ilum->answer       = 0;
+
+    if (G_parser(argc, argv))
+		exit(EXIT_FAILURE);
+
+	if (ilum->answer && azim->answer == NULL)
+		G_fatal_error("Solar azimuth is needed to calculate ilumination terrain model.");
+
+	if (!ilum->answer && input->answer == NULL)
+		G_fatal_error("Reflectance files are needed to topographic correction.");
+
+	zenith  = atof(zeni->answer);
+
+	/* Preparation files */
+	// Abrir el in aqui (común)
+
+	/* Evaluate only IL raster file */
+	// i.topo.corr -i out=cosi.on07 base=SRTM_v2 zenith=33.3631 azimuth=59.8897
+	if (ilum->answer)
+	{
+		azimuth = atof(azim->answer);
+		/* Warning: make buffers and output after set window */
+		full_open_old( &dem, base->answer );
+		/* Set window to DEM file */
+		G_get_window (&window);
+		G_get_cellhd (dem.name, dem.mapset, &hd_dem);
+		G_align_window (&window, &hd_dem);
+		G_set_window (&window);
+		/* Open and buffer of the elevation file */
+		if (G_raster_map_type(dem.name, dem.mapset) != FCELL_TYPE)
+			G_fatal_error(_("Elevation model must be FCELL"));
+		dem.rast  = G_allocate_raster_buf(FCELL_TYPE);
+		/* Open and buffer of the output file */
+		full_open_new( &out, output->answer, DCELL_TYPE );
+		out.rast = G_allocate_raster_buf(out.type);
+		/* -- -- -- */
+		eval_cosi( &out, &dem, zenith, azimuth );
+		/* Close files, buffers, and write history */
+		G_free(dem.rast);
+		G_close_cell(dem.fd);
+		G_free(out.rast);
+		G_close_cell(out.fd);
+		G_short_history(out.name, "raster", &history);
+		G_command_history(&history);
+		G_write_history(out.name, &history);
+	}
+	/* Evaluate topographic correction for all bands */
+	// i.topo.corr input=on07.toar.1 out=tcor base=cosi.on07 zenith=33.3631 method=c-correct
+	else
+	{
+		if (G_strcasecmp(metho->answer, "cosine") == 0)        method = COSINE;
+		else if (G_strcasecmp(metho->answer, "percent") == 0)  method = PERCENT;
+		else if (G_strcasecmp(metho->answer, "minnaert") == 0) method = MINNAERT;
+		else if (G_strcasecmp(metho->answer, "c-factor") == 0) method = C_CORRECT;
+		else G_fatal_error(_("Invalid method: %s"), metho->answer);
+
+		full_open_old( &dem, base->answer );
+		if (dem.type == CELL_TYPE)
+			G_fatal_error(_("Illumination model is CELL"));
+
+		for (i = 0; input->answers[i] != NULL; i++)
+		{
+			G_message("\nBand %s: ", input->answers[i]);
+			/* Abre fichero de bandas y el de salida */
+			full_open_old( &band, input->answers[i] );
+			if (band.type != DCELL_TYPE) {
+				G_warning(_("Reflectance of <%s> is not DCELL"), input->answers[i]);
+				G_close_cell(band.fd);
+				continue;
+			}
+			G_get_cellhd (band.name, band.mapset, &hd_band);
+			G_set_window (&hd_band); /* Antes de out_open y allocate para mismo tamaño */
+			/* ----- */
+			snprintf(bufname, 127, "%s.%s", output->answer, input->answers[i]);
+			////////////////////// QUE COJA SOLO EL NUMERO
+			full_open_new( &out, bufname, DCELL_TYPE );
+			out.rast  = G_allocate_raster_buf(out.type);
+			band.rast = G_allocate_raster_buf(band.type);
+			dem.rast  = G_allocate_raster_buf(dem.type);
+			/* ----- */
+			eval_tcor( method, &out, &dem, &band, zenith );
+			/* ----- */
+			G_free(dem.rast);
+
+			G_free(band.rast);
+			G_close_cell(band.fd);
+
+			G_free(out.rast);
+			G_close_cell(out.fd);
+
+			G_short_history(out.name, "raster", &history);
+			G_command_history(&history);
+			G_write_history(out.name, &history);
+		}
+
+		G_close_cell(dem.fd);
+	}
+
+    exit(EXIT_SUCCESS);
+}
+



More information about the grass-commit mailing list