[GRASS-SVN] r42708 - grass-addons/imagery/i.topo.corr
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jul 7 04:27:30 EDT 2010
Author: neteler
Date: 2010-07-07 08:27:30 +0000 (Wed, 07 Jul 2010)
New Revision: 42708
Modified:
grass-addons/imagery/i.topo.corr/correction.c
grass-addons/imagery/i.topo.corr/illumination.c
grass-addons/imagery/i.topo.corr/main.c
Log:
GRASS style indentation applied
Modified: grass-addons/imagery/i.topo.corr/correction.c
===================================================================
--- grass-addons/imagery/i.topo.corr/correction.c 2010-07-07 06:21:48 UTC (rev 42707)
+++ grass-addons/imagery/i.topo.corr/correction.c 2010-07-07 08:27:30 UTC (rev 42708)
@@ -7,145 +7,138 @@
#include "local_proto.h"
-void eval_tcor(int method, Gfile *out, Gfile *cosi, Gfile *band, double zenith)
+void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
+ double zenith)
{
- int i, row, col, nrows, ncols;
- void *pref, *pcos;
+ 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;
+ 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();
+ nrows = G_window_rows();
+ ncols = G_window_cols();
- cos_z = cos(D2R * zenith);
+ 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);
+ /* 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);
+ 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);
+ 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;
- }
- }
+ 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;
}
- }
- 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);
+ case C_CORRECT:
+ {
+ n++;
+ sx += cos_i;
+ sxx += cos_i * cos_i;
+ sy += ref_i;
+ sxy += cos_i * ref_i;
+ }
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);
+ 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);
+ 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);
+ 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);
+ 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);
+ }
}
Modified: grass-addons/imagery/i.topo.corr/illumination.c
===================================================================
--- grass-addons/imagery/i.topo.corr/illumination.c 2010-07-07 06:21:48 UTC (rev 42707)
+++ grass-addons/imagery/i.topo.corr/illumination.c 2010-07-07 08:27:30 UTC (rev 42708)
@@ -9,262 +9,279 @@
void eval_c_cosi(Gfile * out, Gfile * dem, double zenith, double azimuth)
{
- struct Cell_head window;
+ 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;
+ 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_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;
+ 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;
+ zenith *= D2R;
+ azimuth *= D2R;
- cos_z = cos(zenith);
- sin_z = sin(zenith);
+ cos_z = cos(zenith);
+ sin_z = sin(zenith);
- /* Making cos_i raster ... */
- nrows = G_window_rows();
- ncols = G_window_cols();
+ /* 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);
+ 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);
+ /* 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;
+ 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;
+ 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);
+ cos_i =
+ cos_z * cos(slope) + sin_z * sin(slope) * cos(azimuth -
+ aspect);
- ((DCELL *)out->rast)[col] = (DCELL)cos_i;
- }
- }
+ ((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);
+ }
+ /* 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;
+ 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;
+ 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_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;
+ 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;
+ zenith *= D2R;
+ azimuth *= D2R;
- cos_z = cos(zenith);
- sin_z = sin(zenith);
+ cos_z = cos(zenith);
+ sin_z = sin(zenith);
- /* Making cos_i raster ... */
- nrows = G_window_rows();
- ncols = G_window_cols();
+ /* 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);
+ 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);
+ /* 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;
+ 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;
+ 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);
+ cos_i =
+ cos_z * cos(slope) + sin_z * sin(slope) * cos(azimuth -
+ aspect);
- ((DCELL *)out->rast)[col] = (DCELL)cos_i;
- }
- }
+ ((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);
+ }
+ /* 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;
+ 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;
+ 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_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;
+ 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;
+ zenith *= D2R;
+ azimuth *= D2R;
- cos_z = cos(zenith);
- sin_z = sin(zenith);
+ cos_z = cos(zenith);
+ sin_z = sin(zenith);
- /* Making cos_i raster ... */
- nrows = G_window_rows();
- ncols = G_window_cols();
+ /* 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);
+ 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);
+ /* 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;
+ 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;
+ 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);
+ cos_i =
+ cos_z * cos(slope) + sin_z * sin(slope) * cos(azimuth -
+ aspect);
- ((DCELL *)out->rast)[col] = (DCELL)cos_i;
- }
- }
+ ((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);
+ }
+ /* 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);
}
-
Modified: grass-addons/imagery/i.topo.corr/main.c
===================================================================
--- grass-addons/imagery/i.topo.corr/main.c 2010-07-07 06:21:48 UTC (rev 42707)
+++ grass-addons/imagery/i.topo.corr/main.c 2010-07-07 08:27:30 UTC (rev 42708)
@@ -1,3 +1,4 @@
+
/****************************************************************************
*
* MODULE: i.topo.corr
@@ -25,231 +26,238 @@
extern void eval_f_cosi(Gfile *, Gfile *, double, double);
extern void eval_d_cosi(Gfile *, Gfile *, double, double);
-int full_open_old(Gfile * gf, char * fname)
+int full_open_old(Gfile * gf, char *fname)
{
- gf->fd = -1;
+ gf->fd = -1;
- snprintf(gf->name, 127, "%s", fname);
+ 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->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->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);
+ gf->type = G_raster_map_type(gf->name, gf->mapset);
- return gf->fd;
+ return gf->fd;
}
-int full_open_new(Gfile * gf, char * fname, RASTER_MAP_TYPE ftype)
+int full_open_new(Gfile * gf, char *fname, RASTER_MAP_TYPE ftype)
{
- gf->fd = -1;
+ 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);
+ 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);
+ 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;
+ return gf->fd;
}
int main(int argc, char *argv[])
{
- struct History history;
- struct GModule *module;
- struct Cell_head hd_band, hd_dem, window;
+ struct History history;
+ struct GModule *module;
+ struct Cell_head hd_band, hd_dem, window;
- char bufname[128];
- int nrows, ncols, row, col;
+ char bufname[128];
+ int nrows, ncols, row, col;
- int i;
- struct Option *base, *output, *input, *zeni, *azim, *metho;
- struct Flag *ilum;
+ int i;
+ struct Option *base, *output, *input, *zeni, *azim, *metho;
+ struct Flag *ilum;
- Gfile dem, out, band;
- double zenith, azimuth;
- int method = COSINE;
+ Gfile dem, out, band;
+ double zenith, azimuth;
+ int method = COSINE;
- /* initialize GIS environment */
- G_gisinit(argv[0]);
+ /* 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");
+ /* 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 */
+ /* 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");
+ 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");
+ 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)");
+ 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");
+ 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)");
+ 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";
+ 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;
+ 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);
+ exit(EXIT_FAILURE);
- if (ilum->answer && azim->answer == NULL)
- G_fatal_error("Solar azimuth is necessary to calculate ilumination terrain model.");
+ 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.");
+ if (!ilum->answer && input->answer == NULL)
+ G_fatal_error
+ ("Reflectance files are necessary to make topographic correction.");
- zenith = atof(zeni->answer);
+ 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 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);
}
- /* 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 (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
- {
-// 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);
+ 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"));
- 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);
- 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];
- char command[300];
- sprintf(command, "r.colors map=%s color=grey", out.name);
- system(command);
- }
- G_close_cell(dem.fd);
+ 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