[GRASS-SVN] r45670 - in grass/trunk/raster/simwe: r.sim.sediment
r.sim.water simlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 15 06:16:37 EDT 2011
Author: huhabla
Date: 2011-03-15 03:16:37 -0700 (Tue, 15 Mar 2011)
New Revision: 45670
Added:
grass/trunk/raster/simwe/simlib/observation_points.c
grass/trunk/raster/simwe/simlib/utils.c
Modified:
grass/trunk/raster/simwe/r.sim.sediment/main.c
grass/trunk/raster/simwe/r.sim.water/main.c
grass/trunk/raster/simwe/r.sim.water/r.sim.water.html
grass/trunk/raster/simwe/r.sim.water/spearfish.sh
grass/trunk/raster/simwe/simlib/hydro.c
grass/trunk/raster/simwe/simlib/input.c
grass/trunk/raster/simwe/simlib/output.c
grass/trunk/raster/simwe/simlib/waterglobs.h
Log:
Reimplemented observation vector point support. Code clean up in many
files and partly rewrite of input.c. Docu update related to the
observation points. Update of the spearfish example in r.sim.water.
Modified: grass/trunk/raster/simwe/r.sim.sediment/main.c
===================================================================
--- grass/trunk/raster/simwe/r.sim.sediment/main.c 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/r.sim.sediment/main.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -168,15 +168,20 @@
_("Base name of the output walkers vector points map");
parm.outwalk->guisection = _("Output_options");
- /* needs to be updated to GRASS 6 vector format !!
- parm.sfile = G_define_standard_option(G_OPT_V_INPUT);
- parm.sfile->key = "vector";
- parm.sfile->required = NO;
- parm.sfile->description =
+ parm.observation = G_define_standard_option(G_OPT_V_INPUT);
+ parm.observation->key = "observation";
+ parm.observation->required = NO;
+ parm.observation->description =
_("Name of the sampling locations vector points map");
- parm.sfile->guisection = _("Input_options");
-*/
+ parm.observation->guisection = _("Input_options");
+ parm.logfile = G_define_standard_option(G_OPT_F_OUTPUT);
+ parm.logfile->key = "logfile";
+ parm.logfile->required = NO;
+ parm.logfile->description =
+ _("Name of the sampling points output text file. For each observation vector point the time series of sediment transport is stored.");
+ parm.logfile->guisection = _("Output");
+
parm.tc = G_define_standard_option(G_OPT_R_OUTPUT);
parm.tc->key = "tc";
parm.tc->required = NO;
@@ -306,7 +311,6 @@
wdepth = parm.wdepth->answer;
dxin = parm.dxin->answer;
dyin = parm.dyin->answer;
- /* maskmap = parm.maskmap->answer; */
detin = parm.detin->answer;
tranin = parm.tranin->answer;
tauin = parm.tauin->answer;
@@ -348,9 +352,6 @@
G_message(_("Using metric conversion factor %f, step=%f"), conv,
step);
- /*
- * G_set_embedded_null_value_mode(1);
- */
if ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) &&
(erdep == NULL))
@@ -370,10 +371,6 @@
if (erdep != NULL || et != NULL)
er = G_alloc_fmatrix(my, mx);
- /* if (maskmap != NULL)
- bitmask = BM_create (cols, rows);
- IL_create_bitmask (¶ms, bitmask);
- */
seeds(rand1, rand2);
grad_check();
@@ -388,13 +385,6 @@
G_fatal_error(_("Cannot write raster maps"));
}
-/*
- if (fdwalkers != NULL)
- fclose(fdwalkers);
-
- if (sfile != NULL)
- G_sites_close(fw);
-*/
/* Exit with Success */
exit(EXIT_SUCCESS);
}
Modified: grass/trunk/raster/simwe/r.sim.water/main.c
===================================================================
--- grass/trunk/raster/simwe/r.sim.water/main.c 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/r.sim.water/main.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -91,16 +91,7 @@
#include <grass/waterglobs.h>
-char fncdsm[32];
-char filnam[10];
-/*struct BM *bitmask; */
-/*struct Cell_head cellhd; */
-struct GModule *module;
-struct Map_info Map;
-
-char msg[1024];
-
/****************************************/
/* MAIN */
@@ -112,6 +103,7 @@
double x_orig, y_orig;
static int rand1 = 12345;
static int rand2 = 67891;
+ struct GModule *module;
G_gisinit(argv[0]);
@@ -185,15 +177,20 @@
_("Name of flow controls raster map (permeability ratio 0-1)");
parm.traps->guisection = _("Input");
-/*
- parm.sfile = G_define_standard_option(G_OPT_V_INPUTy);
- parm.sfile->key = "vector";
- parm.sfile->required = NO;
- parm.sfile->description =
+ parm.observation = G_define_standard_option(G_OPT_V_INPUT);
+ parm.observation->key = "observation";
+ parm.observation->required = NO;
+ parm.observation->description =
_("Name of the sampling locations vector points map");
- parm.sfile->guisection = _("Input_options");
-*/
+ parm.observation->guisection = _("Input_options");
+ parm.logfile = G_define_standard_option(G_OPT_F_OUTPUT);
+ parm.logfile->key = "logfile";
+ parm.logfile->required = NO;
+ parm.logfile->description =
+ _("Name of the sampling points output text file. For each observation vector point the time series of water depth is stored.");
+ parm.logfile->guisection = _("Output");
+
parm.depth = G_define_standard_option(G_OPT_R_OUTPUT);
parm.depth->key = "depth";
parm.depth->required = NO;
@@ -333,11 +330,9 @@
disch = parm.disch->answer;
err = parm.err->answer;
outwalk = parm.outwalk->answer;
-/* sfile = parm.sfile->answer; */
sscanf(parm.niter->answer, "%d", ×ec);
sscanf(parm.outiter->answer, "%d", &iterout);
-/* sscanf(parm.density->answer, "%d", &ldemo); */
sscanf(parm.diffc->answer, "%lf", &frac);
sscanf(parm.hmax->answer, "%lf", &hhmax);
sscanf(parm.halpha->answer, "%lf", &halpha);
@@ -454,19 +449,6 @@
G_message(_("Using metric conversion factor %f, step=%f"), conv,
step);
- /*
- * G_set_embedded_null_value_mode(1);
- */
-
-/* replaced with condition that skips outwalk */
-/* if ((depth == NULL) && (disch == NULL) && (err == NULL) &&
- (outwalk == NULL))
- G_warning(_("You are not outputting any raster or vector points maps"));
- ret_val = input_data();
- if (ret_val != 1)
- G_fatal_error(_("Input failed"));
-*/
-
if ((depth == NULL) && (disch == NULL) && (err == NULL))
G_warning(_("You are not outputting any raster maps"));
ret_val = input_data();
@@ -482,10 +464,6 @@
gammas = G_alloc_matrix(my, mx);
dif = G_alloc_fmatrix(my, mx);
- /* if (maskmap != NULL)
- bitmask = BM_create (cols, rows);
- IL_create_bitmask (¶ms, bitmask);
- */
G_debug(2, "seeding randoms");
seeds(rand1, rand2);
grad_check();
@@ -497,13 +475,6 @@
G_fatal_error(_("Cannot write raster maps"));
}
-/*
- if (fdwalkers != NULL)
- fclose(fdwalkers);
-
- if (sfile != NULL)
- fclose(fw);
-*/
/* Exit with Success */
exit(EXIT_SUCCESS);
}
Modified: grass/trunk/raster/simwe/r.sim.water/r.sim.water.html
===================================================================
--- grass/trunk/raster/simwe/r.sim.water/r.sim.water.html 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/r.sim.water/r.sim.water.html 2011-03-15 10:16:37 UTC (rev 45670)
@@ -53,18 +53,23 @@
the <i>err</i> raster map (the resulting water depth is an average, and err is its RMSE).
The output vector points map <i>outwalk</i> can be used to analyze and visualize
spatial distribution of walkers at different simulation times (note that
-the resulting water depth is based on the density of these walkers). Number
-of the output walkers is controled by the <i>density</i> parameter, which controls
-how many walkers used in simulation should be written into the output.
+the resulting water depth is based on the density of these walkers).
+<!--Number of the output walkers is controled by the <i>density</i> parameter, which controls
+how many walkers used in simulation should be written into the output. -->
<!--(<font color="#ff0000"> toto treba upresnit/zmenit, lebo nwalk ide prec</font>). -->
Duration of simulation is controled by the <i>niter</i> parameter. The default value
is 10 minutes, reaching the steady-state may require much longer time,
depending on the time step, complexity of terrain, land cover and size of the area.
-Output water depth and discharge maps can be saved during simulation using
+Output walker, water depth and discharge maps can be saved during simulation using
the time series flag <i>-t</i> and <i>outiter</i> parameter
defining the time step in minutes for writing output files.
Files are saved with a suffix representing time since the start of simulation in seconds
(e.g. wdepth.500, wdepth.1000).
+Monitoring of water depth at specific points is supported. A vector map with observation points and
+a path to a logfile must be provided. For each point in the vector map which is located in
+the computational region the water depth is logged each time step in the logfile. The logfile is
+organized as a table. A single header identifies the category number of the logged vector points.
+In case of invalid water depth data the value -1 is used.
<P>
Overland flow is routed based on partial derivatives of elevation
Modified: grass/trunk/raster/simwe/r.sim.water/spearfish.sh
===================================================================
--- grass/trunk/raster/simwe/r.sim.water/spearfish.sh 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/r.sim.water/spearfish.sh 2011-03-15 10:16:37 UTC (rev 45670)
@@ -5,31 +5,41 @@
# andreas.philipp geo.uni-augsburg.de
# modified by JH, HM, MN, SG
-
dem=elevation.10m
output=simwe
g.region rast=${output}
g.region n=4920800 s=4917800 w=602500 e=606000
g.region -p
-# Manning n
-man05=0.05
-
# rainfall [mm/hr]
-rain01=2.5
+rain01=25
# infitration [mm/hr]
infil0=0.
echo "Preparing input maps..."
r.slope.aspect --o elevation=$dem dx=${output}_dx dy=${output}_dy
-r.mapcalc --o expr="${output}_rain =if(${dem},$rain01,null())"
-r.mapcalc --o expr="${output}_manin=if(${dem},$man05,null())"
-r.mapcalc --o expr="${output}_infil=if(${dem},$infil0,null())"
+r.mapcalc --o expr="${output}_rain =if(${dem}, $rain01, null())"
+r.mapcalc --o expr="${output}_manin=if(${dem}, soils * 0.005, null())"
+r.mapcalc --o expr="${output}_infil=if(${dem}, $infil0, null())"
+r.mapcalc --o expr="${output}_null=if(${dem}, float(0.0), null())"
-echo "r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1"
- r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1
+# Create the observation points
+v.random --o output=observation_points n=5 -z
+echo "r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy rain=${output}_rain man=${output}_manin infil=${output}_infil depth=${output}_depth disch=${output}_disch err=${output}_err outwalk=${output}_walker niter=5 outiter=1 observation=observation_points logfile=simwe.log"
+r.sim.water --o -t elevation=${dem} dx=${output}_dx dy=${output}_dy \
+ rain=${output}_rain man=${output}_manin infil=${output}_infil \
+ depth=${output}_depth disch=${output}_disch err=${output}_err \
+ outwalk=${output}_walker niter=5 outiter=1 observation=observation_points \
+ logfile=simwe.log
+
+echo "Plotting the observation points with R. Result in Rplots.pdf"
+R --vanilla -e "df = read.table(\"simwe.log\", header=1); par(mfrow=c(3,2)); plot(df\$STEP, df\$CAT0001, type=\"l\"); plot(df\$STEP, df\$CAT0002, type=\"l\"); plot(df\$STEP, df\$CAT0003, type=\"l\"); plot(df\$STEP, df\$CAT0004, type=\"l\");plot(df\$STEP, df\$CAT0005, type=\"l\")"
+
+echo "Export of manning, soil and gradients"
+r.out.vtk elevation=${dem} input=${output}_manin,soils vectormaps=${output}_dx,${output}_dy,${output}_null output=manning_soils_gradient.vtk null=0.0
+
echo "Build topology and exporting walker vector points for each time step"
for i in `g.mlist vect | grep ${output}` ; do
v.build $i
@@ -47,4 +57,5 @@
echo "Step throu the time steps and adjust the color tables"
# cleanup
-g.remove --q rast=${output}_dx,${output}_dy,${output}_rain,${output}_manin,${output}_infil
+g.remove --q rast=${output}_dx,${output}_dy,${output}_rain,${output}_manin,${output}_infil,${output}_null
+g.remove --q vect=observation_points
Modified: grass/trunk/raster/simwe/simlib/hydro.c
===================================================================
--- grass/trunk/raster/simwe/simlib/hydro.c 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/simlib/hydro.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -26,29 +26,16 @@
#include <grass/waterglobs.h>
-struct options parm;
-struct flags flag;
-
/*
* Soeren 8. Mar 2011 TODO:
* Put all these global variables into several meaningful structures and
* document use and purpose.
*
- * Example:
- * Put all file descriptors into a input_files struct and rename the variables:
- * input_files.elev
- * input_files.dx
- * input_files.dy
- * input_files.drain
- * ...
- *
*/
-FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps,
- *fdmanin, *fddepth, *fddisch, *fderr;
-FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc,
- *fdflux, *fderdep;
-FILE *fdsfile, *fw;
+struct options parm;
+struct flags flag;
+struct _points points;
char *elevin;
char *dxin;
@@ -57,7 +44,7 @@
char *infil;
char *traps;
char *manin;
-/* char *sfile; */
+/* char *observation; */
char *depth;
char *disch;
char *err;
@@ -161,6 +148,9 @@
nblock = mitfac + 1;
maxwa = maxwa / nblock;
}
+
+ /* Create the observation points */
+ create_observation_points();
G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
@@ -395,8 +385,8 @@
nstack++;
}
- }
- } /* lw loop */
+ } /* lw loop */
+ }
if (i == iter1 && ts == 1) {
/* call output for iteration output */
@@ -409,42 +399,34 @@
if (ii != 1)
G_fatal_error(_("Unable to write raster maps"));
}
-
- /* Soeren 8. Mar 2011 TODO:
- * This hould be replaced by vector functionality and sql database storage */
-/* ascii data site file output for gamma - hydrograph or sediment*/
-/* cchez incl. sqrt(sinsl) */
-/* sediment */
-/*defined area */
-/*
- if (sfile != NULL) {
-
- for (p = 0; p < npoints; p++) {
-
- l = (int)((points[p].east - mixx + stxm) / stepx) - mx -
- 1;
- k = (int)((points[p].north - miyy + stym) / stepy) - my -
- 1;
-
+
+ /* Write the water depth each time step at an observation point */
+ if(points.is_open)
+ {
+ double value = 0.0;
+ int p;
+ fprintf(points.output, "%.6d ", i);
+ /* Write for each point */
+ for(p = 0; p < points.npoints; p++)
+ {
+ l = (int)((points.x[p] - mixx + stxm) / stepx) - mx - 1;
+ k = (int)((points.y[p] - miyy + stym) / stepy) - my - 1;
+
if (zz[k][l] != UNDEF) {
- if (wdepth == NULL) {
- points[p].z1 = step * gama[k][l] * cchez[k][l];
- }
+ if (wdepth == NULL)
+ value = step * gama[k][l] * cchez[k][l];
else
- points[p].z1 = gama[k][l] * slope[k][l];
+ value = gama[k][l] * slope[k][l];
- G_debug(2, " k,l,z1 %d %d %f", k, l, points[p].z1);
-
- fprintf(fw, "%f %f %f\n", points[p].east / conv,
- points[p].north / conv, points[p].z1);
- }
-
- }
-
- }
-*/
-
+ fprintf(points.output, "%2.4f ", value);
+ } else {
+ /* Point is invalid, so a negative value is written */
+ fprintf(points.output, "%2.4f ", -1.0);
+ }
+ }
+ fprintf(points.output, "\n");
+ }
} /* miter */
L_800:
@@ -475,5 +457,11 @@
erod(gama);
}
/* ........ end of iblock loop */
+
+ /* Close the observation logfile */
+ if(points.is_open)
+ fclose(points.output);
+
+ points.is_open = 0;
}
Modified: grass/trunk/raster/simwe/simlib/input.c
===================================================================
--- grass/trunk/raster/simwe/simlib/input.c 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/simlib/input.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -4,14 +4,21 @@
#include <stdlib.h>
#include <math.h>
#include <grass/gis.h>
-/* #include <grass/site.h> */
-#include <grass/bitmap.h>
#include <grass/glocale.h>
#include <grass/linkm.h>
#include <grass/gmath.h>
#include <grass/waterglobs.h>
+/* Local prototypes for raster map reading and array allocation */
+static float ** read_float_raster_map(int rows, int cols, char *name, float unitconv);
+static double ** read_double_raster_map(int rows, int cols, char *name, double unitconv);
+static float ** create_float_matrix(int rows, int cols, float fill_value);
+static double ** create_double_matrix(int rows, int cols, double fill_value);
+static void copy_matrix_undef_double_to_float_values(int rows, int cols, double **source, float **target);
+static void copy_matrix_undef_float_values(int rows, int cols, float **source, float **target);
+
+
/* ************************************************************** */
/* GRASS input procedures, allocations */
/* *************************************************************** */
@@ -20,406 +27,104 @@
* \brief allocate memory, read input rasters, assign UNDEF to NODATA
*
* \return int
- * sites related input/output commented out - needs update to vect, HM nov 2008
*/
-
+/* ************************************************************************* */
+/* Read all input maps and input values into memory ************************ */
int input_data(void)
{
-
- FCELL *elevin_cell, *traps_cell, *manin_cell;
- FCELL *detin_cell, *trainin_cell, *tauin_cell;
- DCELL *dxin_cell, *dyin_cell, *rain_cell, *infil_cell, *wdepth_cell;
- int elevin_fd, dxin_fd, dyin_fd, rain_fd, infil_fd, traps_fd, manin_fd, row, row_rev;
- int detin_fd, trainin_fd, tauin_fd, wdepth_fd;
- int j;
-/* int nn, cc, ii, dd; */
+ int rows = my, cols = mx; /* my and mx are global variables */
double unitconv = 0.0000002; /* mm/hr to m/s */
- const char *mapset;
-/* output water depth and discharge at outlet points given in site file*/
-/* Site *site;
+ int if_rain = 0;
- npoints = 0;
- npoints_alloc = 0;
+ G_debug(1, "Running MAR 2011 version, started modifications on 20080211");
+
+ /* Elevation and gradients are mandatory */
+ zz = read_float_raster_map(rows, cols, elevin, 1.0);
+ v1 = read_double_raster_map(rows, cols, dxin, 1.0);
+ v2 = read_double_raster_map(rows, cols, dyin, 1.0);
- if (sfile != NULL) {
- fw = fopen("simwe_data.txt", "w");
+ /* Update elevation map */
+ copy_matrix_undef_double_to_float_values(rows, cols, v1, zz);
+ copy_matrix_undef_double_to_float_values(rows, cols, v2, zz);
- mapset = G_find_sites(sfile, "");
- if (mapset == NULL)
- G_fatal_error(_("File [%s] not found"), sfile);
+ /* Manning surface roughnes: read map or use a single value */
+ if(manin != NULL) {
+ cchez = read_float_raster_map(rows, cols, manin, 1.0);
+ } else if(manin_val >= 0.0) { /* If no value set its set to -999.99 */
+ cchez = create_float_matrix(rows, cols, manin_val * unitconv);
+ }else{
+ G_fatal_error(_("Raster map <%s> not found, and manin_val undefined, choose one to be allowed to process"), manin);
+ }
+
+ /* Rain: read rain map or use a single value for all cells */
+ if (rain != NULL) {
+ si = read_double_raster_map(rows, cols, rain, unitconv);
+ if_rain = 1;
+ } else if(rain_val >= 0.0) { /* If no value set its set to -999.99 */
+ si = create_double_matrix(rows, cols, rain_val * unitconv);
+ if_rain = 1;
+ } else{
+ si = create_double_matrix(rows, cols, (double)UNDEF);
+ if_rain = 0;
+ }
- if ((fdsfile = G_fopen_sites_old(sfile, mapset)) == NULL)
- G_fatal_error(_("Unable to open file [%s]"), sfile);
+ /* Update elevation map */
+ copy_matrix_undef_double_to_float_values(rows, cols, si, zz);
- if (G_site_describe(fdsfile, &nn, &cc, &ii, &dd) != 0)
- G_fatal_error(_("Failed to guess file format"));
+ /* Load infiltration and traps if rain is present */
+ if(if_rain == 1) {
+ /* Infiltration: read map or use a single value */
+ if (infil != NULL) {
+ inf = read_double_raster_map(rows, cols, infil, unitconv);
+ } else if(infil_val >= 0.0) { /* If no value set its set to -999.99 */
+ inf = create_double_matrix(rows, cols, infil_val * unitconv);
+ } else{
+ inf = create_double_matrix(rows, cols, (double)UNDEF);
+ }
- site = G_site_new_struct(cc, nn, ii, dd);
- G_message(_("Reading sites map (%s) ..."), sfile);
+ /* Traps */
+ if (traps != NULL)
+ trap = read_float_raster_map(rows, cols, traps, 1.0);
+ else
+ trap = create_float_matrix(rows, cols, (double)UNDEF);
+ }
- if (dd==0)
- {
- fprintf(stderr,"\n");
- G_warning("I'm finding records that do not have
- a floating point attributes (fields prefixed with '%').");
- }
+ if (detin != NULL) {
+ dc = read_float_raster_map(rows, cols, detin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, dc, zz);
+ }
- while (G_site_get(fdsfile, site) >= 0) {
- if (npoints_alloc <= npoints) {
- npoints_alloc += 128;
- points =
- (struct Point *)G_realloc(points,
- npoints_alloc *
- sizeof(struct Point));
- }
- points[npoints].east = site->east * conv;
- points[npoints].north = site->north * conv;
- points[npoints].z1 = 0.;
- if ((points[npoints].east / conv <= cellhd.east &&
- points[npoints].east / conv >= cellhd.west) &&
- (points[npoints].north / conv <= cellhd.north &&
- points[npoints].north / conv >= cellhd.south))
- npoints++;
- }
- G_sites_close(fdsfile);
+ if (tranin != NULL) {
+ ct = read_float_raster_map(rows, cols, tranin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, ct, zz);
}
-*/
- /* Allocate raster buffers */
- elevin_cell = Rast_allocate_f_buf();
- dxin_cell = Rast_allocate_d_buf();
- dyin_cell = Rast_allocate_d_buf();
+ if (tauin != NULL) {
+ tau = read_float_raster_map(rows, cols, tauin, 1.0);
+ copy_matrix_undef_float_values(rows, cols, tau, zz);
+ }
- if(manin != NULL)
- manin_cell = Rast_allocate_f_buf();
-
- if (rain != NULL)
- rain_cell = Rast_allocate_d_buf();
-
- if (infil != NULL)
- infil_cell = Rast_allocate_d_buf();
-
- if (traps != NULL)
- traps_cell = Rast_allocate_f_buf();
-
- if (detin != NULL)
- detin_cell = Rast_allocate_f_buf();
-
- if (tranin != NULL)
- trainin_cell = Rast_allocate_f_buf();
-
- if (tauin != NULL)
- tauin_cell = Rast_allocate_f_buf();
-
- if (wdepth != NULL)
- wdepth_cell = Rast_allocate_d_buf();
-
- /* Allocate some double dimension arrays for each input */
- zz = G_alloc_fmatrix(my, mx);
- v1 = G_alloc_matrix(my, mx);
- v2 = G_alloc_matrix(my, mx);
- cchez = G_alloc_fmatrix(my, mx);
+ if (wdepth != NULL) {
+ gama = read_double_raster_map(rows, cols, wdepth, 1.0);
+ copy_matrix_undef_double_to_float_values(rows, cols, gama, zz);
+ }
- if (rain != NULL || rain_val >= 0.0)
- si = G_alloc_matrix(my, mx);
+ /* Array for gradient checking */
+ slope = create_double_matrix(rows, cols, 0.0);
+
+ /* Create the observation points and open the logfile */
+ create_observation_points();
- if (infil != NULL || infil_val >= 0.0)
- inf = G_alloc_matrix(my, mx);
-
- if (traps != NULL)
- trap = G_alloc_fmatrix(my, mx);
-
- if (detin != NULL)
- dc = G_alloc_fmatrix(my, mx);
-
- if (tranin != NULL)
- ct = G_alloc_fmatrix(my, mx);
-
- if (tauin != NULL)
- tau = G_alloc_fmatrix(my, mx);
-
- if (wdepth != NULL)
- gama = G_alloc_matrix(my, mx);
-
- G_debug(3, "Running MAR 2011 version, started modifications on 20080211");
-
- /* Check if data available in mapsets
- * if found, then open the files */
- elevin_fd = Rast_open_old(elevin, "");
-
- /* TO REPLACE BY INTERNAL PROCESSING of dx, dy from Elevin */
- if ((mapset = G_find_raster(dxin, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"), dxin);
-
- dxin_fd = Rast_open_old(dxin, "");
-
- dyin_fd = Rast_open_old(dyin, "");
- /* END OF REPLACEMENT */
-
- /* Rendered Mannings n input map optional to run! */
- /* Careful! (Yann, 20080212) */
- if (manin)
- manin_fd = Rast_open_old(manin, "");
-
- /* Rendered Rainfall input map optional to run! */
- /* Careful! (Yann, 20080212) */
- if (rain)
- rain_fd = Rast_open_old(rain, "");
-
- if (infil)
- infil_fd = Rast_open_old(infil, "");
-
- if (traps)
- traps_fd = Rast_open_old(traps, "");
-
- if (detin)
- detin_fd = Rast_open_old(detin, "");
-
- if (tranin)
- trainin_fd = Rast_open_old(tranin, "");
-
- if (tauin)
- tauin_fd = Rast_open_old(tauin, "");
-
- if (wdepth)
- wdepth_fd = Rast_open_old(wdepth, "");
-
- for (row = 0; row < my; row++) {
- Rast_get_f_row(elevin_fd, elevin_cell, row);
- Rast_get_d_row(dxin_fd, dxin_cell, row);
- Rast_get_d_row(dyin_fd, dyin_cell, row);
-
- if (manin)
- Rast_get_f_row(manin_fd, manin_cell, row);
-
- if (rain)
- Rast_get_d_row(rain_fd, rain_cell, row);
-
- if (infil)
- Rast_get_d_row(infil_fd, infil_cell, row);
-
- if (traps)
- Rast_get_f_row(traps_fd, traps_cell, row);
-
- if (detin)
- Rast_get_f_row(detin_fd, detin_cell, row);
-
- if (tranin)
- Rast_get_f_row(trainin_fd, trainin_cell, row);
-
- if (tauin)
- Rast_get_f_row(tauin_fd, tauin_cell, row);
-
- if (wdepth)
- Rast_get_d_row(wdepth_fd, wdepth_cell, row);
-
- for (j = 0; j < mx; j++) {
- row_rev = my - row - 1;
- /*if elevation data exists store in zz[][] */
- if (!Rast_is_f_null_value(elevin_cell + j))
- zz[row_rev][j] = (float)(conv * elevin_cell[j]);
- else
- zz[row_rev][j] = UNDEF;
-
- if (!Rast_is_d_null_value(dxin_cell + j))
- v1[row_rev][j] = (double)dxin_cell[j];
- else
- v1[row_rev][j] = UNDEF;
-
- if (!Rast_is_d_null_value(dyin_cell + j))
- v2[row_rev][j] = (double)dyin_cell[j];
- else
- v2[row_rev][j] = UNDEF;
-
- /* undef all area if something's missing */
- if (v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF)
- zz[row_rev][j] = UNDEF;
-
- /* should be ?
- * if(v1[row_rev][j] == UNDEF || v2[row_rev][j] == UNDEF ||
- * zz[row_rev][j] == UNDEF) {
- * v1[row_rev][j] == UNDEF;
- * v2[row_rev][j] == UNDEF;
- * zz[row_rev][j] == UNDEF;
- * }
- *//*printout warning? */
-
- /* If Rain Exists, then load data */
- if (rain) {
- if (!Rast_is_d_null_value(rain_cell + j))
- si[row_rev][j] = ((double)rain_cell[j]) * unitconv;
- /*conv mm/hr to m/s */
- /*printf("\n INPUTrain, convert %f %f",si[row_rev][j],unitconv); */
-
- else {
- si[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
-
- /* Load infiltration map too if it exists */
- if (infil) {
- if (!Rast_is_d_null_value(infil_cell + j))
- inf[row_rev][j] = (double)infil_cell[j] * unitconv;
- /*conv mm/hr to m/s */
- /*printf("\nINPUT infilt,convert %f %f",inf[row_rev][j],unitconv); */
- else {
- inf[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- else { /* Added by Yann 20080216 */
- /* If infil==NULL, then use infilval */
- if (infil_val >= 0.0) {
- inf[row_rev][j] = infil_val * unitconv; /*conv mm/hr to m/s */
- /* printf("infil_val = %f \n",inf[row_rev][j]); */
- }
- else {
- inf[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
-
- if (traps) {
- if (!Rast_is_f_null_value(traps_cell + j))
- trap[row_rev][j] = (float)traps_cell[j]; /* no conv, unitless */
- else {
- trap[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- }
- else { /* Added by Yann 20080213 */
- /* If rain==NULL, then use rainval */
- if (rain_val >= 0.0) {
- si[row_rev][j] = rain_val * unitconv; /* conv mm/hr to m/s */
- /*printf("\n INPUTrainval, convert %f %f",si[row_rev][j],unitconv); */
- }
- else {
- si[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
-
- if (infil) {
- if (!Rast_is_d_null_value(infil_cell + j))
- inf[row_rev][j] = (double)infil_cell[j] * unitconv; /*conv mm/hr to m/s */
- /*printf("\nINPUT infilt,convert %f %f",inf[row_rev][j],unitconv); */
- else {
- inf[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- else { /* Added by Yann 20080216 */
- /* If infil==NULL, then use infilval */
- if (infil_val >= 0.0) {
- inf[row_rev][j] = infil_val * unitconv; /*conv mm/hr to m/s */
- /*printf("infil_val = %f \n",inf[row_rev][j]); */
- }
- else {
- inf[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
-
- if (traps) {
- if (!Rast_is_f_null_value(traps_cell + j))
- trap[row_rev][j] = (float)traps_cell[j]; /* no conv, unitless */
- else {
- trap[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- } /* End of added by Yann 20080213 */
- if (manin) {
- if (!Rast_is_f_null_value(manin_cell + j)) {
- cchez[row_rev][j] = (float)manin_cell[j]; /* units in manual */
- }
- else {
- cchez[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- else if (manin_val >= 0.0) { /* Added by Yann 20080213 */
- cchez[row_rev][j] = (float)manin_val;
- }
- else {
- G_fatal_error(_("Raster map <%s> not found, and manin_val undefined, choose one to be allowed to process"),
- manin);
- }
- if (detin) {
- if (!Rast_is_f_null_value(detin_cell + j))
- dc[row_rev][j] = (float)detin_cell[j]; /*units in manual */
- else {
- dc[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
-
- if (tranin) {
- if (!Rast_is_f_null_value(trainin_cell + j))
- ct[row_rev][j] = (float)trainin_cell[j]; /*units in manual */
- else {
- ct[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
-
- if (tauin) {
- if (!Rast_is_f_null_value(tauin_cell + j))
- tau[row_rev][j] = (float)tauin_cell[j]; /*units in manual */
- else {
- tau[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
-
- if (wdepth) {
- if (!Rast_is_d_null_value(wdepth_cell + j))
- gama[row_rev][j] = (double)wdepth_cell[j]; /*units in manual */
- else {
- gama[row_rev][j] = UNDEF;
- zz[row_rev][j] = UNDEF;
- }
- }
- }
- }
- Rast_close(elevin_fd);
- Rast_close(dxin_fd);
- Rast_close(dyin_fd);
-
- if (rain)
- Rast_close(rain_fd);
-
- if (infil)
- Rast_close(infil_fd);
-
- if (traps)
- Rast_close(traps_fd);
- /* Maybe a conditional to manin!=NULL here ! */
- Rast_close(manin_fd);
-
- /****************/
-
- if (detin)
- Rast_close(detin_fd);
-
- if (tranin)
- Rast_close(trainin_fd);
-
- if (tauin)
- Rast_close(tauin_fd);
-
- if (wdepth)
- Rast_close(wdepth_fd);
-
- return 1;
+ return 1;
}
+/* ************************************************************************* */
/* data preparations, sigma, shear, etc. */
int grad_check(void)
{
- int k, l, i, j;
+ int k, l;
double zx, zy, zd2, zd4, sinsl;
double cc, cmul2;
double sheer;
@@ -445,19 +150,6 @@
infsum = 0.;
cmul2 = rhow * gacc;
- /* mandatory alloc. - should be moved to main.c */
- slope = (double **)G_malloc(sizeof(double *) * (my));
-
- for (l = 0; l < my; l++)
- slope[l] = (double *)G_malloc(sizeof(double) * (mx));
-
- for (j = 0; j < my; j++) {
- for (i = 0; i < mx; i++)
- slope[j][i] = 0.;
- }
-
- /*** */
-
for (k = 0; k < my; k++) {
for (l = 0; l < mx; l++) {
if (zz[k][l] != UNDEF) {
@@ -633,7 +325,9 @@
if (sigma[k][l] != 0.)
/* rate of weight loss - w=w*sigma ,
* vaha prechadzky po n-krokoch je sigma^n */
- /* not clear what's here :-\ */
+
+ /*!!!!! not clear what's here :-\ !!!!!*/
+
sigma[k][l] =
exp(-sigma[k][l] * deltap * slope[k][l]);
/* if(sigma[k][l]<0.5) warning, napocitaj,
@@ -645,62 +339,154 @@
return 1;
}
+/* ************************************************************************* */
-double amax1(double arg1, double arg2)
+void copy_matrix_undef_double_to_float_values(int rows, int cols, double **source, float **target)
{
- double res;
+ int col = 0, row = 0;
- if (arg1 >= arg2) {
- res = arg1;
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ if(source[row][col] == UNDEF)
+ target[row][col] = UNDEF;
+ }
}
- else {
- res = arg2;
- }
+}
- return res;
+/* ************************************************************************* */
+
+void copy_matrix_undef_float_values(int rows, int cols, float **source, float **target)
+{
+ int col = 0, row = 0;
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ if(source[row][col] == UNDEF)
+ target[row][col] = UNDEF;
+ }
+ }
}
+/* ************************************************************************* */
-double amin1(double arg1, double arg2)
+float ** create_float_matrix(int rows, int cols, float fill_value)
{
- double res;
+ int col = 0, row = 0;
+ float **matrix = NULL;
- if (arg1 <= arg2) {
- res = arg1;
+ /* Allocate the float marix */
+ matrix = G_alloc_fmatrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ matrix[row][col] = fill_value;
+ }
}
- else {
- res = arg2;
- }
- return res;
+ return matrix;
}
+/* ************************************************************************* */
-int min(int arg1, int arg2)
+double ** create_double_matrix(int rows, int cols, double fill_value)
{
- int res;
+ int col = 0, row = 0;
+ double **matrix = NULL;
- if (arg1 <= arg2) {
- res = arg1;
+ /* Allocate the float marix */
+ matrix = G_alloc_matrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ for(col = 0; col < cols; col++) {
+ matrix[row][col] = fill_value;
+ }
}
- else {
- res = arg2;
- }
- return res;
+ return matrix;
}
+/* ************************************************************************* */
-int max(int arg1, int arg2)
+float ** read_float_raster_map(int rows, int cols, char *name, float unitconv)
{
- int res;
+ FCELL *row_buff = NULL;
+ int fd;
+ int col = 0, row = 0, row_rev = 0;
+ float **matrix = NULL;
- if (arg1 >= arg2) {
- res = arg1;
+ G_message("Reading float map %s into memory", name);
+
+ /* Open raster map */
+ fd = Rast_open_old(name, "");
+
+ /* Allocate the row buffer */
+ row_buff = Rast_allocate_f_buf();
+
+ /* Allocate the float marix */
+ matrix = G_alloc_fmatrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ Rast_get_f_row(fd, row_buff, row);
+
+ for(col = 0; col < cols; col++) {
+ /* we fill the arrays from south to north */
+ row_rev = rows - row - 1;
+ /* Check for null values */
+ if (!Rast_is_f_null_value(row_buff + col))
+ matrix[row_rev][col] = (float)(unitconv * row_buff[col]);
+ else
+ matrix[row_rev][col] = UNDEF;
+ }
}
- else {
- res = arg2;
- }
- return res;
+ /* Free the row buffer */
+ if(row_buff)
+ G_free(row_buff);
+
+ Rast_close(fd);
+
+ return matrix;
}
+
+/* ************************************************************************* */
+
+double ** read_double_raster_map(int rows, int cols, char *name, double unitconv)
+{
+ DCELL *row_buff = NULL;
+ int fd;
+ int col = 0, row = 0, row_rev;
+ double **matrix = NULL;
+
+ G_message("Reading double map %s into memory", name);
+
+ /* Open raster map */
+ fd = Rast_open_old(name, "");
+
+ /* Allocate the row buffer */
+ row_buff = Rast_allocate_d_buf();
+
+ /* Allocate the double marix */
+ matrix = G_alloc_matrix(rows, cols);
+
+ for(row = 0; row < rows; row++) {
+ Rast_get_d_row(fd, row_buff, row);
+
+ for(col = 0; col < cols; col++) {
+ /* we fill the arrays from south to north */
+ row_rev = rows - row - 1;
+ /* Check for null values */
+ if (!Rast_is_d_null_value(row_buff + col))
+ matrix[row_rev][col] = (double)(unitconv * row_buff[col]);
+ else
+ matrix[row_rev][col] = UNDEF;
+ }
+ }
+
+ /* Free the row buffer */
+ if(row_buff)
+ G_free(row_buff);
+
+ Rast_close(fd);
+
+ return matrix;
+}
\ No newline at end of file
Added: grass/trunk/raster/simwe/simlib/observation_points.c
===================================================================
--- grass/trunk/raster/simwe/simlib/observation_points.c (rev 0)
+++ grass/trunk/raster/simwe/simlib/observation_points.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -0,0 +1,133 @@
+/* obervation_points.c (simlib), 14.mar.2011, SG */
+
+#include <grass/glocale.h>
+#include <grass/vector.h>
+#include <grass/waterglobs.h>
+
+
+static void init_points(struct _points *, int);
+static void realloc_points(struct _points *, int);
+static void insert_next_point(struct _points *p, double x, double y, int cat);
+
+/* ************************************************************************* */
+
+void create_observation_points()
+{
+ int if_log = 0;
+ struct Map_info Map;
+ struct line_pnts *pts;
+ struct line_cats *cts;
+ double x, y;
+ int type, cat, i;
+
+ if(parm.observation->answer != NULL)
+ if_log += 1;
+
+ if(parm.logfile->answer != NULL)
+ if_log += 1;
+
+ /* Nothing to do */
+ if(if_log == 0)
+ return;
+
+ if(if_log == 1)
+ G_fatal_error("Observation vector map and logfile must be provided");
+
+ Vect_set_open_level(1);
+
+ if (Vect_open_old(&Map, parm.observation->answer, "") < 0)
+ G_fatal_error(_("Unable to open vector map <%s>"), parm.observation->answer);
+
+ Vect_rewind(&Map);
+
+ pts = Vect_new_line_struct();
+ cts = Vect_new_cats_struct();
+
+ /* Initialize point structure */
+ init_points(&points, 128);
+
+ /* Read all vector points */
+ while(1) {
+ type = Vect_read_next_line(&Map, pts, cts);
+
+ if(type == -2) {
+ break;
+ }
+
+ if(type == -1) {
+ Vect_close(&Map);
+ G_fatal_error(_("Unable to read points from map %s"), parm.observation->answer);
+ }
+
+ if(type == GV_POINT) {
+ x = pts->x[0];
+ y = pts->y[0];
+ cat = cts->cat[0];
+
+ /* Check region bounds befor inserting point */
+ if(x <= cellhd.east && x >= cellhd.west && y <= cellhd.north && y >= cellhd.south) {
+ insert_next_point(&points, x, y, cat);
+ }
+ }
+ }
+
+ Vect_close(&Map);
+
+ /* Open the logfile */
+ points.output = fopen(parm.logfile->answer, "w");
+
+ if(points.output == NULL)
+ G_fatal_error(_("Unable to open observation logfile %s for writing"), parm.logfile->answer);
+
+ points.is_open = 1;
+
+ fprintf(points.output, "STEP ");
+
+ /* Write the vector cats as header information in the logfile */
+ for(i = 0; i < points.npoints; i++)
+ {
+ fprintf(points.output, "CAT%.4d ", points.cats[i]);
+ }
+ fprintf(points.output, "\n");
+
+
+ return;
+}
+
+/* ************************************************************************* */
+
+void init_points(struct _points *p, int size)
+{
+ p->x = (double*)G_calloc(size, sizeof(double));
+ p->y = (double*)G_calloc(size, sizeof(double));
+ p->cats = (int*)G_calloc(size, sizeof(int));
+ p->npoints = 0;
+ p->npoints_alloc = size;
+ p->output = NULL;
+ p->is_open = 0;
+}
+
+/* ************************************************************************* */
+
+void realloc_points(struct _points *p, int add_size)
+{
+ p->x = (double*)G_realloc(p->x, (p->npoints_alloc + add_size) * sizeof(double));
+ p->y = (double*)G_realloc(p->y, (p->npoints_alloc + add_size) * sizeof(double));
+ p->cats = (int*)G_realloc(p->cats, (p->npoints_alloc + add_size) * sizeof(int));
+ p->npoints_alloc += add_size;
+}
+
+/* ************************************************************************* */
+
+void insert_next_point(struct _points *p, double x, double y, int cat)
+{
+ if(p->npoints == p->npoints_alloc)
+ realloc_points(p, 128);
+
+ G_debug(3, "Insert point %g %g %i id %i\n", x, y, cat, p->npoints);
+
+ p->x[p->npoints] = x;
+ p->y[p->npoints] = y;
+ p->cats[p->npoints] = cat;
+ p->npoints++;
+}
Modified: grass/trunk/raster/simwe/simlib/output.c
===================================================================
--- grass/trunk/raster/simwe/simlib/output.c 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/simlib/output.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -15,9 +15,9 @@
static void output_walker_as_vector(int tt, int ndigit);
-
-/* This function was added by Soeren 8. Mar 2011
- * It replaces the site walker output implementation */
+/* This function was added by Soeren 8. Mar 2011 */
+/* It replaces the site walker output implementation */
+/* Only the 3d coordinates of the walker are stored. */
void output_walker_as_vector(int tt, int ndigit)
{
char buf[256];
@@ -32,7 +32,7 @@
/* In case of time series we extent the output name with the time value */
if (ts == 1) {
- sprintf(buf, "%s_%.*d", outwalk, ndigit, tt);
+ G_snprintf(buf, 256, "%s_%.*d", outwalk, ndigit, tt);
outwalk_time = G_store(buf);
Vect_open_new(&Out, outwalk_time, WITH_Z);
G_message("Writing %i walker into vector file %s", nstack, outwalk_time);
@@ -104,10 +104,18 @@
output_walker_as_vector(tt, ndigit);
Rast_set_window(&cellhd);
+
+ if (my != Rast_window_rows())
+ G_fatal_error("OOPS: rows changed from %d to %d\n", mx,
+ Rast_window_rows());
+ if (mx != Rast_window_cols())
+ G_fatal_error("OOPS: cols changed from %d to %d\n", my,
+ Rast_window_cols());
+
if (depth) {
depth_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", depth, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", depth, ndigit, tt);
depth0 = G_store(buf);
depth_fd = Rast_open_fp_new(depth0);
}
@@ -118,7 +126,7 @@
if (disch) {
disch_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", disch, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", disch, ndigit, tt);
disch0 = G_store(buf);
disch_fd = Rast_open_fp_new(disch0);
}
@@ -129,7 +137,7 @@
if (err) {
err_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", err, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", err, ndigit, tt);
err0 = G_store(buf);
err_fd = Rast_open_fp_new(err0);
}
@@ -140,7 +148,7 @@
if (conc) {
conc_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", conc, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", conc, ndigit, tt);
conc0 = G_store(buf);
conc_fd = Rast_open_fp_new(conc0);
}
@@ -151,7 +159,7 @@
if (flux) {
flux_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", flux, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", flux, ndigit, tt);
flux0 = G_store(buf);
flux_fd = Rast_open_fp_new(flux0);
}
@@ -162,7 +170,7 @@
if (erdep) {
erdep_cell = Rast_allocate_f_buf();
if (ts == 1) {
- sprintf(buf, "%s.%.*d", erdep, ndigit, tt);
+ G_snprintf(buf, 256,"%s.%.*d", erdep, ndigit, tt);
erdep0 = G_store(buf);
erdep_fd = Rast_open_fp_new(erdep0);
}
@@ -170,13 +178,6 @@
erdep_fd = Rast_open_fp_new(erdep);
}
- if (my != Rast_window_rows())
- G_fatal_error("OOPS: rows changed from %d to %d\n", mx,
- Rast_window_rows());
- if (mx != Rast_window_cols())
- G_fatal_error("OOPS: cols changed from %d to %d\n", my,
- Rast_window_cols());
-
for (iarc = 0; iarc < my; iarc++) {
i = my - iarc - 1;
if (depth) {
Added: grass/trunk/raster/simwe/simlib/utils.c
===================================================================
--- grass/trunk/raster/simwe/simlib/utils.c (rev 0)
+++ grass/trunk/raster/simwe/simlib/utils.c 2011-03-15 10:16:37 UTC (rev 45670)
@@ -0,0 +1,60 @@
+/* output.c (simlib), 20.nov.2002, JH */
+
+#include <grass/waterglobs.h>
+
+double amax1(double arg1, double arg2)
+{
+ double res;
+
+ if (arg1 >= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+double amin1(double arg1, double arg2)
+{
+ double res;
+
+ if (arg1 <= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+int min(int arg1, int arg2)
+{
+ int res;
+
+ if (arg1 <= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
+int max(int arg1, int arg2)
+{
+ int res;
+
+ if (arg1 >= arg2) {
+ res = arg1;
+ }
+ else {
+ res = arg2;
+ }
+
+ return res;
+}
+
Modified: grass/trunk/raster/simwe/simlib/waterglobs.h
===================================================================
--- grass/trunk/raster/simwe/simlib/waterglobs.h 2011-03-15 05:14:23 UTC (rev 45669)
+++ grass/trunk/raster/simwe/simlib/waterglobs.h 2011-03-15 10:16:37 UTC (rev 45670)
@@ -7,16 +7,6 @@
#include <grass/raster.h>
-/*
-extern FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps,
- *fdmanin, *fddepth, *fddisch, *fderr, *fdoutwalk, *fdwalkers;
-*/
-extern FILE *fdelevin, *fddxin, *fddyin, *fdrain, *fdinfil, *fdtraps,
- *fdmanin, *fddepth, *fddisch, *fderr;
-extern FILE *fdwdepth, *fddetin, *fdtranin, *fdtauin, *fdtc, *fdet, *fdconc,
- *fdflux, *fderdep;
-extern FILE *fdsfile, *fw;
-
extern char *elevin;
extern char *dxin;
extern char *dyin;
@@ -24,7 +14,7 @@
extern char *infil;
extern char *traps;
extern char *manin;
-/* extern char *sfile; */
+/* extern char *observation; */
extern char *depth;
extern char *disch;
extern char *err;
@@ -50,10 +40,10 @@
struct options
{
struct Option *elevin, *dxin, *dyin, *rain, *infil, *traps, *manin,
- *sfile, *depth, *disch, *err, *outwalk, *nwalk, *niter, *outiter,
+ *observation, *depth, *disch, *err, *outwalk, *nwalk, *niter, *outiter,
*density, *diffc, *hmax, *halpha, *hbeta, *wdepth, *detin, *tranin,
*tauin, *tc, *et, *conc, *flux, *erdep, *rainval, *maninval,
- *infilval;
+ *infilval, *logfile;
};
extern struct options parm;
@@ -75,18 +65,18 @@
extern struct Cell_head cellhd;
-struct Point
+struct _points
{
- double north, east;
- double z1;
+ double *x; /* x coor for each point */
+ double *y; /* y coor for each point*/
+ int *cats; /* Category for each point */
+ int npoints; /* Number of observation points */
+ int npoints_alloc; /* Number of allocated points */
+ FILE *output; /* Output file descriptor */
+ int is_open; /* Set to 1 if open, 0 if closed */
};
-/*
-extern struct Point *points;
-extern int npoints;
-extern int npoints_alloc;
-*/
-
+extern struct _points points;
extern int input_data(void);
extern int seeds(long int, long int);
extern int seedg(long int, long int);
@@ -101,6 +91,7 @@
extern double amin1(double, double);
extern int min(int, int);
extern int max(int, int);
+extern void create_observation_points();
extern double xmin, ymin, xmax, ymax;
extern double mayy, miyy, maxx, mixx;
More information about the grass-commit
mailing list