r.topmodel: create properly sorted topidxstats file (merge from https://trac.osgeo.org/grass/changeset/60793)

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@60794 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.topmodel/file_io.c

@@ -21,8 +21,7 @@ void read_input(void)
 {
     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 @@ void read_input(void)
 
     topidxstats.atb = NULL;
     topidxstats.Aatb_r = NULL;
-    misc.ncells = 0;
 
     for (i = 0; !feof(fp);) {
 	double atb;
@@ -43,30 +41,13 @@ void read_input(void)
 	    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 @@ void write_output(void)
     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 @@ void write_output(void)
     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 @@ void write_output(void)
     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);

raster/r.topmodel/main.c

@@ -145,8 +145,8 @@ int main(int argc, char **argv)
 	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) {

raster/r.topmodel/topmodel.c

@@ -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);
+
+    args[0] = "r.stats";
+    args[1] = "-nc";
+    args[2] = input;
+    args[3] = nsteps;
+    args[4] = NULL;
 
-    if (G_spawn("r.stats", "r.stats", "-Anc", input, output, nsteps, NULL) != 0)
+    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 @@ double calculate_lambda(void)
 /* Initialize the flows */
 void initialize(void)
 {
-    int i, j, t;
+    int i;
     double A1, A2;
 
     /* average topographic index */
@@ -81,6 +124,8 @@ void initialize(void)
     /* 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 @@ void calculate_flows(void)
 		    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];