[GRASS-SVN] r30465 - in grass-addons: . i.topo.corr
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 4 09:54:24 EST 2008
Author: ejtizado
Date: 2008-03-04 09:54:24 -0500 (Tue, 04 Mar 2008)
New Revision: 30465
Added:
grass-addons/i.topo.corr/
grass-addons/i.topo.corr/Makefile
grass-addons/i.topo.corr/correction.c
grass-addons/i.topo.corr/description.html
grass-addons/i.topo.corr/illumination.c
grass-addons/i.topo.corr/local_proto.h
grass-addons/i.topo.corr/main.c
Log:
Topographic corrections
Added: grass-addons/i.topo.corr/Makefile
===================================================================
--- grass-addons/i.topo.corr/Makefile (rev 0)
+++ grass-addons/i.topo.corr/Makefile 2008-03-04 14:54:24 UTC (rev 30465)
@@ -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/i.topo.corr/correction.c
===================================================================
--- grass-addons/i.topo.corr/correction.c (rev 0)
+++ grass-addons/i.topo.corr/correction.c 2008-03-04 14:54:24 UTC (rev 30465)
@@ -0,0 +1,151 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include <unistd.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 = cos_z;
+ ckb = 0.;
+ kk = m;
+ G_message("Minnaert constant = %lf", kk);
+ break;
+ case C_CORRECT:
+ cka = cos_z + a/m;
+ ckb = a/m;
+ kk = 1.;
+ G_message("C-factor constant = %lf (a=%.4f; m=%.4f)", a/m, a, m);
+ break;
+ case PERCENT:
+ cka = 2.;
+ ckb = 1.;
+ kk = 1.;
+ break;
+ default: /* COSINE */
+ cka = 0.;
+ 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( cka/(cos_i + ckb), kk ) );
+ }
+ }
+ G_put_raster_row(out->fd, out->rast, out->type);
+ }
+}
Added: grass-addons/i.topo.corr/description.html
===================================================================
--- grass-addons/i.topo.corr/description.html (rev 0)
+++ grass-addons/i.topo.corr/description.html 2008-03-04 14:54:24 UTC (rev 30465)
@@ -0,0 +1,94 @@
+<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:
+<ul>
+ <li> cos_i = cos(s) · cos(z) + sin(s) · sin(z) · sin(a - o) </li>
+</ul>
+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>
+
+<p>
+<ul>
+ <li> ref_c = ref_o · cos_z / cos_i </li>
+</ul>
+</p>
+
+
+<H3>Method: minnaert</H3>
+
+<p>
+<ul>
+ <li>ref_c = ref_o · (cos_z / cos_i) ^ k</li>
+</ul>
+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>
+
+<p>
+<ul>
+ <li>ref_c = ref_o · (cos_z + c)/ (cos_i + c)</li>
+</ul>
+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
+<ul>
+ <li>ref_c = ref_o · 2 / (cos_i + 1)</li>
+</ul>
+</p>
+
+
+<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>
+
+<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>
+
+
+<H2>AUTHOR</H2>
+
+E. Jorge Tizado (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of León, Spain<BR>
+
+<p>
+<i>Last changed: $Date: 2008/02/10 00:00:00 $</i>
Added: grass-addons/i.topo.corr/illumination.c
===================================================================
--- grass-addons/i.topo.corr/illumination.c (rev 0)
+++ grass-addons/i.topo.corr/illumination.c 2008-03-04 14:54:24 UTC (rev 30465)
@@ -0,0 +1,270 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+void eval_c_cosi(Gfile * out, Gfile * dem, double zenith, double azimuth)
+{
+ struct Cell_head window;
+
+ int i, row, col, nrows, ncols;
+ CELL *cell[3], *temp;
+ CELL *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] = (CELL *) G_calloc (ncols + 1, sizeof(CELL));
+ G_set_c_null_value(cell[0], ncols);
+ cell[1] = (CELL *) G_calloc (ncols + 1, sizeof(CELL));
+ G_set_c_null_value(cell[1], ncols);
+ cell[2] = (CELL *) G_calloc (ncols + 1, sizeof(CELL));
+ G_set_c_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_c_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_c_null_value(c1) || G_is_c_null_value(c2) ||
+ G_is_c_null_value(c3) || G_is_c_null_value(c4) ||
+ G_is_c_null_value(c5) || G_is_c_null_value(c6) ||
+ G_is_c_null_value(c7) || G_is_c_null_value(c8) ||
+ G_is_c_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);
+}
+
+void eval_f_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);
+}
+
+void eval_d_cosi(Gfile * out, Gfile * dem, double zenith, double azimuth)
+{
+ struct Cell_head window;
+
+ int i, row, col, nrows, ncols;
+ DCELL *cell[3], *temp;
+ DCELL *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] = (DCELL *) G_calloc (ncols + 1, sizeof(DCELL));
+ G_set_d_null_value(cell[0], ncols);
+ cell[1] = (DCELL *) G_calloc (ncols + 1, sizeof(DCELL));
+ G_set_d_null_value(cell[1], ncols);
+ cell[2] = (DCELL *) G_calloc (ncols + 1, sizeof(DCELL));
+ G_set_d_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_d_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_d_null_value(c1) || G_is_d_null_value(c2) ||
+ G_is_d_null_value(c3) || G_is_d_null_value(c4) ||
+ G_is_d_null_value(c5) || G_is_d_null_value(c6) ||
+ G_is_d_null_value(c7) || G_is_d_null_value(c8) ||
+ G_is_d_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/i.topo.corr/local_proto.h
===================================================================
--- grass-addons/i.topo.corr/local_proto.h (rev 0)
+++ grass-addons/i.topo.corr/local_proto.h 2008-03-04 14:54:24 UTC (rev 30465)
@@ -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/i.topo.corr/main.c
===================================================================
--- grass-addons/i.topo.corr/main.c (rev 0)
+++ grass-addons/i.topo.corr/main.c 2008-03-04 14:54:24 UTC (rev 30465)
@@ -0,0 +1,251 @@
+/****************************************************************************
+ *
+ * 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 of reflectance");
+ module->keywords = _("Topographic correction, Cosine, Minnaert, C-Factor, Percent");
+
+ /* 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 to correct topographically");
+
+ 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 (flag -i) or prefix of output files");
+
+ 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,float");
+ zeni->description = _("Solar zenith in degrees");
+
+ azim = G_define_option();
+ azim->key = _("azimuth");
+ azim->type = TYPE_DOUBLE;
+ azim->required = NO;
+ azim->gisprompt = _("azimuth,float");
+ 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 = _("To 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 necessary to calculate ilumination terrain model.");
+
+ if (!ilum->answer && input->answer == NULL)
+ G_fatal_error("Reflectance files are necessary to make topographic correction.");
+
+ zenith = atof(zeni->answer);
+
+ /* Evaluate only cos_i 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 output file */
+ full_open_new( &out, output->answer, DCELL_TYPE );
+ out.rast = G_allocate_raster_buf(out.type);
+ /* Open and buffer of the elevation file */
+ if (dem.type == CELL_TYPE)
+ {
+ dem.rast = G_allocate_raster_buf(CELL_TYPE);
+ eval_c_cosi( &out, &dem, zenith, azimuth );
+ }
+ else if (dem.type == FCELL_TYPE)
+ {
+ dem.rast = G_allocate_raster_buf(FCELL_TYPE);
+ eval_f_cosi( &out, &dem, zenith, azimuth );
+ }
+ else if (dem.type == DCELL_TYPE)
+ {
+ dem.rast = G_allocate_raster_buf(DCELL_TYPE);
+ eval_d_cosi( &out, &dem, zenith, azimuth );
+ }
+ else
+ {
+ G_fatal_error(_("Elevation file of unknown type"));
+ }
+ /* 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-factor
+ 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);
+
+ if (metho->answer[1] == 'o') method = COSINE;
+ else if (metho->answer[1] == 'e') method = PERCENT;
+ else if (metho->answer[1] == 'i') method = MINNAERT;
+ else if (metho->answer[1] == '-') 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]);
+ 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);
+
+ char command[300];
+ sprintf(command, "r.colors map=%s color=grey", out.name);
+ system(command);
+ }
+ G_close_cell(dem.fd);
+ }
+
+ exit(EXIT_SUCCESS);
+}
+
More information about the grass-commit
mailing list