[GRASS-SVN] r60793 - grass/trunk/raster/r.topmodel
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 11 13:57:21 PDT 2014
Author: hcho
Date: 2014-06-11 13:57:20 -0700 (Wed, 11 Jun 2014)
New Revision: 60793
Modified:
grass/trunk/raster/r.topmodel/file_io.c
grass/trunk/raster/r.topmodel/main.c
grass/trunk/raster/r.topmodel/topmodel.c
Log:
r.topmodel: create properly sorted topidxstats file
Modified: grass/trunk/raster/r.topmodel/file_io.c
===================================================================
--- grass/trunk/raster/r.topmodel/file_io.c 2014-06-11 14:01:54 UTC (rev 60792)
+++ grass/trunk/raster/r.topmodel/file_io.c 2014-06-11 20:57:20 UTC (rev 60793)
@@ -21,8 +21,7 @@
{
char buf[1024];
FILE *fp;
- int i, j;
- double x;
+ int i;
/* Read topographic index statistics file */
if ((fp = fopen(file.topidxstats, "r")) == NULL)
@@ -30,7 +29,6 @@
topidxstats.atb = NULL;
topidxstats.Aatb_r = NULL;
- misc.ncells = 0;
for (i = 0; !feof(fp);) {
double atb;
@@ -43,30 +41,13 @@
topidxstats.Aatb_r = (double *)G_realloc(topidxstats.Aatb_r,
(i + 1) * sizeof(double));
topidxstats.atb[i] = atb;
- topidxstats.Aatb_r[i] = Aatb_r;
- misc.ncells += (int)topidxstats.Aatb_r[i++];
+ topidxstats.Aatb_r[i++] = Aatb_r;
}
}
misc.ntopidxclasses = i;
fclose(fp);
- for (i = 0; i < misc.ntopidxclasses; i++)
- topidxstats.Aatb_r[i] /= (double)misc.ncells;
-
- for (i = 0; i < misc.ntopidxclasses; i++) {
- for (j = i; j < misc.ntopidxclasses; j++) {
- if (topidxstats.atb[i] < topidxstats.atb[j]) {
- x = topidxstats.atb[i];
- topidxstats.atb[i] = topidxstats.atb[j];
- topidxstats.atb[j] = x;
- x = topidxstats.Aatb_r[i];
- topidxstats.Aatb_r[i] = topidxstats.Aatb_r[j];
- topidxstats.Aatb_r[j] = x;
- }
- }
- }
-
/* Read parameters file */
if ((fp = fopen(file.params, "r")) == NULL)
G_fatal_error(_("Unable to open input file <%s>"), file.params);
@@ -225,7 +206,7 @@
ltime->tm_mon++;
if ((fp = fopen(file.output, "w")) == NULL)
- G_fatal_error(_("Unable to open output file <%s>"), file.output);
+ G_fatal_error(_("Unable to create output file <%s>"), file.output);
fprintf(fp, "# r.topmodel output file for %s\n", params.name);
fprintf(fp, "# Run time: %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
@@ -241,7 +222,6 @@
fprintf(fp, "# lambda [ln(m^2)]: Average topographic index\n");
fprintf(fp, "# qss [m/timestep]: Saturated subsurface flow per unit area\n");
fprintf(fp, "# qs0 [m/timestep]: Initial subsurface flow per unit area\n");
- fprintf(fp, "# ncells: Number of valid cells\n");
fprintf(fp, "# ntopidxclasses: Number of topographic index classes\n");
fprintf(fp, "# dt [h]: Time step\n");
fprintf(fp, "# ntimesteps: Number of time steps\n");
@@ -282,7 +262,6 @@
fprintf(fp, "lambda: %10.3e\n", misc.lambda);
fprintf(fp, "qss: %10.3e\n", misc.qss);
fprintf(fp, "qs0: %10.3e\n", misc.qs0);
- fprintf(fp, "ncells: %10d\n", misc.ncells);
fprintf(fp, "ntopidxclasses: %10d\n", misc.ntopidxclasses);
fprintf(fp, "dt: %10.3e\n", input.dt);
fprintf(fp, "ntimesteps: %10d\n", input.ntimesteps);
Modified: grass/trunk/raster/r.topmodel/main.c
===================================================================
--- grass/trunk/raster/r.topmodel/main.c 2014-06-11 14:01:54 UTC (rev 60792)
+++ grass/trunk/raster/r.topmodel/main.c 2014-06-11 20:57:20 UTC (rev 60793)
@@ -145,8 +145,8 @@
ntopidxclasses = atoi(params.ntopidxclasses->answer);
outtopidxstats = params.outtopidxstats->answer;
- if (ntopidxclasses <= 0)
- G_fatal_error(_("%s must be positive"), "ntopidxclasses");
+ if (ntopidxclasses <= 1)
+ G_fatal_error(_("%s must be greater than 1"), "ntopidxclasses");
create_topidxstats(topidx, ntopidxclasses, outtopidxstats);
} else if (params.topidx->answer) {
Modified: grass/trunk/raster/r.topmodel/topmodel.c
===================================================================
--- grass/trunk/raster/r.topmodel/topmodel.c 2014-06-11 14:01:54 UTC (rev 60792)
+++ grass/trunk/raster/r.topmodel/topmodel.c 2014-06-11 20:57:20 UTC (rev 60793)
@@ -6,19 +6,62 @@
void create_topidxstats(char *topidx, int ntopidxclasses, char *outtopidxstats)
{
- char input[GPATH_MAX];
- char output[GPATH_MAX];
- char nsteps[32];
+ char input[GPATH_MAX], nsteps[32];
+ const char *args[5];
+ struct Popen child;
+ FILE *fp;
+ double *atb, *Aatb_r;
+ int i;
+ int total_ncells;
sprintf(input, "input=%s", topidx);
- sprintf(output, "output=%s", outtopidxstats);
- sprintf(nsteps, "nsteps=%d", ntopidxclasses);
+ sprintf(nsteps, "nsteps=%d", ntopidxclasses - 1);
G_message("Creating topographic index statistics file...");
- G_verbose_message("r.stats -Anc %s %s %s ...", input, output, nsteps);
+ G_verbose_message("r.stats -nc %s %s ...", input, nsteps);
- if (G_spawn("r.stats", "r.stats", "-Anc", input, output, nsteps, NULL) != 0)
+ args[0] = "r.stats";
+ args[1] = "-nc";
+ args[2] = input;
+ args[3] = nsteps;
+ args[4] = NULL;
+
+ if ((fp = G_popen_read(&child, "r.stats", args)) == NULL)
G_fatal_error(_("Unable to run %s"), "r.stats");
+
+ atb = (double *)G_malloc(ntopidxclasses * sizeof(double));
+ Aatb_r = (double *)G_malloc(ntopidxclasses * sizeof(double));
+
+ total_ncells = 0;
+ for (i = 0; i < ntopidxclasses - 1 && !feof(fp);) {
+ double atb1, atb2;
+ int ncells;
+
+ get_line(fp, buf);
+ if (sscanf(buf, "%lf-%lf %d", &atb1, &atb2, &ncells) == 3) {
+ atb[i] = atb1;
+ Aatb_r[i] = (double)ncells;
+ total_ncells += ncells;
+
+ if (++i == ntopidxclasses - 1) {
+ atb[i] = atb2;
+ Aatb_r[i] = 0.0;
+ }
+ }
+ }
+
+ G_popen_close(&child);
+
+ if (i < ntopidxclasses - 1)
+ G_fatal_error(_("Invalid %s output"), "r.stats");
+
+ if ((fp = fopen(outtopidxstats, "w")) == NULL)
+ G_fatal_error(_("Unable to create output file <%s>"), outtopidxstats);
+
+ for (i = ntopidxclasses - 1; i >= 0; i--)
+ fprintf(fp, "%10.3e %10.3e\n", atb[i], Aatb_r[i] / total_ncells);
+
+ fclose(fp);
}
/* Calculate the areal average of topographic index */
@@ -38,7 +81,7 @@
/* Initialize the flows */
void initialize(void)
{
- int i, j, t;
+ int i;
double A1, A2;
/* average topographic index */
@@ -81,6 +124,8 @@
/* cumulative ratio of the contribution area for each time step */
misc.Ad = (double *)G_malloc(misc.tcsub * sizeof(double));
for (i = 0; i < misc.tcsub; i++) {
+ int j, t;
+
t = misc.delay + i + 1;
if (t > misc.tch[params.nch - 1])
misc.Ad[i] = 1.0;
@@ -266,9 +311,7 @@
qo = topidxstats.Aatb_r[j] *
(misc.ex[i][j - 1] + misc.ex[i][j]) / 2.0;
else if (misc.ex[i][j - 1] > 0.0)
- qo = Aatb_r * misc.ex[i][j - 1] /
- (misc.ex[i][j - 1] -
- misc.ex[i][j]) * misc.ex[i][j - 1] / 2.0;
+ qo = Aatb_r * misc.ex[i][j - 1] / 2.0;
}
misc.qo[i][j] = qo;
misc.qo[i][misc.ntopidxclasses] += misc.qo[i][j];
More information about the grass-commit
mailing list