[GRASS-SVN] r51454 - grass/trunk/imagery/i.topo.corr
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 17 08:48:46 EDT 2012
Author: mmetz
Date: 2012-04-17 05:48:46 -0700 (Tue, 17 Apr 2012)
New Revision: 51454
Modified:
grass/trunk/imagery/i.topo.corr/correction.c
grass/trunk/imagery/i.topo.corr/local_proto.h
grass/trunk/imagery/i.topo.corr/main.c
Log:
add scale option to i.topo.corr
Modified: grass/trunk/imagery/i.topo.corr/correction.c
===================================================================
--- grass/trunk/imagery/i.topo.corr/correction.c 2012-04-16 20:31:58 UTC (rev 51453)
+++ grass/trunk/imagery/i.topo.corr/correction.c 2012-04-17 12:48:46 UTC (rev 51454)
@@ -11,6 +11,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
+#include <float.h>
#include <unistd.h>
#include <grass/gis.h>
#include <grass/raster.h>
@@ -19,20 +20,25 @@
#include "local_proto.h"
void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
- double zenith)
+ double zenith, int do_scale)
{
int row, col, nrows, ncols;
void *pref, *pcos;
- double cos_z, cos_i, ref_i;
+ double cos_z, cos_i, ref_i, result;
double n, sx, sxx, sy, sxy, tx, ty;
double a, m, cka, ckb, kk;
+ double imin, imax, omin, omax, factor; /* for scaling to input */
nrows = Rast_window_rows();
ncols = Rast_window_cols();
cos_z = cos(D2R * zenith);
+
+ imin = omin = DBL_MAX;
+ imax = omax = -DBL_MAX;
+ factor = 1;
/* Calculating regression */
if (method > NON_LAMBERTIAN) {
@@ -53,6 +59,12 @@
if (!Rast_is_null_value(pref, band->type) &&
!Rast_is_null_value(pcos, cosi->type)) {
ref_i = Rast_get_d_value(pref, band->type);
+
+ if (imin > ref_i)
+ imin = ref_i;
+ if (imax < ref_i)
+ imax = ref_i;
+
switch (method) {
case MINNAERT:
if (cos_i > 0. && cos_z > 0. && ref_i > 0.) {
@@ -106,6 +118,45 @@
cka = ckb = 0.;
kk = 1.;
}
+
+ if (do_scale) {
+ /* Topographic correction, pass 1 */
+ 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);
+ result = (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, result);
+
+ if (omin > result)
+ omin = result;
+ if (omax < result)
+ omax = result;
+
+ }
+ pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
+ pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
+ }
+ }
+ G_percent(1, 1, 1);
+ factor = (imax - imin) / (omax - omin);
+ }
/* Topographic correction */
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
@@ -126,8 +177,13 @@
!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));
+ result = (DCELL) (ref_i * pow((cos_z + cka) /
+ (cos_i + ckb), kk));
+
+ if (do_scale)
+ result = (result - omin) * factor + imin;
+
+ ((DCELL *) out->rast)[col] = result;
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,
@@ -138,4 +194,5 @@
}
Rast_put_row(out->fd, out->rast, out->type);
}
+ G_percent(1, 1, 1);
}
Modified: grass/trunk/imagery/i.topo.corr/local_proto.h
===================================================================
--- grass/trunk/imagery/i.topo.corr/local_proto.h 2012-04-16 20:31:58 UTC (rev 51453)
+++ grass/trunk/imagery/i.topo.corr/local_proto.h 2012-04-17 12:48:46 UTC (rev 51454)
@@ -23,6 +23,6 @@
#define C_CORRECT 12
void eval_cosi(Gfile *, Gfile *, double, double);
-void eval_tcor(int, Gfile *, Gfile *, Gfile *, double);
+void eval_tcor(int, Gfile *, Gfile *, Gfile *, double, int);
#endif
Modified: grass/trunk/imagery/i.topo.corr/main.c
===================================================================
--- grass/trunk/imagery/i.topo.corr/main.c 2012-04-16 20:31:58 UTC (rev 51453)
+++ grass/trunk/imagery/i.topo.corr/main.c 2012-04-17 12:48:46 UTC (rev 51454)
@@ -32,11 +32,12 @@
int i;
struct Option *base, *output, *input, *zeni, *azim, *metho;
- struct Flag *ilum;
+ struct Flag *ilum, *scl;
Gfile dem, out, band;
double zenith, azimuth;
int method = COSINE;
+ int do_scale;
/* initialize GIS environment */
G_gisinit(argv[0]);
@@ -88,6 +89,10 @@
ilum->key = 'i';
ilum->description = _("Output sun illumination terrain model");
+ scl = G_define_flag();
+ scl->key = 's';
+ scl->description = _("Scale output to input and copy color rules");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -100,6 +105,7 @@
zenith = atof(zeni->answer);
out.type = DCELL_TYPE;
+ do_scale = scl->answer;
/* Evaluate only cos_i raster file */
/* i.topo.corr -i out=cosi.on07 base=SRTM_v2 zenith=33.3631 azimuth=59.8897 */
@@ -180,7 +186,7 @@
band.rast = Rast_allocate_buf(band.type);
dem.rast = Rast_allocate_buf(dem.type);
/* ----- */
- eval_tcor(method, &out, &dem, &band, zenith);
+ eval_tcor(method, &out, &dem, &band, zenith, do_scale);
/* ----- */
G_free(dem.rast);
Rast_close(dem.fd);
@@ -196,10 +202,18 @@
struct FPRange range;
DCELL min, max;
struct Colors grey;
+ int make_colors = 1;
- 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);
+ if (do_scale) {
+ if (Rast_read_colors(band.name, "", &grey) >= 0)
+ make_colors = 0;
+ }
+
+ if (make_colors) {
+ 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);
}
}
More information about the grass-commit
mailing list