move i.topo.corr to trunk (sync to 6.x)

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

imagery/i.topo.corr/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.topo.corr
+
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

imagery/i.topo.corr/correction.c

@@ -0,0 +1,141 @@
+/* File: correction.c
+ *
+ *  AUTHOR:    E. Jorge Tizado, Spain 2010
+ *
+ *  COPYRIGHT: (c) 2007-10 E. Jorge Tizado
+ *             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 <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
+	       double zenith)
+{
+    int 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 = Rast_window_rows();
+    ncols = Rast_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);
+
+	    Rast_get_row(band->fd, band->rast, row, band->type);
+	    Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
+	    
+	    pref = band->rast;
+	    pcos = cosi->rast;
+
+	    for (col = 0; col < ncols; col++) {
+		
+		cos_i = Rast_get_d_value(pcos, cosi->type);
+
+		if (!Rast_is_null_value(pref, band->type) &&
+		    !Rast_is_null_value(pcos, cosi->type)) {
+		    ref_i = Rast_get_d_value(pref, band->type);
+		    switch (method) {
+		    case MINNAERT:
+			if (cos_i > 0. && cos_z > 0. && ref_i > 0.) {
+			    n++;
+                            /* tx = log(cos_i / cos_z) */
+                            /* cos_z is constant then m not changes */
+                            tx = log(cos_i);
+			    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;
+		    }
+		}
+		pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
+		pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
+	    }
+	}
+	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; /* Richter changes to m/a */
+	kk = 1.;
+	G_message("C-factor constant = %lf (a=%.4f; m=%.4f)", 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);
+
+	Rast_get_row(band->fd, band->rast, row, band->type);
+	Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
+
+	pref = band->rast;
+	pcos = cosi->rast;
+	
+	Rast_set_null_value(out->rast, ncols, DCELL_TYPE);
+
+	for (col = 0; col < ncols; col++) {
+
+	    cos_i = Rast_get_d_value(pcos, cosi->type);
+
+	    if (!Rast_is_null_value(pref, band->type) &&
+		!Rast_is_null_value(pcos, cosi->type)) {
+
+		ref_i = Rast_get_d_value(pref, band->type);
+		((DCELL *) out->rast)[col] =
+		    (DCELL) (ref_i * pow((cos_z + cka) / (cos_i + ckb), kk));
+		G_debug(3,
+			"Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
+			ref_i, cka, cos_i, ckb, kk,
+			((DCELL *) out->rast)[col]);
+	    }
+	    pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
+	    pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
+	}
+	Rast_put_row(out->fd, out->rast, out->type);
+    }
+}

imagery/i.topo.corr/i.topo.corr.html

@@ -0,0 +1,118 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.topo.corr</em> is used to topographically correct reflectance
+from imagery files, e.g. obtained with <em>i.landsat.toar</em>, using 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. It can be obtained with <em>r.sun</em>
+(parameter incidout), and then calculating its cosine with float precision.
+
+<p>
+Using flag -i and given an elevation map as basemap (metric),
+<em>i.topo.corr</em> permits to obtain a simple illumination model from the formula:
+<ul>
+	<li> cos_i = cos(s) * cos(z) + sin(s) * sin(z) * cos(a - o) </li>
+</ul>
+where,
+<em>s</em> is 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>
+For each band file, the corrected reflectance (ref_c) is calculate from
+the original reflectance (ref_o) using one of the four offered methods
+<!-- TODO: fix next numbers -->
+(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>
+where,
+<em>k</em> is obtained by linear regression of<br>
+ln(ref_o) = ln(ref_c) - k ln(cos_i/cos_z)
+
+<h3>Method: c-factor</h3>
+
+<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
+
+<h3>Method: percent</h3>
+
+We can use cos_i to estimate the percent of solar incidence on the surface,
+then the transformation (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 the corrected reflectance can
+be calculated as
+<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 (null remain
+    null) and its resolution.</li>
+<li>The elevation map to calculate the illumination model should be metric.</li>
+</ol>
+
+<h2>EXAMPLES</h2>
+
+First, make a illumination model from the elevation map (here, SRTM). Then
+make perform the topographic correction of e.g. 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=SRTM.illumination
+i.topo.corr base=SRTM.illumination 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, 2004. Topographic Correction For Differential
+    Illumination Effects On Ikonos Satellite Imagery. International Archives of
+    Photogrammetry Remote Sensing and Spatial Information, pp. 641-646.</li>
+<li>Meyer, P. and Itten, K.I. and Kellenberger, KJ and Sandmeier, S. and
+    Sandmeier, R., 1993. Radiometric corrections of topographically induced
+    effects on Landsat TM data in alpine terrain. Photogrammetric Engineering
+    and Remote Sensing 48(17).</li>
+<li>Ria&ntilde;o, D. and Chuvieco, E. and Salas, J. and Aguado, I., 2003.
+    Assessment of Different Topographic Corrections in Landsat-TM
+    Data for Mapping Vegetation Types. IEEE Transactions On Geoscience
+    And Remote Sensing, Vol. 41, No. 5</li>
+<li>Twele A. and Erasmi S, 2005. Evaluating topographic correction algorithms
+    for improved land cover discrimination in mountainous areas of
+    Central Sulawesi. G&ouml;ttinger Geographische Abhandlungen, vol. 113.</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.html">r.sun</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+E. Jorge Tizado  (ej.tizado unileon es)<br>
+Dept. Biodiversity and Environmental Management, University of Le&oacute;n, Spain
+
+<p>
+<i>Last changed: $Date$</i>

imagery/i.topo.corr/illumination.c

@@ -0,0 +1,113 @@
+/* File: illumination.c
+ *
+ *  AUTHOR:    E. Jorge Tizado, Spain 2010
+ *
+ *  COPYRIGHT: (c) 2007-10 E. Jorge Tizado
+ *             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 <math.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+void eval_cosi(Gfile * out, Gfile * dem, double zenith, double azimuth)
+{
+    struct Cell_head window;
+
+    int 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;
+
+    Rast_get_window(&window);
+
+    G_begin_distance_calculations();
+    north = Rast_row_to_northing(0.5, &window);
+    center = Rast_row_to_northing(1.5, &window);
+    south = Rast_row_to_northing(2.5, &window);
+    east = Rast_col_to_easting(2.5, &window);
+    west = Rast_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 = Rast_window_rows();
+    ncols = Rast_window_cols();
+
+    cell[0] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
+    Rast_set_d_null_value(cell[0], ncols);
+    cell[1] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
+    Rast_set_d_null_value(cell[1], ncols);
+    cell[2] = (DCELL *) G_calloc(ncols + 1, sizeof(DCELL));
+    Rast_set_d_null_value(cell[2], ncols);
+
+    /* First row is null */
+    Rast_set_null_value((DCELL *) out->rast, Rast_window_cols(), DCELL_TYPE);
+    Rast_put_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;
+	Rast_get_d_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 (Rast_is_d_null_value(c1) || Rast_is_d_null_value(c2) ||
+		Rast_is_d_null_value(c3) || Rast_is_d_null_value(c4) ||
+		Rast_is_d_null_value(c5) || Rast_is_d_null_value(c6) ||
+		Rast_is_d_null_value(c7) || Rast_is_d_null_value(c8) ||
+		Rast_is_d_null_value(c9)) {
+		Rast_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;
+	    }
+	}
+
+	Rast_put_row(out->fd, out->rast, DCELL_TYPE);
+    }
+    /* Last row is null */
+    Rast_set_null_value((DCELL *) out->rast, Rast_window_cols(), DCELL_TYPE);
+    Rast_put_row(out->fd, out->rast, DCELL_TYPE);
+}
+

imagery/i.topo.corr/local_proto.h

@@ -0,0 +1,28 @@
+#include <grass/raster.h>
+
+#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[GNAME_MAX];
+    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

imagery/i.topo.corr/main.c

@@ -0,0 +1,209 @@
+
+/****************************************************************************
+ *
+ * MODULE:       i.topo.corr
+ *
+ * AUTHOR(S):    E. Jorge Tizado - ej.tizado@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 <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+int main(int argc, char *argv[])
+{
+    struct History history;
+    struct GModule *module;
+    struct Cell_head hd_band, hd_dem, window;
+
+    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();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("terrain"));
+    G_add_keyword(_("topographic correction"));
+    module->description = _("Computes topographic correction of reflectance.");
+
+    /* It defines the different parameters */
+
+    input = G_define_standard_option(G_OPT_R_INPUTS);
+    input->required = NO;
+    input->multiple = YES;
+    input->description =
+	_("Name of reflectance raster maps to be corrected topographically");
+
+    output = G_define_standard_option(G_OPT_R_OUTPUT);
+    output->description =
+	_("Name (flag -i) or prefix for output raster maps");
+
+    base = G_define_standard_option(G_OPT_R_MAP);
+    base->key = "basemap";
+    base->description = _("Name of input base raster map (elevation or illumination)");
+
+    zeni = G_define_option();
+    zeni->key = "zenith";
+    zeni->type = TYPE_DOUBLE;
+    zeni->required = YES;
+    zeni->description = _("Solar zenith in degrees");
+
+    azim = G_define_option();
+    azim->key = "azimuth";
+    azim->type = TYPE_DOUBLE;
+    azim->required = NO;
+    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->description = _("Topographic correction method");
+    metho->answer = "c-factor";
+
+    ilum = G_define_flag();
+    ilum->key = 'i';
+    ilum->description = _("Output sun illumination terrain model");
+    
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    if (ilum->answer && azim->answer == NULL)
+	G_fatal_error(_("Solar azimuth is necessary to calculate illumination terrain model"));
+
+    if (!ilum->answer && input->answer == NULL)
+	G_fatal_error
+	    (_("Reflectance maps are necessary to make topographic correction"));
+
+    zenith = atof(zeni->answer);
+    out.type = DCELL_TYPE;
+
+    /* 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) {
+	Rast_get_window(&window);
+	azimuth = atof(azim->answer);
+	/* Warning: make buffers and output after set window */
+	strcpy(dem.name, base->answer);
+	/* Set window to DEM file */
+	Rast_get_window(&window);
+	Rast_get_cellhd(dem.name, "", &hd_dem);
+	Rast_align_window(&window, &hd_dem);
+	dem.fd = Rast_open_old(dem.name, "");
+	dem.type = Rast_get_map_type(dem.fd);
+	/* Open and buffer of the output file */
+	strcpy(out.name, output->answer);
+	out.fd = Rast_open_new(output->answer, DCELL_TYPE);
+	out.rast = Rast_allocate_buf(out.type);
+	/* Open and buffer of the elevation file */
+	dem.rast = Rast_allocate_buf(dem.type);
+	eval_cosi(&out, &dem, zenith, azimuth);
+	/* Close files, buffers, and write history */
+	G_free(dem.rast);
+	Rast_close(dem.fd);
+	G_free(out.rast);
+	Rast_close(out.fd);
+	Rast_short_history(out.name, "raster", &history);
+	Rast_command_history(&history);
+	Rast_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);
+
+	dem.fd = Rast_open_old(base->answer, "");
+	dem.type = Rast_get_map_type(dem.fd);
+	Rast_close(dem.fd);
+	if (dem.type == CELL_TYPE)
+	    G_fatal_error(_("Illumination model is of CELL type"));
+
+	for (i = 0; input->answers[i] != NULL; i++) {
+	    G_message(_("Band %s: "), input->answers[i]);
+	    /* Abre fichero de bandas y el de salida */
+	    strcpy(band.name, input->answers[i]);
+	    Rast_get_cellhd(band.name, "", &hd_band);
+	    Rast_set_window(&hd_band);	/* Antes de out_open y allocate para mismo tamaño */
+	    band.fd = Rast_open_old(band.name, "");
+	    band.type = Rast_get_map_type(band.fd);
+	    if (band.type != DCELL_TYPE) {
+		G_warning(_("Reflectance of <%s> is not of DCELL type - ignored."),
+			  input->answers[i]);
+		Rast_close(band.fd);
+		continue;
+	    }
+	    /* ----- */
+	    dem.fd = Rast_open_old(base->answer, "");
+	    snprintf(out.name, GNAME_MAX - 1, "%s.%s", output->answer,
+		     input->answers[i]);
+	    out.fd = Rast_open_new(out.name, DCELL_TYPE);
+	    out.rast = Rast_allocate_buf(out.type);
+	    band.rast = Rast_allocate_buf(band.type);
+	    dem.rast = Rast_allocate_buf(dem.type);
+	    /* ----- */
+	    eval_tcor(method, &out, &dem, &band, zenith);
+	    /* ----- */
+	    G_free(dem.rast);
+	    Rast_close(dem.fd);
+	    G_free(band.rast);
+	    Rast_close(band.fd);
+	    G_free(out.rast);
+	    Rast_close(out.fd);
+	    Rast_short_history(out.name, "raster", &history);
+	    Rast_command_history(&history);
+	    Rast_write_history(out.name, &history);
+
+	    {
+		struct FPRange range;
+		DCELL min, max;
+		struct Colors grey;
+
+		Rast_read_fp_range(out.name, G_mapset(), &range);
+		Rast_get_fp_range_min_max(&range, &min, &max);
+		Rast_make_grey_scale_colors(&grey, min, max);
+		Rast_write_colors(out.name, G_mapset(), &grey);
+	    }
+	}
+    }
+
+    exit(EXIT_SUCCESS);
+}

imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh

@@ -0,0 +1,108 @@
+#!/bin/sh
+
+# Script to test i.topo.corr with a synthetic map
+#
+# Use North Carolina location to test:
+#   grass64 ~/grassdata/nc_spm_08/user1
+
+if test "$GISBASE" = ""; then
+ echo "You must be in GRASS to run this program."
+ exit
+fi
+
+export GRASS_OVERWRITE=1
+
+# we use the NC location, time zone UTM-5
+g.region n=308500 s=215000 w=630000 e=704800 nsres=250 ewres=250 -pa
+
+# note: no daylight saving time in summer 1954!
+DATETIME=`echo "2 Aug 1954 09:34:53 -0500"`
+
+DAY=`echo $DATETIME | cut -d' ' -f1 | awk '{printf "%d", $1}'`
+# need to translate word month to numeric month for r.sunmask
+MON=`echo $DATETIME | cut -d' ' -f2 | sed 's+Jan+1+g' | sed 's+Feb+2+g' | sed 's+Mar+3+g' | sed 's+Apr+4+g' | sed 's+May+5+g' | sed 's+Jun+6+g' | sed 's+Jul+7+g' | sed 's+Aug+8+g' |sed 's+Sep+9+g' | sed 's+Oct+10+g' | sed 's+Nov+11+g' | sed 's+Dec+12+g'`
+YEAR=`echo $DATETIME | cut -d' ' -f3 | awk '{printf "%d", $1}'`
+TMPTIME=`echo $DATETIME | cut -d' ' -f4 | awk '{printf "%d", $1}'`
+HOUR=`echo $TMPTIME | cut -d':' -f1 | awk '{printf "%d", $1}'`
+MIN=`echo $TMPTIME | cut -d':' -f2 | awk '{printf "%d", $1}'`
+SEC=`echo $TMPTIME | cut -d':' -f3 | awk '{printf "%d", $1}'`
+TIMEZ=`echo $DATETIME | cut -d' ' -f5 | awk '{printf "%d", $1/100}'`
+unset TMPTIME
+
+# create synthetic DEM (kind of roof)
+r.plane --o myplane0 dip=45 az=0 east=637500 north=221750 elev=1000 type=float
+r.plane --o myplane90 dip=45 az=90 east=684800 north=221750 elev=1000 type=float
+r.plane --o myplane180 dip=45 az=180 east=684800 north=260250 elev=1000 type=float
+r.plane --o myplane270 dip=45 az=270 east=684800 north=221750 elev=1000 type=float
+r.mapcalc "myplane_pyr = double(min(myplane90,myplane270,myplane0,myplane180)/10. + 8600.)"
+
+# nviz
+# nviz myplane_pyr
+
+# get sun position
+eval `r.sunmask -s -g output=dummy elev=myplane_pyr year=$YEAR month=8 day=$DAY hour=$HOUR minute=$MIN second=$SEC timezone=$TIMEZ`
+
+solarzenith=`echo $sunangleabovehorizon | awk '{printf "%f", 90. - $1}'`
+echo "Sun position ($DATETIME): solarzenith: $solarzenith, sunazimuth: $sunazimuth"
+
+# shade relief
+r.shaded.relief input=myplane_pyr output=myplane_pyr_shaded altitude=$sunangleabovehorizon azimuth=$sunazimuth
+d.mon stop=wx0 2> /dev/null
+d.mon wx0
+d.rast myplane_pyr_shaded
+
+# pre-run: illumination map
+i.topo.corr -i input=myplane_pyr output=myplane_pyr_illumination \
+	    basemap=myplane_pyr zenith=$solarzenith azimuth=$sunazimuth 
+r.colors myplane_pyr_illumination color=gyr
+
+# show original
+d.erase -f
+r.colors myplane_pyr color=grey
+d.rast.leg myplane_pyr
+echo "Original" | d.text color=black
+
+# nviz with illumination draped over
+# nviz myplane_pyr color=myplane_pyr_illumination
+
+# making the 'band' reflectance file from the shade map
+r.mapcalc "myplane_pyr_band = double((myplane_pyr_shaded - 60.)/18.)"
+r.colors myplane_pyr_band color=gyr
+d.mon stop=wx0 2> /dev/null
+d.mon wx0
+d.rast myplane_pyr_band
+d.legend myplane_pyr_band
+echo "Band reflectance" | d.text color=black
+
+## test it:
+# percent
+METHOD=percent
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+d.mon stop=wx1 2> /dev/null
+d.mon wx1
+d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
+echo "METHOD=percent" | d.text color=black
+
+# minnaert
+METHOD=minnaert
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+d.mon stop=wx2 2> /dev/null
+d.mon wx2
+d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
+echo "METHOD=minnaert" | d.text color=black
+
+# c-factor
+METHOD=c-factor
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+d.mon stop=wx3 2> /dev/null
+d.mon wx3
+d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
+echo "METHOD=c-factor" | d.text color=black
+
+# cosine
+METHOD=cosine
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+d.mon stop=wx4 2> /dev/null
+d.mon wx4
+d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
+echo "METHOD=cosine" | d.text color=black