|
|
@@ -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;
|
|
|
-
|
|
|
- npoints = 0;
|
|
|
- npoints_alloc = 0;
|
|
|
-
|
|
|
- if (sfile != NULL) {
|
|
|
- fw = fopen("simwe_data.txt", "w");
|
|
|
-
|
|
|
- mapset = G_find_sites(sfile, "");
|
|
|
- if (mapset == NULL)
|
|
|
- G_fatal_error(_("File [%s] not found"), sfile);
|
|
|
-
|
|
|
- if ((fdsfile = G_fopen_sites_old(sfile, mapset)) == NULL)
|
|
|
- G_fatal_error(_("Unable to open file [%s]"), sfile);
|
|
|
-
|
|
|
- if (G_site_describe(fdsfile, &nn, &cc, &ii, &dd) != 0)
|
|
|
- G_fatal_error(_("Failed to guess file format"));
|
|
|
-
|
|
|
- site = G_site_new_struct(cc, nn, ii, dd);
|
|
|
- G_message(_("Reading sites map (%s) ..."), sfile);
|
|
|
-
|
|
|
- if (dd==0)
|
|
|
- {
|
|
|
- fprintf(stderr,"\n");
|
|
|
- G_warning("I'm finding records that do not have
|
|
|
- a floating point attributes (fields prefixed with '%').");
|
|
|
- }
|
|
|
-
|
|
|
- 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);
|
|
|
- }
|
|
|
-*/
|
|
|
-
|
|
|
- /* Allocate raster buffers */
|
|
|
- elevin_cell = Rast_allocate_f_buf();
|
|
|
- dxin_cell = Rast_allocate_d_buf();
|
|
|
- dyin_cell = Rast_allocate_d_buf();
|
|
|
-
|
|
|
- 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();
|
|
|
+ int if_rain = 0;
|
|
|
|
|
|
- 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);
|
|
|
+ G_debug(1, "Running MAR 2011 version, started modifications on 20080211");
|
|
|
|
|
|
- if (rain != NULL || rain_val >= 0.0)
|
|
|
- si = G_alloc_matrix(my, mx);
|
|
|
-
|
|
|
- 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;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
+ /* 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);
|
|
|
+
|
|
|
+ /* 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);
|
|
|
+
|
|
|
+ /* 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;
|
|
|
}
|
|
|
- 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);
|
|
|
|
|
|
- /****************/
|
|
|
+ /* Update elevation map */
|
|
|
+ copy_matrix_undef_double_to_float_values(rows, cols, si, zz);
|
|
|
+
|
|
|
+ /* 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);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Traps */
|
|
|
+ if (traps != NULL)
|
|
|
+ trap = read_float_raster_map(rows, cols, traps, 1.0);
|
|
|
+ else
|
|
|
+ trap = create_float_matrix(rows, cols, (double)UNDEF);
|
|
|
+ }
|
|
|
|
|
|
- if (detin)
|
|
|
- Rast_close(detin_fd);
|
|
|
+ if (detin != NULL) {
|
|
|
+ dc = read_float_raster_map(rows, cols, detin, 1.0);
|
|
|
+ copy_matrix_undef_float_values(rows, cols, dc, zz);
|
|
|
+ }
|
|
|
|
|
|
- if (tranin)
|
|
|
- Rast_close(trainin_fd);
|
|
|
+ if (tranin != NULL) {
|
|
|
+ ct = read_float_raster_map(rows, cols, tranin, 1.0);
|
|
|
+ copy_matrix_undef_float_values(rows, cols, ct, zz);
|
|
|
+ }
|
|
|
|
|
|
- if (tauin)
|
|
|
- Rast_close(tauin_fd);
|
|
|
+ if (tauin != NULL) {
|
|
|
+ tau = read_float_raster_map(rows, cols, tauin, 1.0);
|
|
|
+ copy_matrix_undef_float_values(rows, cols, tau, zz);
|
|
|
+ }
|
|
|
|
|
|
- if (wdepth)
|
|
|
- Rast_close(wdepth_fd);
|
|
|
+ if (wdepth != NULL) {
|
|
|
+ gama = read_double_raster_map(rows, cols, wdepth, 1.0);
|
|
|
+ copy_matrix_undef_double_to_float_values(rows, cols, gama, zz);
|
|
|
+ }
|
|
|
+
|
|
|
+ /* Array for gradient checking */
|
|
|
+ slope = create_double_matrix(rows, cols, 0.0);
|
|
|
+
|
|
|
+ /* Create the observation points and open the logfile */
|
|
|
+ create_observation_points();
|
|
|
|
|
|
- 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 @@ int grad_check(void)
|
|
|
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 @@ int grad_check(void)
|
|
|
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 @@ int grad_check(void)
|
|
|
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;
|
|
|
- }
|
|
|
- else {
|
|
|
- res = arg2;
|
|
|
+ for(row = 0; row < rows; row++) {
|
|
|
+ for(col = 0; col < cols; col++) {
|
|
|
+ if(source[row][col] == UNDEF)
|
|
|
+ target[row][col] = UNDEF;
|
|
|
+ }
|
|
|
}
|
|
|
-
|
|
|
- return res;
|
|
|
}
|
|
|
|
|
|
+/* ************************************************************************* */
|
|
|
|
|
|
-double amin1(double arg1, double arg2)
|
|
|
+void copy_matrix_undef_float_values(int rows, int cols, float **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;
|
|
|
+}
|
|
|
+
|
|
|
+/* ************************************************************************* */
|
|
|
+
|
|
|
+float ** create_float_matrix(int rows, int cols, float fill_value)
|
|
|
+{
|
|
|
+ int col = 0, row = 0;
|
|
|
+ float **matrix = NULL;
|
|
|
+
|
|
|
+ /* 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;
|
|
|
+ }
|
|
|
}
|
|
|
|
|
|
- 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;
|
|
|
- }
|
|
|
- else {
|
|
|
- res = arg2;
|
|
|
+ /* 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;
|
|
|
+ }
|
|
|
}
|
|
|
|
|
|
- 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;
|
|
|
- }
|
|
|
- else {
|
|
|
- res = arg2;
|
|
|
+ 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;
|
|
|
+ }
|
|
|
}
|
|
|
|
|
|
- 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;
|
|
|
+}
|
Soeren Gebbert