Added dt, ntimesteps, and nch in output.

Removed the obsflow parameter and let the user evaluate the model performance.


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

raster/r.topmodel/file_io.c

@@ -1,7 +1,6 @@
 #include <stdlib.h>
 #include <string.h>
 #include <time.h>
-#include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/glocale.h>
 #include "global.h"
@@ -16,8 +15,6 @@ void get_line(FILE * fp, char *buffer)
 
     if ((str = (char *)strchr(buffer, '#')))
 	*str = 0;
-
-    return;
 }
 
 void read_input(void)
@@ -208,29 +205,10 @@ void read_input(void)
     input.ntimesteps = i;
     fclose(fp);
 
-    /* Read obsflow file */
-    if (file.obsflow) {
-	if ((fp = fopen(file.obsflow, "r")) == NULL)
-	    G_fatal_error(_("Unable to open input file <%s>"), file.obsflow);
-
-	misc.Qobs = (double *)G_malloc(input.ntimesteps * sizeof(double));
-
-	for (i = 0; i < input.ntimesteps && !feof(fp);) {
-	    get_line(fp, buf);
-	    if (sscanf(buf, "%lf", &(misc.Qobs[i])) == 1)
-		i++;
-	}
-
-	input.ntimesteps = (i < input.ntimesteps ? i : input.ntimesteps);
-	fclose(fp);
-    }
-
     if (!(misc.timestep > 0 && misc.timestep < input.ntimesteps + 1))
 	misc.timestep = 0;
     if (!(misc.topidxclass > 0 && misc.topidxclass < misc.ntopidxclasses + 1))
 	misc.topidxclass = 0;
-
-    return;
 }
 
 void write_output(void)
@@ -255,113 +233,72 @@ void write_output(void)
 	    ltime->tm_year, ltime->tm_mon, ltime->tm_mday,
 	    ltime->tm_hour, ltime->tm_min, ltime->tm_sec);
     fprintf(fp, "#\n");
-    if (file.obsflow) {
-	fprintf(fp, "# %-15s Model efficiency\n", "Em:");
-	fprintf(fp, "# %-15s Peak observed Q\n"
-		"# %77s\n", "Qobs_peak:", "[m^3/timestep]");
-	fprintf(fp, "# %-15s Peak time for observed Q\n"
-		"# %77s\n", "tobs_peak:", "[timestep]");
-	fprintf(fp, "# %-15s Mean observed Q\n"
-		"# %77s\n", "Qobs_mean:", "[m^3/timestep]");
-    }
-    fprintf(fp, "# %-15s Peak simulated Q\n"
-	    "# %77s\n", "Qt_peak:", "[m^3/timestep]");
-    fprintf(fp, "# %-15s Peak time for simulated Q\n"
-	    "# %77s\n", "tt_peak:", "[timestep]");
-    fprintf(fp, "# %-15s Mean simulated Q\n"
-	    "# %77s\n", "Qt_mean:", "[m^3/timestep]");
-    fprintf(fp, "# %-15s Number of time steps\n", "ntimesteps:");
-    fprintf(fp, "# %-15s Number of non-NULL cells\n", "ncells:");
-    fprintf(fp, "# %-15s Number of topographic index classes\n",
-	    "ntopidxclasses:");
-    fprintf(fp, "# %-15s Number of delay time steps (delay time between "
-	    "rainfall and\n#\t\t\tflow response)\n", "ndelays:");
-    fprintf(fp, "# %-15s Number of reach time steps "
-	    "(time of concentration)\n", "nreaches:");
-    fprintf(fp, "# %-15s Areal average of ln(T0) = ln(Te)\n"
-	    "# %77s\n", "lnTe:", "[ln(m^2/timestep)]");
-    fprintf(fp, "# %-15s Main channel routing velocity\n"
-	    "# %77s\n", "vch:", "[m/timestep]");
-    fprintf(fp, "# %-15s Internal subcatchment routing velocity\n"
-	    "# %77s\n", "vr:", "[m/timestep]");
-    fprintf(fp, "# %-15s Areal average of topographic index\n"
-	    "# %77s\n", "lambda:", "[ln(m^2)]");
-    fprintf(fp, "# %-15s Subsurface flow per unit area at a soil surface\n"
-	    "# %77s\n", "qss:", "[m/timestep]");
-    fprintf(fp, "# %-15s Initial subsurface flow per unit area\n"
-	    "# %77s\n", "qs0:", "[m/timestep]");
+    fprintf(fp, "# Qt_peak [m^3/timestep]:  Total peak flow\n");
+    fprintf(fp, "# tt_peak [timestep]:      Peak time for total flow\n");
+    fprintf(fp, "# Qt_mean [m^3/timestep]:  Mean total flow\n");
+    fprintf(fp, "# lnTe [ln(m^2/timestep)]: Areal average of ln(T0)\n");
+    fprintf(fp, "# vch [m/timestep]:        Main channel routing velocity\n");
+    fprintf(fp, "# vr [m/timestep]:         Internal subcatchment routing velocity\n");
+    fprintf(fp, "# lambda [ln(m^2)]:        Areal average of topographic index\n");
+    fprintf(fp, "# qss [m/timestep]:        Subsurface flow per unit area at a soil surface\n");
+    fprintf(fp, "# qs0 [m/timestep]:        Initial subsurface flow per unit area\n");
+    fprintf(fp, "# ncells:                  Number of non-NULL 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");
+    fprintf(fp, "# nch:                     Number of routing time steps\n");
+    fprintf(fp, "# nreaches:                Number of reach time steps (time of concentration)\n");
+    fprintf(fp, "# ndelays:                 Delay time steps between rainfall and flow response\n");
     fprintf(fp, "#\n");
-    fprintf(fp, "# %-15s Routing time step\n"
-	    "# %77s\n", "tch:", "[timestep]");
-    fprintf(fp, "# %-15s Difference in contribution area for each reach "
-	    "time step\n" "# %77s\n", "Ad:", "[m^2]");
-    fprintf(fp, "# %-15s Total flow\n" "# %77s\n", "Qt:", "[m^3/timestep]");
-    fprintf(fp, "# %-15s Total flow per unit area\n"
-	    "# %77s\n", "qt:", "[m/timestep]");
-    fprintf(fp, "# %-15s Saturation overland flow per unit area\n"
-	    "# %77s\n", "qo:", "[m/timestep]");
-    fprintf(fp, "# %-15s Subsurface flow per unit area\n"
-	    "# %77s\n", "qs:", "[m/timestep]");
-    fprintf(fp, "# %-15s Vertical flux (or drainage flux)\n"
-	    "# %77s\n", "qv:", "[m/timestep]");
-    fprintf(fp, "# %-15s Mean saturation deficit in the watershed\n"
-	    "# %77s\n", "S_mean:", "[m]");
+    fprintf(fp, "# tch [timestep]:          Routing time\n");
+    fprintf(fp, "# Ad [m^2]:                Difference in contribution area for reach time steps\n");
+    fprintf(fp, "# Qt [m^3/timestep]:       Total flow\n");
+    fprintf(fp, "# qt [m/timestep]:         Total flow per unit area\n");
+    fprintf(fp, "# qo [m/timestep]:         Saturation overland flow per unit area\n");
+    fprintf(fp, "# qs [m/timestep]:         Subsurface flow per unit area\n");
+    fprintf(fp, "# qv [m/timestep]:         Vertical flux (or drainage flux)\n");
+    fprintf(fp, "# S_mean [m]:              Mean saturation deficit in the watershed\n");
     if (params.infex) {
-	fprintf(fp, "# %-15s Infiltration rate\n"
-		"# %30s\n", "f:", "[m/timestep]");
-	fprintf(fp, "# %-15s Infiltration excess runoff\n"
-		"# %77s\n", "fex:", "[m/timestep]");
+	fprintf(fp, "# f [m/timestep]:          Infiltration rate\n");
+	fprintf(fp, "# fex [m/timestep]:        Infiltration excess runoff\n");
     }
 
     if (misc.timestep || misc.topidxclass) {
 	fprintf(fp, "#\n");
-	fprintf(fp, "# %-15s Root zone storage deficit\n"
-		"# %77s\n", "Srz:", "[m]");
-	fprintf(fp, "# %-15s Unsaturated (gravity drainage) zone "
-		"storage\n" "# %77s\n", "Suz:", "[m]");
-	fprintf(fp, "# %-15s Local saturated zone deficit due to "
-		"gravity drainage\n" "# %77s\n", "S:", "[m]");
-	fprintf(fp, "# %-15s Actual evapotranspiration\n"
-		"# %77s\n", "Ea:", "[m/timestep]");
-	fprintf(fp, "# %-15s Excess flow from a fully saturated "
-		"area per unit area\n" "# %77s\n", "ex:", "[m/timestep]");
+	fprintf(fp, "# Srz [m]:                 Root zone storage deficit\n");
+	fprintf(fp, "# Suz [m]:                 Unsaturated (gravity drainage) zone storage\n");
+	fprintf(fp, "# S [m]:                   Local saturated zone deficit due to gravity drainage\n");
+	fprintf(fp, "# Ea [m/timestep]:         Actual evapotranspiration\n");
+	fprintf(fp, "# ex [m/timestep]:         Excess flow from a fully saturated area per unit area\n");
     }
     fprintf(fp, "\n");
 
-    if (file.obsflow) {
-	fprintf(fp, "%-16s ", "Em:");
-	if (!Rast_is_d_null_value(&misc.Em))
-	    fprintf(fp, "%10.5lf\n", misc.Em);
-	else
-	    fprintf(fp, "Not resolved due to constant observed Q\n");
-	fprintf(fp, "%-16s %10.3le\n", "Qobs_peak:", misc.Qobs_peak);
-	fprintf(fp, "%-16s %10d\n", "tobs_peak:", misc.tobs_peak);
-	fprintf(fp, "%-16s %10.3le\n", "Qobs_mean:", misc.Qobs_mean);
-    }
-    fprintf(fp, "%-16s %10.3le\n", "Qt_peak:", misc.Qt_peak);
-    fprintf(fp, "%-16s %10d\n", "tt_peak:", misc.tt_peak);
-    fprintf(fp, "%-16s %10.3le\n", "Qt_mean:", misc.Qt_mean);
-    fprintf(fp, "%-16s %10d\n", "ntimesteps:", input.ntimesteps);
-    fprintf(fp, "%-16s %10d\n", "ncells:", misc.ncells);
-    fprintf(fp, "%-16s %10d\n", "ntopidxclasses:", misc.ntopidxclasses);
-    fprintf(fp, "%-16s %10d\n", "ndelays:", misc.ndelays);
-    fprintf(fp, "%-16s %10d\n", "nreaches:", misc.nreaches);
-    fprintf(fp, "%-16s %10.3le\n", "lnTe:", misc.lnTe);
-    fprintf(fp, "%-16s %10.3le\n", "vch:", misc.vch);
-    fprintf(fp, "%-16s %10.3le\n", "vr:", misc.vr);
-    fprintf(fp, "%-16s %10.3le\n", "lambda:", misc.lambda);
-    fprintf(fp, "%-16s %10.3le\n", "qss:", misc.qss);
-    fprintf(fp, "%-16s %10.3le\n", "qs0:", misc.qs0);
+    fprintf(fp, "Qt_peak:        %10.3e\n", misc.Qt_peak);
+    fprintf(fp, "tt_peak:        %10d\n", misc.tt_peak);
+    fprintf(fp, "Qt_mean:        %10.3e\n", misc.Qt_mean);
+    fprintf(fp, "lnTe:           %10.3e\n", misc.lnTe);
+    fprintf(fp, "vch:            %10.3e\n", misc.vch);
+    fprintf(fp, "vr:             %10.3e\n", misc.vr);
+    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);
+    fprintf(fp, "nch:            %10d\n", params.nch);
+    fprintf(fp, "nreaches:       %10d\n", misc.nreaches);
+    fprintf(fp, "ndelays:        %10d\n", misc.ndelays);
     fprintf(fp, "\n");
 
     fprintf(fp, "%10s\n", "tch");
     for (i = 0; i < params.nch; i++)
-	fprintf(fp, "%10.3le\n", misc.tch[i]);
+	fprintf(fp, "%10.3e\n", misc.tch[i]);
     fprintf(fp, "\n");
 
     fprintf(fp, "%10s\n", "Ad");
     for (i = 0; i < misc.nreaches; i++)
-	fprintf(fp, "%10.3le\n", misc.Ad[i]);
+	fprintf(fp, "%10.3e\n", misc.Ad[i]);
     fprintf(fp, "\n");
 
     st = et = si = ei = 0;
@@ -385,20 +322,20 @@ void write_output(void)
 	}
     }
 
-    fprintf(fp, "%10s %10s %10s %10s %10s %10s %10s", "timestep", "Qt",
-	    "qt", "qo", "qs", "qv", "S_mean");
+    fprintf(fp, "%10s %10s %10s %10s %10s %10s %10s",
+	   "timestep", "Qt", "qt", "qo", "qs", "qv", "S_mean");
     if (params.infex)
 	fprintf(fp, " %10s %10s", "f", "fex");
     fprintf(fp, "\n");
 
     for (i = 0; i < input.ntimesteps; i++) {
-	fprintf(fp, "%10d %10.3le %10.3le %10.3le %10.3le %10.3le "
-		"%10.3le", i + 1, misc.Qt[i],
+	fprintf(fp, "%10d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e",
+		i + 1, misc.Qt[i],
 		misc.qt[i][misc.ntopidxclasses],
 		misc.qo[i][misc.ntopidxclasses], misc.qs[i],
 		misc.qv[i][misc.ntopidxclasses], misc.S_mean[i]);
 	if (params.infex)
-	    fprintf(fp, " %10.3le %10.3le", misc.f[i], misc.fex[i]);
+	    fprintf(fp, " %10.3e %10.3e", misc.f[i], misc.fex[i]);
 	fprintf(fp, "\n");
     }
 
@@ -432,9 +369,9 @@ void write_output(void)
 		    fprintf(fp, "%10d ", i + 1);
 		}
 
-		fprintf(fp, "%10.3le %10.3le %10.3le "
-			"%10.3le %10.3le %10.3le "
-			"%10.3le %10.3le %10.3le\n",
+		fprintf(fp, "%10.3e %10.3e %10.3e "
+			"%10.3e %10.3e %10.3e "
+			"%10.3e %10.3e %10.3e\n",
 			misc.qt[i][j], misc.qo[i][j],
 			misc.qs[i], misc.qv[i][j],
 			misc.Srz[i][j], misc.Suz[i][j],
@@ -443,6 +380,4 @@ void write_output(void)
     }
 
     fclose(fp);
-
-    return;
 }

raster/r.topmodel/global.h

@@ -80,7 +80,6 @@ struct file
     char *topidxstats;
     char *input;
     char *output;
-    char *obsflow;
 };
 
 /* Miscellaneous TOPMODEL variables */
@@ -90,8 +89,6 @@ struct misc
     int ncells;
     /* Number of topographic index classes */
     int ntopidxclasses;
-    /* Model efficiency */
-    double Em;
     int ndelays;
     int nreaches;
     double lnTe;
@@ -100,9 +97,7 @@ struct misc
     double lambda;
     double qss;
     double qs0;
-    double Qobs_peak;
     double Qt_peak;
-    double Qobs_mean;
     double Qt_mean;
     int tobs_peak;
     int tt_peak;
@@ -111,7 +106,6 @@ struct misc
     /* misc.nreach's */
     double *Ad;
     /* input.ntimestep's */
-    double *Qobs;
     double *Qt;
     double *qs;			/* spatially constant? */
     double *S_mean;

raster/r.topmodel/main.c

@@ -30,7 +30,6 @@ int main(int argc, char **argv)
 	struct Option *topidxstats;
 	struct Option *input;
 	struct Option *output;
-	struct Option *obsflow;
 	struct Option *timestep;
 	struct Option *topidxclass;
 	struct Option *topidx;
@@ -68,11 +67,6 @@ int main(int argc, char **argv)
     params.output = G_define_standard_option(G_OPT_F_OUTPUT);
     params.output->description = _("Name for output file");
 
-    params.obsflow = G_define_standard_option(G_OPT_F_INPUT);
-    params.obsflow->key = "obsflow";
-    params.obsflow->description = _("Name of observed flow file");
-    params.obsflow->required = NO;
-
     params.timestep = G_define_option();
     params.timestep->key = "timestep";
     params.timestep->label = _("Time step");
@@ -130,7 +124,6 @@ int main(int argc, char **argv)
     file.topidxstats = params.topidxstats->answer;
     file.input = params.input->answer;
     file.output = params.output->answer;
-    file.obsflow = params.obsflow->answer;
 
     if (!params.timestep->answer)
 	params.timestep->answer = "0";

raster/r.topmodel/r.topmodel.html

@@ -101,22 +101,6 @@ This file contains observed weather data.
 </pre></div>
 </dd>
 
-<dt><b>obsflow</b></dt>
-<dd>
-Compare simulated flows with observed flows and calculate the model
-efficiency. This file contains observed flow data and the number of records
-should match the number of time steps (ntimesteps in the input file).
-<div class="code"><pre>
-# Qobs [m^3/dt]: Observed flow per time step
-2568918.24
-1668573.562
-1039800.24
-.
-.
-.
-</pre></div>
-</dd>
-
 <dt><b>timestep</b></dt>
 <dd>
 If a time step is specified, output will be generated for the specific time

raster/r.topmodel/topmodel.c

@@ -1,5 +1,4 @@
 #include <math.h>
-#include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/spawn.h>
 #include "global.h"
@@ -24,16 +23,14 @@ void create_topidxstats(char *topidx, int ntopidxclasses, char *outtopidxstats)
 double calculate_lambda(void)
 {
     int i;
-    double retval;
+    double lambda;
 
-
-    retval = 0.0;
+    lambda = 0.0;
     for (i = 1; i < misc.ntopidxclasses; i++)
-	retval += topidxstats.Aatb_r[i] *
+	lambda += topidxstats.Aatb_r[i] *
 	    (topidxstats.atb[i] + topidxstats.atb[i - 1]) / 2.0;
 
-
-    return retval;
+    return lambda;
 }
 
 void initialize(void)
@@ -41,7 +38,6 @@ void initialize(void)
     int i, j, t;
     double A1, A2;
 
-
     misc.lambda = calculate_lambda();
     misc.lnTe = params.lnTe + log(input.dt);
     misc.vch = params.vch * input.dt;
@@ -117,9 +113,6 @@ void initialize(void)
 	A1 += misc.Ad[i];
 	misc.Qt[misc.ndelays + i] = misc.qs0 * (params.A - A1);
     }
-
-
-    return;
 }
 
 void calculate_flows(void)
@@ -129,7 +122,6 @@ void calculate_flows(void)
     double R;
     double qo, qv;
 
-
     misc.S = (double **)G_malloc(input.ntimesteps * sizeof(double *));
     misc.Ea = (double **)G_malloc(input.ntimesteps * sizeof(double *));
     misc.ex = (double **)G_malloc(input.ntimesteps * sizeof(double *));
@@ -261,45 +253,6 @@ void calculate_flows(void)
 	}
     }
 
-
-    return;
-}
-
-/* Objective function for hydrograph suggested by Servet and Dezetter(1991) */
-double calculate_efficiency(void)
-{
-    int i;
-    double Em, numerator, denominator;
-
-
-    misc.Qobs_mean = 0.0;
-    numerator = 0.0;
-    for (i = 0; i < input.ntimesteps; i++) {
-	misc.Qobs_mean += misc.Qobs[i];
-	numerator += pow(misc.Qobs[i] - misc.Qt[i], 2.0);
-    }
-    misc.Qobs_mean /= input.ntimesteps;
-
-    denominator = 0.0;
-    for (i = 0; i < input.ntimesteps; i++)
-	denominator += pow(misc.Qobs[i] - misc.Qobs_mean, 2.0);
-
-    if (denominator == 0.0) {
-	G_warning("Em can not be resolved due to constant observed Q");
-	Rast_set_d_null_value(&Em, 1);
-    }
-    else {
-	Em = 1.0 - numerator / denominator;
-    }
-
-    return Em;
-}
-
-void calculate_others(void)
-{
-    int i;
-
-
     misc.Qt_mean = 0.0;
     for (i = 0; i < input.ntimesteps; i++) {
 	misc.Qt_mean += misc.Qt[i];
@@ -309,25 +262,10 @@ void calculate_others(void)
 	}
     }
     misc.Qt_mean /= input.ntimesteps;
-
-    if (file.obsflow) {
-	misc.Em = calculate_efficiency();
-	for (i = 0; i < input.ntimesteps; i++) {
-	    if (!i || misc.Qobs_peak < misc.Qobs[i]) {
-		misc.Qobs_peak = misc.Qobs[i];
-		misc.tobs_peak = i + 1;
-	    }
-	}
-    }
-
-    return;
 }
 
 void run_topmodel(void)
 {
     initialize();
     calculate_flows();
-    calculate_others();
-
-    return;
 }