[GRASS-SVN] r71171 - grass-addons/grass7/vector/v.kriging
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 8 03:32:50 PDT 2017
Author: evas
Date: 2017-06-08 03:32:49 -0700 (Thu, 08 Jun 2017)
New Revision: 71171
Modified:
grass-addons/grass7/vector/v.kriging/Makefile
grass-addons/grass7/vector/v.kriging/geostat.c
grass-addons/grass7/vector/v.kriging/getval.c
grass-addons/grass7/vector/v.kriging/local_proto.h
grass-addons/grass7/vector/v.kriging/main.c
grass-addons/grass7/vector/v.kriging/quantile.c
grass-addons/grass7/vector/v.kriging/spatial_index.c
grass-addons/grass7/vector/v.kriging/utils.c
grass-addons/grass7/vector/v.kriging/utils_kriging.c
grass-addons/grass7/vector/v.kriging/utils_raster.c
grass-addons/grass7/vector/v.kriging/utils_write.c
grass-addons/grass7/vector/v.kriging/v.kriging.html
Log:
v.kriging: optimized kriging algorithm (2D and 3D); added pdf with preliminary case studies
Modified: grass-addons/grass7/vector/v.kriging/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/Makefile 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/Makefile 2017-06-08 10:32:49 UTC (rev 71171)
@@ -1,12 +1,12 @@
-MODULE_TOPDIR = ../..
+MODULE_TOPDIR = ../../../..
PGM = v.kriging
-
+
LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(GMATHLIB) $(IOSTREAMLIB) $(MATHLIB) $(RTREELIB) $(RASTER3DLIB) $(RASTERLIB)
DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(GISDEP) $(GMATHDEP) $(IOSTREAMDEP) $(RTREEDEP) $(RASTER3DDEP) $(RASTERDEP)
EXTRA_INC = $(VECT_INC) $(PROJINC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
-include $(MODULE_TOPDIR)/include/Make/Module.make
+include $(MODULE_TOPDIR)/include/Make/Module.make
default: cmd
Modified: grass-addons/grass7/vector/v.kriging/geostat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/geostat.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/geostat.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -7,15 +7,18 @@
{
// Variogram properties
struct parameters *var_pars;
+ int direction;
double max_dist, max_dist_vert; // max distance of interpolated points
double radius, radius_vert; // radius of interpolation in the horizontal plane
switch (type) {
case 0: // horizontal variogram
var_pars = &pars->hz; // 0: horizontal variogram
+ direction = 12;
break;
- case 1: // vertica variogram
- var_pars = &pars->vert; // 1: vertical variogram
+ case 1: // vertical variogram
+ var_pars = &pars->vert; // 1: vertical variogram
+ direction = 3;
break;
case 2: // bivariate variogram
var_pars = &pars->fin;
@@ -28,6 +31,7 @@
var_pars = &pars->fin; // default: final variogram
max_dist = max_dist_vert = var_pars->max_dist;
radius = radius_vert = SQUARE(max_dist);
+ direction = 0;
break;
}
@@ -73,7 +77,7 @@
int nrows = nLag;
int ncols = type == 2 ? nLag_vert : 1;
- struct write *report = &xD->report;
+ struct write *report = xD->report;
// Variogram processing variables
int s; // index of horizontal segment
@@ -133,11 +137,7 @@
// control initialization:
if (dr == NULL || search == NULL || var_pars->h == NULL ||
(type == 2 && var_pars->vert == NULL)) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
- "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
- }
+ report_error(xD->report);
G_fatal_error(_("Memory allocation failed..."));
}
@@ -146,11 +146,7 @@
gamma_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix (vector) of gammas
if (c_M == NULL || gamma_M == NULL) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
- "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
- }
+ report_error(xD->report);
G_fatal_error(_("Memory allocation failed..."));
}
@@ -180,9 +176,9 @@
for (s = 0; s < nrows; s++) { // for each horizontal lag (isotrophy!!!)...
*h = (s + 1) * lag; // lag distance
- r = &pnts->r[0]; // pointer to the input coordinates
// for every bs cycle ...
+ r = &pnts->r[0]; // pointer to the input coordinates
i_vals = &vals[0]; // pointer to input values to be interpolated
gamma_lag = 0.; // gamma in dir direction and h distance
cpls = 0.; // count of couples in dir direction and h distance
@@ -206,13 +202,13 @@
n_vals = list->n_values; // # of input values located on NN
if (n_vals > 0) {
- correct_indices(list, r, pnts, var_pars);
+ correct_indices(direction, list, r, pnts, var_pars);
ii = &list->value[0]; // indices of these input values (note: increased by 1)
j_vals = &vals[*ii]; // pointer to input values
for (j = 1; j < n_vals; j++) { // for each point within overlapping rectangle
if (*ii > i) { // use the points just once
- coord_diff(i, *ii, pnts_r, dr); // compute coordinate differences
+ coord_diff(i, *ii, pnts_r, dr); // compute coordinate differences
// Variogram processing
if (type == 1) { // vertical variogram:
@@ -221,15 +217,12 @@
}
else { // hz / aniso / bivar variogram:
tv = atan2(*(dr + 1), *dr); // bearing
- if (tv < 0.) {
- tv += 2. * PI;
- }
}
- ddir1 = dir - tv; // difference between bearing and azimuth
- ddir2 = (dir + PI) - tv;
+ ddir1 = tv - dir; // difference between bearing and azimuth
+ ddir2 = tv + (PI - dir); // reverse
- if (fabs(ddir1) <= td || fabs(ddir2) <= td) { // angle test: compare the diff with critical value
+ if (fabsf(ddir1) <= td || fabsf(ddir2) <= td) { // angle test: compare the diff with critical value
// test squared distance: vertical variogram => 0., ...
rv = type == 1 ? 0. : radius_hz_diff(dr); // ... otherwise horizontal distance
@@ -237,9 +230,11 @@
rv += SQUARE(*(dr + 2)); // consider also vertical direction
}
- rvh = sqrt(rv) - *h; // the difference between distance and lag
- if (rv <= radius && fabs(rvh) <= lag) { // distance test: compare the distance with critical value and find out if the j-point is located within i-lag
- if (type == 2) { // vertical test for bivariate variogram:
+ rvh = *h - sqrt(rv); // the difference between distance and lag
+ // distance test: compare the distance with critical value: find out if the j-point is located within i-lag
+ if (rv <= radius && 0. <= rvh && rvh <= lag) { // 0. <= rvh && rvh <= lag
+ // vertical test for bivariate variogram:
+ if (type == 2) {
rvh = *(dr + 2) - *vert; // compare vertical
if (fabs(rvh) <= lag_vert) { // elevation test: vertical lag
@@ -248,7 +243,7 @@
else {
goto end;
}
- }
+ } // end if: type == 2
delta_V:
dv = *j_vals - *i_vals; // difference of values located on pair of points i, j
gamma_lag += SQUARE(dv); // sum of squared differences
@@ -262,16 +257,12 @@
} // end j for loop: points within overlapping rectangles
} // end test: n_vals > 0
else {
- if (report->name) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("This point does not have neighbours in given radius..."));
}
-
r += 3; // go to the next search point
i_vals++; // and to its value
+ G_free_ilist(list); // free list memory
} // end for loop i: variogram computation for each i-th search point
if (isnan(gamma_lag) || cpls == 0.0) { // empty lags:
@@ -289,6 +280,7 @@
h++;
c++;
gamma++;
+ //G_fatal_error(_("end 1st loop"));
} // end for loop s
if (type == 2) { // vertical variogram:
@@ -297,11 +289,7 @@
} // end for loop b
if (err0 == nLag) { // todo> kedy nie je riesitelny teoreticky variogram?
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Unable to compute experimental variogram..."));
} // end error
@@ -315,26 +303,32 @@
sill(var_pars); // compute sill
}
- if (report->write2file == TRUE) { // report file available:
+ if (report->name) { // report file available:
write2file_variogram_E(xD, var_pars, c_M); // write to file
}
write_temporary2file(xD, var_pars);
- G_free_ilist(list); // free list memory
G_matrix_free(c_M);
G_matrix_free(gamma_M);
}
/* theoretical variogram */
-void T_variogram(int type, int i3, struct opts opt,
- struct parameters *var_pars, struct write *report)
+void T_variogram(int type, struct opts opt,
+ struct parameters *var_pars, struct int_par *xD)
{
+ int i3 = xD->i3;
+ struct write *report;
+
+ report = xD->report;
+
char *variogram;
// report
- if (report->write2file == TRUE) { // report file available:
+ if (report->name) { // report file available:
+ report->fp = fopen(report->name, "a");
time(&report->now); // write down time of start
+ // just for the beginning of the phase => type 0 or 2
if (type != 1) {
fprintf(report->fp,
"\nComputation of theoretical variogram started on %s\n",
@@ -347,27 +341,34 @@
var_pars->const_val = 0; // input values are not constants
switch (type) {
- case 0: // horizontal variogram
- if (i3 == TRUE) { // 3D interpolation (middle phase)
+ // horizontal variogram
+ case 0:
+ // 3D interpolation (middle phase)
+ if (i3 == TRUE) {
variogram = opt.function_var_hz->answer; // function type available:
var_pars->function = set_function(variogram, report);
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
+ // nonlinear or not parabolic
+ if (strcmp(variogram, "linear") != 0) {
var_pars->nugget = atof(opt.nugget_hz->answer);
var_pars->h_range = atof(opt.range_hz->answer);
if (opt.sill_hz->answer) {
var_pars->sill = atof(opt.sill_hz->answer);
}
}
- else { // function type not available:
+ // function type not available:
+ else {
LMS_variogram(var_pars, report);
}
}
- else { // 2D interpolation (final phase)
+
+ // 2D interpolation (final phase)
+ else {
variogram = opt.function_var_final->answer; // function type available:
var_pars->function = set_function(variogram, report);
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
+ // nonlinear or not parabolic
+ if (strcmp(variogram, "linear") != 0) {
var_pars->nugget = atof(opt.nugget_final->answer);
var_pars->h_range = atof(opt.range_final->answer);
if (opt.sill_final->answer) {
@@ -379,7 +380,8 @@
}
}
- if (report->name) {
+ if (report->name && strcmp(variogram, "linear") != 0) {
+ report->fp = fopen(report->name, "a");
fprintf(report->fp, "Parameters of horizontal variogram:\n");
fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
fprintf(report->fp, "Sill: %f\n", var_pars->sill);
@@ -387,23 +389,28 @@
}
break;
- case 1: // vertical variogram
+ // vertical variogram
+ case 1:
var_pars->nugget = atof(opt.nugget_vert->answer);
var_pars->h_range = atof(opt.range_vert->answer);
+
if (opt.sill_vert->answer) {
var_pars->sill = atof(opt.sill_vert->answer);
}
+
variogram = opt.function_var_vert->answer;
var_pars->function = set_function(variogram, report);
- if (report->name) {
+ if (report->name && strcmp(variogram, "linear") != 0) {
fprintf(report->fp, "Parameters of vertical variogram:\n");
fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
fprintf(report->fp, "Sill: %f\n", var_pars->sill);
fprintf(report->fp, "Range: %f\n", var_pars->h_range);
}
break;
- case 2: // bivariate variogram (just final phase)
+
+ // bivariate variogram (just final phase)
+ case 2:
if (!(opt.function_var_final->answer && opt.function_var_final_vert->answer) || strcmp(opt.function_var_final->answer, "linear") == 0) { // planar function:
var_pars->function = 5; // planar variogram (3D)
LMS_variogram(var_pars, report);
@@ -448,21 +455,22 @@
}
}
- plot_var(i3, TRUE, var_pars); // Plot variogram using gnuplot
+ plot_var(xD, TRUE, var_pars); // Plot variogram using gnuplot
break;
case 3: // univariate (just final phase)
variogram = opt.function_var_final->answer;
var_pars->function = set_function(variogram, report);
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear and not parabolic variogram:
+ // nonlinear and not parabolic variogram:
+ if (strcmp(variogram, "linear") != 0) {
var_pars->nugget = atof(opt.nugget_final->answer);
var_pars->h_range = atof(opt.range_final->answer);
if (opt.sill_final->answer) {
var_pars->sill = atof(opt.sill_final->answer);
}
variogram = opt.function_var_final->answer;
- if (report->write2file == TRUE) {
+ if (report->name) {
if (i3 == TRUE) { // 3D interpolation:
fprintf(report->fp,
"Parameters of anisotropic variogram:\n");
@@ -475,6 +483,7 @@
fprintf(report->fp, "Range: %f\n", var_pars->h_range);
}
}
+ // linear variogram:
else { // linear or parabolic variogram
LMS_variogram(var_pars, report);
}
@@ -482,7 +491,7 @@
}
if (type != 2) {
- plot_var(i3, FALSE, var_pars); // Plot variogram using gnuplot
+ plot_var(xD, FALSE, var_pars); // Plot variogram using gnuplot
}
}
@@ -490,12 +499,10 @@
struct points *pnts, struct var_par *pars,
struct output *out)
{
- G_fatal_error(_("Interpolating values is currently under maintenance (optimization). Theoretical variogram of your data has been computed."));
// Local variables
int i3 = xD->i3;
- double *vals = pnts->invals; // values to be used for interpolation
- struct write *report = &xD->report;
- struct write *crossvalid = &xD->crossvalid;
+ struct write *report = xD->report;
+ struct write *crossvalid = xD->crossvalid;
struct parameters *var_par = &pars->fin;
int type = var_par->type;
@@ -508,30 +515,37 @@
unsigned int percents = 50; // counter
unsigned int row, col, dep; // indices of cells/voxels
- double rslt_OK; // interpolated value located on r0
+ int resample;
+ struct krig_pars krig;
+ int add_trend = (out->trend[0] == 0. && out->trend[1] == 0. &&
+ out->trend[2] == 0. &&
+ out->trend[3] == 0.) ? FALSE : TRUE;
pnts->max_dist = var_par->lag;
- struct ilist *list;
+ struct ilist *list, *list_new;
- double *r0; // xyz coordinates of cell/voxel centre
- mat_struct *GM;
- mat_struct *GM_sub; // submatrix of selected points
- mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
- mat_struct *g0; // diffences between known and unknown values = theor_var(dist)
- mat_struct *w0; // weights of values located on the input points
+ double *r0, rslt, trend_val; // xyz coordinates of cell/voxel centre
// Cell/voxel center coords (location of interpolated value)
r0 = (double *)G_malloc(3 * sizeof(double));
- if (report->write2file) { // report file available:
+ int *row_init, *row0, count, next = 1, total, complete, new_matrix =
+ 0, setup, mat_row;
+ row_init = (int *)G_malloc(reg->ncols * sizeof(int));
+ row0 = &row_init[0];
+
+ int i, ndeps = reg->ndeps, nrows = reg->nrows, ncols = reg->ncols;
+
+ krig.rslt = G_matrix_init(nrows * ndeps, ncols, nrows * ndeps);
+ krig.first = TRUE;
+
+ if (report->name) { // report file available:
time(&report->now);
fprintf(report->fp, "Interpolating values started on %s\n\n",
ctime(&report->now));
fflush(report->fp);
}
- G_message(_("Interpolating unknown values..."));
- G_fatal_error(_("... is currently under maintenance (optimization). Theoretical variogram of your data has been computed."));
if (percents) {
G_percent_reset();
}
@@ -542,111 +556,252 @@
goto constant_voxel_centre;
}
- GM = set_up_G(pnts, var_par, &xD->report); // set up matrix of dissimilarities of input points
- var_par->GM = G_matrix_copy(GM); // copy matrix because of cross validation
+ set_up_G(pnts, var_par, xD->report, &krig); // set up matrix of dissimilarities of input points
+ var_par->GM = G_matrix_copy(krig.GM); // copy matrix because of cross validation
// perform cross validation...
- if (crossvalid->write2file) { // ... if desired
- crossvalidation(xD, pnts, var_par);
+ if (crossvalid->name) { // ... if desired
+ crossvalidation(xD, pnts, var_par, reg);
}
+ G_message(_("Interpolating unknown values..."));
constant_voxel_centre:
- for (dep = 0; dep < reg->ndeps; dep++) {
- if (xD->i3 == TRUE) {
- if (percents) {
- G_percent(dep, reg->ndeps, 1);
+ count = 0;
+ col = row = dep = 0;
+ total = ndeps * nrows * ncols;
+
+ while (count <= total) {
+ if (percents) {
+ G_percent(count, total, 1);
+ }
+
+ if (var_par->const_val == 1) { // constant input values
+ goto constant_voxel_val;
+ }
+
+ // coordinates of output point (center of the pixel / voxel)
+ if (next < 2) {
+ cell_centre(col, row, dep, xD, reg, r0, var_par);
+
+ // initial subsample
+ if (col == 0 && row == 0 && dep == 0) {
+ list = list_NN(xD, r0, pnts, max_dist, max_dist_vert);
+ make_subsamples(xD, list, r0, dep * nrows + row, col, pnts,
+ var_par, &krig);
}
+ // lists to compare
+ else if (krig.new == FALSE) {
+ list_new = list_NN(xD, r0, pnts, max_dist, max_dist_vert);
+ next = compare_NN(list, list_new, krig.modified);
+ }
}
- for (row = 0; row < reg->nrows; row++) {
- if (xD->i3 == FALSE) {
- if (percents) {
- G_percent(row, reg->nrows, 1);
- }
+
+ if (next > 0) {
+ rslt = interpolate(xD, list, r0, pnts, var_par, &krig);
+ if (add_trend == TRUE) {
+ trend_val = trend(r0, out, var_par->function, xD);
+ rslt += trend_val;
}
- //#pragma omp parallel for private(col, r0, GM, GM_Inv, g0, w0, rslt_OK)
- for (col = 0; col < reg->ncols; col++) {
+ mat_row = dep * nrows + row;
- if (var_par->const_val == 1) { // constant input values
- goto constant_voxel_val;
- }
+ setup = G_matrix_set_element(krig.rslt, mat_row, col, rslt);
- cell_centre(col, row, dep, xD, reg, r0, var_par); // coords of output point
+ if (G_matrix_get_element(krig.rslt, 0, 0) == 0. && (row != 0))
+ G_fatal_error(_("%d %d %d %d %f %f"), dep, row, col,
+ mat_row, G_matrix_get_element(krig.rslt, 0, 0),
+ krig.rslt->vals[(dep * nrows + row) * ncols +
+ col]);
- // add cell centre to the R-tree
- list = G_new_ilist(); // create list of overlapping rectangles
+ if (setup < 0) {
+ report_error(report);
+ G_fatal_error(_("The value %f was not written to the cell %d %d %d"),
+ rslt, dep, row, col);
+ }
+ count++;
+ krig.new = FALSE;
- if (i3 == TRUE) { // 3D kriging:
- list =
- find_NNs_within(3, r0, pnts, max_dist, max_dist_vert);
+ next = next == 2 ? 1 : next; // mark it as coincident
+
+ // Create output
+ constant_voxel_val:
+ if (var_par->const_val == 1) { // constant input values:
+ rslt = (double)*pnts->invals; // setup input as output
+ }
+ } // end if next > 0
+
+ // go to the next cell point:
+ switch (next) {
+ // the cell was not interpolated:
+ case 0:
+ // new column:
+ if (row == 0 || row == *row0) {
+ Vect_reset_list(list); // refresh list
+ Vect_list_append_list(list, list_new); // add new list
+ G_free_ilist(list_new); // free list memory
+ krig.modified = 1;
+ krig.new = TRUE;
+ next = 2; // interpolate using new subsample
+ G_matrix_free(krig.GM_Inv); // free old matrix
+ new_matrix++; // counter of skipped matrices
+ make_subsamples(xD, list, r0, dep * nrows + row, col, pnts,
+ var_par, &krig);
+ }
+ // the same column:
+ else {
+ *row0 = row; // save row index to continue
+
+ // general cells:
+ if (col < ncols - 1) {
+ col++; // go to the next row
+ row0++; // go to the col index in the next row
+ while (*row0 == nrows) {
+ if (col < ncols - 1) {
+ col++; // go to the next row
+ row0++; // go to the col index in the next row
+ }
+ else {
+ krig.first = FALSE;
+ col = 0;
+ row0 = &row_init[0];
+ }
+ }
}
- else { // 2D kriging:
- list =
- find_NNs_within(2, r0, pnts, max_dist, max_dist_vert);
- }
+ // last column:
+ else {
+ krig.first = FALSE;
+ col = 0; // start new sampling
+ row0 = &row_init[0]; // from the beginning;
- if (list->n_values > 1) { // positive # of selected points:
- correct_indices(list, r0, pnts, var_par);
+ // skip full cols:
+ while (*row0 == nrows) { // full row
+ if (col < ncols - 1) {
+ col++; // go to the next row
+ row0++; // to test if it is completed
+ }
+ } // end while: full row
+ } // end else: last column
+ } // end else: the same column
- GM_sub = submatrix(list, GM, report); // make submatrix for selected points
- GM_Inv = G_matrix_inverse(GM_sub); // invert submatrix
- G_matrix_free(GM_sub);
+ row = krig.first == TRUE ? 0 : *row0; // setup starting cell
+ break; // count does not rise
- g0 = set_up_g0(xD, pnts, list, r0, var_par); // Diffs inputs - unknowns (incl. cond. 1))
- w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
+ // the cell was interpolated -> examine the next one
+ case 1:
+ // save index to continue:
+ row++; // continue in the row
- G_matrix_free(GM_Inv);
- G_matrix_free(g0);
+ // last row (full column):
+ if (row == nrows) {
+ *row0 = row; // save the number to distinguish full columns
- rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
- G_matrix_free(w0);
- }
- else if (list->n_values == 1) {
- rslt_OK = vals[list->value[0] - 1]; // Estimated cell/voxel value rslt_OK = w x inputs
- }
- else if (list->n_values == 0) {
- if (report->write2file) { // report file available:
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("This point does not have neighbours in given radius..."));
- } // end else: error
+ // the last column
+ if (col == ncols - 1) {
+ // start from the beginning
+ krig.first = FALSE;
+ col = 0;
+ row0 = &row_init[0];
- G_free_ilist(list); // free list memory
+ // skip full columns
+ while (*row0 == nrows) {
+ if (col < ncols - 1) {
+ col++;
+ row0++;
+ }
+ else {
+ if (i3 == FALSE || dep == ndeps - 1) {
+ if (dep > 1) {
+ G_message(_("Vertical level %d has been processed..."),
+ dep + 1);
+ }
+ goto accomplished;
+ }
+ else {
+ dep++; // next vertical level
+ col = 0; // the first column
- // Create output
- constant_voxel_val:
- if (var_par->const_val == 1) { // constant input values:
- rslt_OK = (double)*vals; // setup input as output
- }
+ // refresh rows
+ new_vertical(row_init, ncols); // row0 = 0
+ row0 = &row_init[0];
+ row = *row0;
- // write output to the (3D) raster layer
- if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
- if (report->write2file) { // report file available
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp); // close report file
- }
- G_fatal_error(_("Error writing result into output layer..."));
- }
- } // end col
- } // end row
- } // end dep
+ next =
+ new_sample(xD, list, list_new, pnts, dep,
+ row, col, r0, max_dist,
+ max_dist_vert, reg, var_par,
+ &krig, &new_matrix);
+ goto new_dep;
+ } // end else: complete level
+ }
+ } // end if: *row0 == nrows
+ } // end if: last column
- if (report->write2file) {
+ else {
+ col++;
+ row0++; // read starting row
+ complete = 0;
+
+ while (*row0 == nrows) {
+ if (col < ncols - 1) {
+ col++;
+ row0++;
+ complete++;
+ }
+ else {
+ krig.first = FALSE;
+ col = 0;
+ row0 = &row_init[0];
+ complete = 0;
+ }
+ if (complete == ncols - 1) {
+ if (i3 == FALSE || dep == ndeps - 1) {
+ if (dep > 1) {
+ G_message(_("Vertical level %d has been processed..."),
+ dep + 1);
+ }
+ goto accomplished;
+ }
+ else {
+ dep++; // next vertical level
+ col = 0; // the first column
+
+ // refresh rows
+ new_vertical(row_init, ncols); // row0 = 0
+ row0 = &row_init[0];
+ row = *row0;
+
+ next =
+ new_sample(xD, list, list_new, pnts, dep,
+ row, col, r0, max_dist,
+ max_dist_vert, reg, var_par,
+ &krig, &new_matrix);
+ goto new_dep;
+ } // end else: complete level
+ } // if: complete == ncols - 1
+ } // while: *row0 == nrows
+ } // else: general row
+ row = krig.first == TRUE ? 0 : *row0;
+ } // end else: go to the new column
+
+ break;
+ } // end switch next
+ new_dep:
+ if (row == 0 && col == 0) {
+ G_message(_("Vertical level %d has been processed..."), dep);
+ }
+ } // end while
+
+ accomplished:
+ G_message(_("# of points: %d # of matrices: %d diff: %d"), total,
+ new_matrix, total - new_matrix);
+
+ // write output to the (3D) raster layer
+ write2layer(xD, reg, out, krig.rslt);
+
+ if (report->name) {
fprintf(report->fp,
"\n************************************************\n\n");
time(&report->now);
fprintf(report->fp, "v.kriging completed on %s", ctime(&report->now));
fclose(report->fp);
}
-
- switch (xD->i3) {
- case TRUE:
- Rast3d_close(out->fd_3d); // Close 3D raster map
- break;
- case FALSE:
- Rast_close(out->fd_2d); // Close 2D raster map
- break;
- }
}
Modified: grass-addons/grass7/vector/v.kriging/getval.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/getval.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/getval.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -7,8 +7,7 @@
/* get array of values from attribute column (function based on part of v.buffer2 (Radim Blazek, Rosen Matev)) */
double *get_col_values(struct Map_info *map, struct int_par *xD,
- struct points *pnts, int field, const char *column,
- int detrend)
+ struct points *pnts, int field, const char *column)
{
struct select *in_reg = &pnts->in_reg;
@@ -35,10 +34,10 @@
db_CatValArray_init(&cvarr); /* array of categories and values initialised */
Fi = Vect_get_field(map, field); /* info about call of DB */
if (Fi == NULL) {
- if (save == 1 && xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
+ if (save == 1 && xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error(_("Database connection not defined for layer %d"),
field);
@@ -46,10 +45,10 @@
Driver = db_start_driver_open_database(Fi->driver, Fi->database); /* started connection to DB */
if (Driver == NULL) {
- if (save == 1 && xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
+ if (save == 1 && xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
Fi->database, Fi->driver);
@@ -64,20 +63,20 @@
db_select_CatValArray(Driver, Fi->table, Fi->key, column, NULL,
&cvarr);
if (nrec < 0) {
- if (save == 1 && xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
+ if (save == 1 && xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error(_("Unable to select data from table <%s>"), Fi->table);
}
ctype = cvarr.ctype;
if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE) {
- if (save == 1 && xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
+ if (save == 1 && xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error(_("Column must be numeric"));
}
@@ -116,65 +115,11 @@
}
}
- if (detrend == 1) {
- double *x, *y, *z;
-
- x = (double *)G_malloc(n * sizeof(double));
- y = (double *)G_malloc(n * sizeof(double));
- z = (double *)G_malloc(n * sizeof(double));
-
- mat_struct *A, *gR, *T, *res;
-
- x = &pnts->r[0];
- y = &pnts->r[1];
- z = &pnts->r[2];
- vals = &values[0];
- // Number of columns of design matrix A
- A = G_matrix_init(n, 2, n); // initialise design matrix, normal grav: 7
- gR = G_matrix_init(n, 1, n); // initialise vector of observations
-
- for (i = 0; i < n; i++) {
- G_matrix_set_element(A, i, 0, *x);
- G_matrix_set_element(A, i, 1, *y);
- G_matrix_set_element(gR, i, 0, *vals);
- x += 3;
- y += 3;
- z += 3;
- vals++;
- }
- T = LSM(A, gR); // Least Square Method
- res = G_matrix_product(A, T);
-
- FILE *fp;
-
- x = &pnts->r[0];
- y = &pnts->r[1];
- z = &pnts->r[2];
- fp = fopen("trend.txt", "w");
-
- vals = &values[0];
- double *resid;
-
- resid = &res->vals[0];
- for (i = 0; i < n; i++) {
- *vals = *vals - *resid;
- fprintf(fp, "%f %f %f %f\n", *x, *y, *z, *vals);
- x += 3;
- y += 3;
- z += 3;
- vals++;
- resid++;
- }
- fclose(fp);
- //G_debug(0,"a=%f b=%f c=%f d=%f", T->vals[0], T->vals[1], T->vals[2], T->vals[3]);
- pnts->trend = T;
+ if (xD->phase == 0 && xD->report->name) {
+ write2file_values(xD->report, column);
+ test_normality(n, values, xD->report);
}
- if (xD->report.write2file && xD->phase == 0) {
- write2file_values(&xD->report, column);
- test_normality(n, values, &xD->report);
- }
-
db_CatValArray_free(&cvarr); // free memory of the array of categories and values
return values;
@@ -227,7 +172,7 @@
// Get 3rd coordinate of 2D points from attribute column -> 3D interpolation
if (xD->v3 == FALSE && zcol != NULL) {
- z_attr = get_col_values(map, NULL, point, field, zcol, FALSE);
+ z_attr = get_col_values(map, NULL, point, field, zcol);
}
nskipped = 0; // # of skipped elements (lines, areas)
@@ -236,7 +181,7 @@
while (TRUE) {
type = Vect_read_next_line(map, Points, NULL);
if (type == -1) {
- if (report->write2file == TRUE) { // close report file
+ if (report->name) { // close report file
fprintf(report->fp,
"Error (see standard output). Process killed...");
fclose(report->fp);
@@ -254,7 +199,7 @@
}
if (isnan(Points->x[0]) || isnan(Points->y[0]) || isnan(Points->z[0])) {
- if (report->write2file == TRUE) { // close report file
+ if (report->name) { // close report file
fprintf(report->fp,
"Error (see standard output). Process killed...");
fclose(report->fp);
@@ -277,11 +222,7 @@
*rz = Points->z[0]; // set up z coordinate
if (*rz != Points->z[0]) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Error reading input coordinates z..."));
}
}
@@ -356,7 +297,7 @@
Vect_destroy_line_struct(Points);
if (n_in_reg == 0) {
- if (report->write2file == TRUE) { // close report file
+ if (report->name) { // close report file
fprintf(report->fp,
"Error (see standard output). Process killed...");
fclose(report->fp);
@@ -368,7 +309,7 @@
G_message(_("Unused points: %d (out of region)"), out_reg);
}
- if (xD->report.write2file == TRUE && xD->phase == 0) { // initial phase:
+ if (xD->phase == 0 && xD->report->name) { // initial phase:
write2file_vector(xD, point); // describe properties
}
}
@@ -395,10 +336,10 @@
Rast3d_get_window(®->reg_3d); /* stores the current default window in region */
if (reg->reg_3d.bottom == 0 &&
(reg->reg_3d.tb_res == 0 || reg->reg_3d.depths == 0)) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp,
+ if (xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error
("To process 3D interpolation, please set 3D region settings.");
@@ -425,54 +366,55 @@
{
FILE *fp;
- int j, nLag_vert;
- double lag_vert, max_dist_vert;
- double *v_elm, *g_elm;
- double sill_hz, sill_vert;
-
-
- fp = fopen(file_name, "r");
+ fp = fopen(file_name, "r"); // open file
if (fp == NULL) {
G_fatal_error(_("Temporary file <%s> is missing, please repeat an initial phase..."),
file_name);
}
- else { // file exists:
- int i, type;
- int nLag;
- double lag, max_dist, td_hz, sill;
- double *h_elm, *gamma;
- int file, file_length;
+ // file exists:
+ else {
+ int i, type; // index, type of variogram
+ int nLag; // # of bins
+ double lag, max_dist, td_hz, sill; // size of the bin, maximum horizontal distance
+ double *h, *h_elm, *gamma; // elements of horizontal bins and gamma matrices
+ int file, file_length; // type of file (9/8), length of the filename
- for (i = 0; i < 2; i++) {
- if (fscanf(fp, "%d", &file_length) == 0) { // filename length
+ int j, nLag_vert; // index and # of vertical bins
+ double lag_vert, max_dist_vert; // size of the lag and maximum vertical distance
+ double *v_elm, *g_elm; // elements of vertical bins and gamma matrices
+ double sill_hz, sill_vert; // horizontal and vertical sill
+
+ // check filename length (1st variable in the file)
+ if (fscanf(fp, "%d", &file_length) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ // valid filename
+ if (file_length > 3) {
+ // read file code (Y/N report)
+ if (fscanf(fp, "%d", &file) == 0) {
G_fatal_error(_("Nothing to scan..."));
}
- if (file_length > 3) {
- if (fscanf(fp, "%d", &file) == 0) {
+
+ // read report name
+ if (file == 9) { // report
+ xD->report->name =
+ (char *)G_malloc((file_length + 1) * sizeof(char));
+ if (fscanf(fp, "%s", xD->report->name) == 0) { // filename
G_fatal_error(_("Nothing to scan..."));
}
- if (file == 9) { // filetype code (9 - report, 8 - crossval)
- xD->report.name =
- (char *)G_malloc(file_length * sizeof(char));
- if (fscanf(fp, "%s", xD->report.name) == 0) { // filename
- G_fatal_error(_("Nothing to scan..."));
- }
- continue;
+
+ // middle / final phase: check if file does not exist
+ if (access(xD->report->name, F_OK) != 0) {
+ G_fatal_error(_("Report file does not exist; please check the name or repeat initial phase..."));
}
- else if (file == 8) {
- xD->crossvalid.name =
- (char *)G_malloc(file_length * sizeof(char));
- if (fscanf(fp, "%s", xD->crossvalid.name) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- continue;
- }
+ xD->report->fp = fopen(xD->report->name, "a");
}
- else {
- type = file_length;
- goto no_file;
- }
+ }
+ else {
+ type = file_length;
+ goto no_file;
} // todo: test without report and crossval
if (fscanf(fp, "%d", &type) == 0) { // read type
@@ -480,9 +422,9 @@
}
no_file:
-
switch (type) {
- case 2: // bivariate variogram
+ // bivariate variogram
+ case 2:
var_par->type = type;
if (fscanf
(fp, "%d %lf %lf %d %lf %lf %lf", &nLag_vert, &lag_vert,
@@ -536,26 +478,42 @@
break;
default:
- if (type == 3) { // anisotropic variogram:
- if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0) { // anisotropic ratio
+ // anisotropic variogram:
+ if (type == 3) {
+ // anisotropic ratio
+ if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0) {
G_fatal_error(_("Nothing to scan..."));
}
}
+ // distance settings
if (fscanf(fp, "%d %lf %lf", &nLag, &lag, &max_dist) < 3) {
G_fatal_error(_("Nothing to scan..."));
}
+ // angular settings
if (type != 1) {
if (fscanf(fp, "%lf", &td_hz) == 0) {
G_fatal_error(_("Nothing to scan..."));
}
}
- var_par->h = (double *)G_malloc(nLag * sizeof(double));
- h_elm = &var_par->h[0];
+ // empirical variogram:
+ h = (double *)G_malloc(nLag * sizeof(double));
+ // control initialization:
+ if (h == NULL) {
+ report_error(xD->report);
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+ h_elm = &h[0];
var_par->gamma = G_matrix_init(nLag, 1, nLag);
+ // control initialization:
+ if (var_par->gamma == NULL) {
+ report_error(xD->report);
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
gamma = &var_par->gamma->vals[0];
for (i = 0; i < nLag; i++) {
@@ -573,6 +531,10 @@
break;
}
+ var_par->h = (double *)G_malloc(nLag * sizeof(double));
+ memcpy(var_par->h, h, nLag * sizeof(double));
+ G_free(h);
+
if (type != 2) {
var_par->type = type;
var_par->nLag = nLag;
Modified: grass-addons/grass7/vector/v.kriging/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.kriging/local_proto.h 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/local_proto.h 2017-06-08 10:32:49 UTC (rev 71171)
@@ -38,161 +38,185 @@
struct opts
{
- struct Option *input, *output, *phase, *report, *crossvalid, *function_var_hz, *function_var_vert, *function_var_final, *function_var_final_vert, *form_file, *field, *intpl, *zcol, *var_dir_hz, *var_dir_vert, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert, *nugget_final, *nugget_final_vert, *sill_hz, *sill_vert, *sill_final, *sill_final_vert, *range_hz, *range_vert, *range_final, *range_final_vert;
+ struct Option *input, *output, *phase, *report, *crossvalid,
+ *function_var_hz, *function_var_vert, *function_var_final,
+ *function_var_final_vert, *form_file, *field, *intpl, *zcol, *trend_a,
+ *trend_b, *trend_c, *trend_d, *var_dir_hz, *var_dir_vert, *maxL,
+ *maxZ, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert,
+ *nugget_final, *nugget_final_vert, *sill_hz, *sill_vert, *sill_final,
+ *sill_final_vert, *range_hz, *range_vert, *range_final,
+ *range_final_vert;
};
struct flgs
{
- struct Flag *d23; /* 2D/3D interpolation */
- struct Flag *bivariate;
- struct Flag *univariate;
- struct Flag *detrend;
+ struct Flag *d23; /* 2D/3D interpolation */
+ struct Flag *bivariate;
+ struct Flag *univariate;
+ struct Flag *detrend;
};
struct select
{
- int n; // # of selected
- int out; // # of the others
- int total; // total number
- int *indices; // indices of selected
+ int n; // # of selected
+ int out; // # of the others
+ int total; // total number
+ int *indices; // indices of selected
};
-struct points // inputs
+struct points // inputs
{
- int n; // number of points
- double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
- struct RTree *R_tree; // spatial index
- struct RTree *Rtree_hz; // spatial index
- struct RTree *Rtree_vert; // spatial index
- double *r_min; // min coords
- double *r_max; // max coords
- double max_dist; // maximum distance
- double center[3]; // center
- double *invals; // values to be interpolated
- struct select in_reg; // points in region
- mat_struct *trend;
+ int n; // number of points
+ double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
+ struct RTree *R_tree; // spatial index
+ struct RTree *Rtree_hz; // spatial index
+ struct RTree *Rtree_vert; // spatial index
+ double *r_min; // min coords
+ double *r_max; // max coords
+ double max_dist; // maximum distance
+ double center[3]; // center
+ double *invals; // values to be interpolated
+ struct select in_reg; // points in region
+ mat_struct *trend;
};
struct bivar
{
- int vert;
- char *variogram;
+ int vert;
+ char *variogram;
- double *h;
- double max_dist;
- int nLag;
- double lag;
+ double *h;
+ double max_dist;
+ int nLag;
+ double lag;
- int function;
- double sill;
- double nugget;
- double h_range;
+ int function;
+ double sill;
+ double nugget;
+ double h_range;
};
struct parameters
{
- int function; // variogram function: lin, exp, spher, Gauss, bivar
- int type; // variogram type: hz / vert / aniso / bivar
- int const_val;
- double dir; // azimuth for variogram computing
- double td; // maximum azimuth
+ int function; // variogram function: lin, exp, spher, Gauss, bivar
+ int type; // variogram type: hz / vert / aniso / bivar
+ int const_val;
+ double dir; // azimuth for variogram computing
+ double td; // maximum azimuth
- double radius; // radius (squared maximum distance)
- double max_dist; // maximum distance - hz / aniso
- double max_dist_vert; // maximum distance - vert
+ double radius; // radius (squared maximum distance)
+ double max_dist; // maximum distance - hz / aniso
+ double max_dist_vert; // maximum distance - vert
- int nLag; // number of lags - hz / aniso
- double lag; // lag size
+ int nLag; // number of lags - hz / aniso
+ double lag; // lag size
- int nLag_vert; // number of lags - vert
- double lag_vert; // lag size
+ int nLag_vert; // number of lags - vert
+ double lag_vert; // lag size
- double *h; // lag distance from search point - hz / aniso
- double *vert; // lag distance from search point - vert
+ double *h; // lag distance from search point - hz / aniso
+ double *vert; // lag distance from search point - vert
- int gamma_n; // # of dissimilarities between input points
- mat_struct *gamma; // experimental variogram matrix
- double gamma_sum; // sum of gamma values
+ int gamma_n; // # of dissimilarities between input points
+ mat_struct *gamma; // experimental variogram matrix
+ double gamma_sum; // sum of gamma values
- double nugget;
- double sill;
- double part_sill;
- double h_range;
+ double nugget;
+ double sill;
+ double part_sill;
+ double h_range;
- struct bivar horizontal; // horizontal variogram properties
- struct bivar vertical; // vertical variogram properties
+ struct bivar horizontal; // horizontal variogram properties
+ struct bivar vertical; // vertical variogram properties
- mat_struct *A; // plan matrix
- mat_struct *T; // coefficients of theoretical variogram
- mat_struct *GM; // GM = theor_var(dist: input, output points)
+ mat_struct *A; // plan matrix
+ mat_struct *T; // coefficients of theoretical variogram
+ mat_struct *GM; // GM = theor_var(dist: input, output points)
- char *name; // name of input vector layer
- char term[12]; // output format - gnuplot terminal
- char ext[4]; // output format - file format
+ char *name; // name of input vector layer
+ char term[12]; // output format - gnuplot terminal
+ char ext[4]; // output format - file format
};
-struct var_par // parameters of experimental variogram
+struct var_par // parameters of experimental variogram
{
- struct parameters hz;
- struct parameters vert;
- struct parameters fin;
+ struct parameters hz;
+ struct parameters vert;
+ struct parameters fin;
+};
- //double Lmax; // maximal length in horizontal direction
- //double dzmax; // maximal elevation difference
- //double radius; // radius for kd-tree
+struct krig_pars // parameters of ordinary kriging
+{
+ int new;
+ int first;
+ int modified;
+ mat_struct *GM;
+ mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
+ mat_struct *rslt;
};
struct write
{
- int write2file; // write to file or not
- char *name; // filename
- FILE *fp;
- time_t now;
+ char *name; // filename
+ FILE *fp;
+ time_t now;
};
-struct int_par // Interpolation settings
+struct int_par // Interpolation settings
{
- int i3; // TRUE = 3D interpolation, FALSE = 2D interpolation (user sets by flag)
- char v3; // TRUE = 3D layer, FALSE = 2D layer (according to Vect_is_3d())
- int phase;
- int bivar;
- int univar;
- double aniso_ratio;
- int const_val;
- struct write report;
- struct write crossvalid;
+ int i3; // TRUE = 3D interpolation, FALSE = 2D interpolation (user sets by flag)
+ char v3; // TRUE = 3D layer, FALSE = 2D layer (according to Vect_is_3d())
+ int phase;
+ int bivar;
+ int univar;
+ double aniso_ratio;
+ int const_val;
+ struct write *report;
+ struct write *crossvalid;
};
-struct reg_par // Region settings -> output extent and resolution
+struct reg_par // Region settings -> output extent and resolution
{
- struct Cell_head reg_2d; // region for 2D interpolation
- RASTER3D_Region reg_3d; // region for 3D interpolation
- double west; // region.west
- double east; // region.east
- double north; // region.north
- double south; // region.south
- double bot; // region.bottom
- double top; // region.top
- double ew_res; // east-west resolution
- double ns_res; // north-south resolution
- double bt_res; // bottom-top resolution
- int nrows; // # of rows
- int ncols; // # of cols
- int ndeps; // # of deps
- int nrcd; // # of cells
+ struct Cell_head reg_2d; // region for 2D interpolation
+ RASTER3D_Region reg_3d; // region for 3D interpolation
+ double west; // region.west
+ double east; // region.east
+ double north; // region.north
+ double south; // region.south
+ double bot; // region.bottom
+ double top; // region.top
+ double ew_res; // east-west resolution
+ double ns_res; // north-south resolution
+ double bt_res; // bottom-top resolution
+ int nrows; // # of rows
+ int ncols; // # of cols
+ int ndeps; // # of deps
+ int nrcd; // # of cells
};
+struct trend
+{
+ double a;
+ double b;
+ double c;
+ double d;
+};
+
struct output
{
- DCELL *dcell;
- int fd_2d;
- RASTER3D_Map *fd_3d;
- char *name; // name of output 2D/3D raster
+ DCELL *dcell;
+ int fd_2d;
+ RASTER3D_Map *fd_3d;
+ char *name; // name of output 2D/3D raster
+ int add_trend;
+ double trend[4];
};
-double *get_col_values(struct Map_info *, struct int_par *, struct points *, int, const char *, int);
+double *get_col_values(struct Map_info *, struct int_par *, struct points *,
+ int, const char *);
void test_normality(int, double *, struct write *);
-void read_points(struct Map_info *, struct reg_par *, struct points *, struct int_par *, const char *, int, struct write *);
+void read_points(struct Map_info *, struct reg_par *, struct points *,
+ struct int_par *, const char *, int, struct write *);
double min(double *, struct points *);
double max(double *, struct points *);
@@ -202,31 +226,38 @@
struct ilist *find_n_NNs(int, int, struct points *, int);
double sum_NN(int, int, struct ilist *, struct points *);
-void correct_indices(struct ilist *, double *, struct points *, struct parameters *);
+void correct_indices(int, struct ilist *, double *, struct points *,
+ struct parameters *);
int cmpVals(const void *, const void *);
void coord_diff(int, int, double *, double *);
double distance_diff(double *);
double radius_hz_diff(double *);
double zenith_angle(double *);
void triple(double, double, double, double *);
-double lag_size(int, struct points *, struct parameters *, struct write *);
+double lag_size(int, int, struct points *, struct parameters *,
+ struct write *);
int lag_number(double, double *);
void optimize(double *, int *, double);
-void variogram_restricts(struct int_par *, struct points *, struct parameters *);
+void variogram_restricts(struct int_par *, struct points *,
+ struct parameters *);
void geometric_anisotropy(struct int_par *, struct points *);
-double find_intersect_x(double *, double *, double *, double *, struct write *);
-double find_intersect_y(double *, double *, double *, double *, double , struct write *);
+double find_intersect_x(double *, double *, double *, double *,
+ struct write *);
+double find_intersect_y(double *, double *, double *, double *, double,
+ struct write *);
mat_struct *LSM(mat_struct *, mat_struct *);
-mat_struct *nonlin_LMS(int , double *, double *);
+mat_struct *nonlin_LMS(int, double *, double *);
void E_variogram(int, struct int_par *, struct points *, struct var_par *);
-void T_variogram(int, int, struct opts, struct parameters *, struct write *);
-void ordinary_kriging(struct int_par *, struct reg_par *, struct points *, struct var_par *, struct output *);
+void T_variogram(int, struct opts, struct parameters *, struct int_par *);
+void ordinary_kriging(struct int_par *, struct reg_par *, struct points *,
+ struct var_par *, struct output *);
void LMS_variogram(struct parameters *, struct write *);
double bivar_sill(int, mat_struct *);
void sill(struct parameters *);
-void sill_compare(struct int_par *, struct flgs *, struct var_par *, struct points *);
+void sill_compare(struct int_par *, struct flgs *, struct var_par *,
+ struct points *);
int set_function(char *, struct write *);
double RBF(double *);
double linear(double, double, double);
@@ -236,29 +267,53 @@
double variogram_fction(struct parameters *, double *);
void set_gnuplot(char *, struct parameters *);
void plot_experimental_variogram(struct int_par *, struct parameters *);
-void plot_var(int, int, struct parameters *);
+void plot_var(struct int_par *, int, struct parameters *);
+void new_vertical(int *, int);
+
void variogram_type(int, char *);
void write2file_basics(struct int_par *, struct opts *);
void write2file_vector(struct int_par *, struct points *);
void write2file_values(struct write *, const char *);
void write2file_varSetsIntro(int, struct write *);
void write2file_varSets(struct write *, struct parameters *);
-void write2file_variogram_E(struct int_par *, struct parameters *, mat_struct *);
+void write2file_variogram_E(struct int_par *, struct parameters *,
+ mat_struct *);
void write2file_variogram_T(struct write *);
void write_temporary2file(struct int_par *, struct parameters *);
+void report_error(struct write *);
void read_tmp_vals(const char *, struct parameters *, struct int_par *);
-
-mat_struct *set_up_G(struct points *, struct parameters *, struct write *);
-mat_struct *set_up_g0(struct int_par *, struct points *, struct ilist *, double *, struct parameters *);
+
+void set_up_G(struct points *, struct parameters *, struct write *,
+ struct krig_pars *);
+mat_struct *set_up_g0(struct int_par *, struct points *, struct ilist *,
+ double *, struct parameters *);
mat_struct *submatrix(struct ilist *, mat_struct *, struct write *);
-double result(struct points *, struct ilist *, mat_struct *);
+double result(struct points *, struct ilist *, mat_struct *);
-void crossvalidation(struct int_par *, struct points *, struct parameters *);
-void cell_centre(unsigned int, unsigned int, unsigned int, struct int_par *, struct reg_par *, double *, struct parameters *);
+double find_center(double, double);
+void adjacent_cell(int, double *, struct reg_par *, double *);
+void crossvalidation(struct int_par *, struct points *, struct parameters *,
+ struct reg_par *);
+void cell_centre(unsigned int, unsigned int, unsigned int, struct int_par *,
+ struct reg_par *, double *, struct parameters *);
+struct ilist *list_NN(struct int_par *, double *, struct points *, double,
+ double);
+int compare_NN(struct ilist *, struct ilist *, int);
+void make_subsamples(struct int_par *, struct ilist *, double *, int, int,
+ struct points *, struct parameters *,
+ struct krig_pars *);
+double interpolate(struct int_par *, struct ilist *, double *,
+ struct points *, struct parameters *, struct krig_pars *);
+double trend(double *, struct output *, int, struct int_par *);
+int new_sample(struct int_par *, struct ilist *, struct ilist *,
+ struct points *, int, int, int, double *, double, double,
+ struct reg_par *, struct parameters *, struct krig_pars *,
+ int *);
void get_region_pars(struct int_par *, struct reg_par *);
void open_layer(struct int_par *, struct reg_par *, struct output *);
-int write2layer(struct int_par *, struct reg_par *, struct output *, unsigned int, unsigned int, unsigned int, double);
+void write2layer(struct int_par *, struct reg_par *, struct output *,
+ mat_struct *);
double get_quantile(int);
void get_slot_counts(int, double *);
Modified: grass-addons/grass7/vector/v.kriging/main.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/main.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/main.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -41,6 +41,13 @@
// Geostatistical parameters
struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
+
+ xD.report = (struct write *)G_malloc(sizeof(struct write));
+ struct write *report = xD.report;
+
+ xD.crossvalid = (struct write *)G_malloc(sizeof(struct write));
+ struct write *crossvalid = xD.crossvalid;
+
struct var_par var_pars; // Variogram (experimental and theoretical)
// Outputs
@@ -140,8 +147,40 @@
flg.detrend = G_define_flag();
flg.detrend->key = 't';
flg.detrend->description = _("Eliminate trend if variogram is parabolic");
- flg.detrend->guisection = _("Initial");
+ flg.detrend->guisection = _("Final");
+ opt.trend_a = G_define_option();
+ opt.trend_a->key = "atrend";
+ opt.trend_a->type = TYPE_DOUBLE;
+ opt.trend_a->required = NO;
+ opt.trend_a->answer = "0.0";
+ opt.trend_a->description = _("Trend: f(x,y,z) = a*x + b*y + c*z + d");
+ opt.trend_a->guisection = _("Final");
+
+ opt.trend_b = G_define_option();
+ opt.trend_b->key = "btrend";
+ opt.trend_b->type = TYPE_DOUBLE;
+ opt.trend_b->required = NO;
+ opt.trend_b->answer = "0.0";
+ opt.trend_b->description = _("Trend: f(x,y,z) = a*x + b*y + c*z + d");
+ opt.trend_b->guisection = _("Final");
+
+ opt.trend_c = G_define_option();
+ opt.trend_c->key = "ctrend";
+ opt.trend_c->type = TYPE_DOUBLE;
+ opt.trend_c->required = NO;
+ opt.trend_c->answer = "0.0";
+ opt.trend_c->description = _("Trend: f(x,y,z) = a*x + b*y + c*z + d");
+ opt.trend_c->guisection = _("Final");
+
+ opt.trend_d = G_define_option();
+ opt.trend_d->key = "dtrend";
+ opt.trend_d->type = TYPE_DOUBLE;
+ opt.trend_d->required = NO;
+ opt.trend_d->answer = "0.0";
+ opt.trend_d->description = _("Trend: f(x,y,z) = a*x + b*y + c*z + d");
+ opt.trend_d->guisection = _("Final");
+
opt.form_file = G_define_option(); // Variogram plot - various output formats
opt.form_file->key = "fileformat";
opt.form_file->options = "cdr,dxf,eps,tex,pdf,png,svg";
@@ -180,6 +219,21 @@
_("Zenith angle of variogram computing (isotrophic)");
opt.var_dir_vert->guisection = _("Initial");
+ opt.maxL = G_define_option();
+ opt.maxL->key = "lmax";
+ opt.maxL->type = TYPE_DOUBLE;
+ opt.maxL->required = NO;
+ opt.maxL->description = _("Maximum distance in horizontal direction");
+ opt.maxL->guisection = _("Initial");
+
+ opt.maxZ = G_define_option();
+ opt.maxZ->key = "vmax";
+ opt.maxZ->type = TYPE_DOUBLE;
+ opt.maxZ->required = NO;
+ opt.maxZ->description =
+ _("Maximum distance in horizontal direction (only for 3D variogram)");
+ opt.maxZ->guisection = _("Initial");
+
opt.nL = G_define_option();
opt.nL->key = "lpieces";
opt.nL->type = TYPE_INTEGER;
@@ -316,39 +370,30 @@
}
// Open report file if desired
- if (opt.report->answer) {
- xD.report.write2file = TRUE;
- xD.report.name = opt.report->answer;
+ if (xD.phase == 0) {
+ if (opt.report->answer) {
+ report->name = opt.report->answer;
- // initial phase: check if the file exists
- if (xD.phase == 0 && access(xD.report.name, F_OK) != -1) {
- G_fatal_error(_("Report file exists; please set up different name..."));
- }
+ // initial phase: check if the file exists
+ if (access(report->name, F_OK) != -1) {
+ G_fatal_error(_("Report file exists; please set up different name..."));
+ }
- // middle / final phase: check if file does not exist
- if (xD.phase != 0 && access(xD.report.name, F_OK) != 0) {
- G_fatal_error(_("Report file does not exist; please check the name or repeat initial phase..."));
+ report->fp = fopen(report->name, "a");
+ time(&report->now);
+ fprintf(report->fp, "v.kriging started on %s\n\n",
+ ctime(&report->now));
+ G_message(_("Report is being written to %s..."), report->name);
}
-
- xD.report.fp = fopen(xD.report.name, "a");
- time(&xD.report.now);
- fprintf(xD.report.fp, "v.kriging started on %s\n\n",
- ctime(&xD.report.now));
- G_message(_("Report is being written to %s..."), xD.report.name);
+ else {
+ G_warning(_("The name of report file missing..."));
+ }
}
- else {
- xD.report.write2file = FALSE;
- G_warning(_("The name of report file missing..."));
- }
- if (opt.crossvalid->answer) {
- xD.crossvalid.write2file = TRUE;
- xD.crossvalid.name = opt.crossvalid->answer;
- xD.crossvalid.fp = fopen(xD.crossvalid.name, "w");
+ if (xD.phase == 2 && opt.crossvalid->answer) {
+ crossvalid->name = opt.crossvalid->answer;
+ crossvalid->fp = fopen(crossvalid->name, "w");
}
- else {
- xD.crossvalid.write2file = FALSE;
- }
var_pars.hz.name = var_pars.vert.name = var_pars.fin.name = opt.input->answer; // set name of variogram
var_pars.hz.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
@@ -357,14 +402,43 @@
var_pars.hz.td = DEG2RAD(atof(opt.td_hz->answer)); // Angle of variogram processing
- if (opt.nL->answer) { // Test if nL have been set up (optional)
- var_pars.hz.nLag = atoi(opt.nL->answer);
- if (var_pars.hz.nLag < 1) { // Invalid value
- G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
- var_pars.hz.nLag = 20;
+ if (xD.phase == 0) {
+ if (opt.maxL->answer) { // Test if nL have been set up (optional)
+ var_pars.hz.max_dist = atoi(opt.maxL->answer);
+ if (var_pars.hz.max_dist < 0.) { // Invalid value
+ G_fatal_error(_("Maximum distance must be positive..."));
+ }
}
+ else {
+ var_pars.hz.max_dist = -1.;
+ }
}
+ if (xD.phase == 0) {
+ if (opt.maxZ->answer) { // Test if nL have been set up (optional)
+ var_pars.vert.max_dist = atoi(opt.maxZ->answer);
+ if (var_pars.vert.max_dist < 0.) { // Invalid value
+ G_fatal_error(_("Maximum distance must be positive..."));
+ }
+ }
+ else {
+ var_pars.vert.max_dist = -1.;
+ }
+ }
+
+ if (xD.phase == 0) {
+ if (opt.nL->answer) { // Test if nL have been set up (optional)
+ var_pars.hz.nLag = atoi(opt.nL->answer);
+ if (var_pars.hz.nLag < 1) { // Invalid value
+ G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
+ var_pars.hz.nLag = 20;
+ }
+ }
+ else {
+ var_pars.hz.nLag = -1.;
+ }
+ }
+
if (opt.form_file->answer) { // Plotting variogram
set_gnuplot(opt.form_file->answer, &var_pars.hz);
set_gnuplot(opt.form_file->answer, &var_pars.vert);
@@ -379,9 +453,9 @@
xD.bivar = flg.bivariate->answer == TRUE ? TRUE : FALSE;
xD.univar = flg.univariate->answer == TRUE ? TRUE : FALSE;
if (xD.bivar == TRUE && xD.univar == TRUE) {
- if (xD.report.write2file == TRUE) {
- fclose(xD.report.fp);
- remove(xD.report.name);
+ if (report->name) {
+ fclose(report->fp);
+ remove(report->name);
}
G_fatal_error(_("You should mark either univariate, or bivariate variogram, not both of them..."));
} // error
@@ -390,9 +464,9 @@
Vect_set_open_level(2); // Open input vector map
if (0 > Vect_open_old2(&map, opt.input->answer, "", opt.field->answer)) {
- if (xD.report.write2file == TRUE) {
- fclose(xD.report.fp);
- remove(xD.report.name);
+ if (xD.phase == 0 && report->name) {
+ fclose(report->fp);
+ remove(report->name);
}
G_fatal_error(_("Unable to open vector map <%s>"), opt.input->answer);
} // error
@@ -414,46 +488,51 @@
if (xD.i3 == TRUE) { // 3D interpolation:
if (xD.v3 == FALSE) { // 2D input:
if (!opt.zcol->answer) { // zcolumn not available:
- if (xD.report.write2file == TRUE) { // close report file
- fprintf(xD.report.fp,
+ if (report->name) { // close report file
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("To process 3D interpolation based on 2D input, please set attribute column containing z coordinates or switch to 2D interpolation."));
}
} // end if zcol == NULL
// 3D or 2,5D input
- if (opt.nZ->answer) { // Test if nZ have been set up (optional)
- if (var_pars.vert.nLag < 1) { // Invalid value
- G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
- }
+ if (xD.phase == 0) {
+ if (opt.nZ->answer) { // Test if nZ have been set up (optional)
+ if (var_pars.vert.nLag < 1) { // Invalid value
+ G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
+ }
+ else {
+ var_pars.vert.nLag = atof(opt.nZ->answer);
+ }
+ } // end if nZ->answer
else {
- var_pars.vert.nLag = atof(opt.nZ->answer);
+ var_pars.vert.nLag = -1.;
}
- } // end if nZ->answer
+ }
} // end if 3D interpolation
else { // 2D interpolation:
var_pars.vert.nLag = -1; // abs will be used in next steps
if (xD.v3 == TRUE) {
- if (xD.report.write2file == TRUE) { // close report file
- fprintf(xD.report.fp,
+ if (report->name) { // close report file
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("2D interpolation based on 3D input has been temporarily disabled. Please select another option..."));
}
}
field = Vect_get_field_number(&map, opt.field->answer);
- if (xD.report.write2file == TRUE)
+ if (xD.phase == 0 && report->name)
write2file_basics(&xD, &opt);
/* ---------------------------------------------------------- */
/* Get... */
get_region_pars(&xD, ®); // ... region parameters
- read_points(&map, ®, &pnts, &xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
- pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
+ read_points(&map, ®, &pnts, &xD, opt.zcol->answer, field, report); // ... coordinates of points
+ pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer); // ... values for interpolation
/* Estimate 2D/3D variogram */
switch (xD.phase) {
@@ -474,7 +553,7 @@
}
E_variogram(0, &xD, &pnts, &var_pars); // horizontal variogram (for both 2D and 3D interpolation)
- if (xD.report.write2file == FALSE) {
+ if (!report->name) {
G_message(_("\nExperimental variogram of your data has been computed. To continue interpolation performance, please repeat initial phase with non-empty <report> parameter..."));
}
@@ -490,11 +569,13 @@
goto end;
case 1:
- read_tmp_vals("variogram_hz_tmp.txt", &var_pars.hz, &xD); // read properties of horizontal variogram from temp file
- read_tmp_vals("variogram_vert_tmp.txt", &var_pars.vert, &xD); // read properties of vertical variogram from temp file
+ // read properties of variogram from temporary file
+ read_tmp_vals("variogram_hz_tmp.txt", &var_pars.hz, &xD); // horizontal
+ read_tmp_vals("variogram_vert_tmp.txt", &var_pars.vert, &xD); // vertical
- T_variogram(0, TRUE, opt, &var_pars.hz, &xD.report); // compute theoretical variogram - hz
- T_variogram(1, TRUE, opt, &var_pars.vert, &xD.report); // compute theoretical variogram - vert
+ // compute theoretical variogram
+ T_variogram(0, opt, &var_pars.hz, &xD); // horizontal
+ T_variogram(1, opt, &var_pars.vert, &xD); // vertical
/* compare range of hz and vert variogram:
- if the difference is greater than 5% -> bivariate variogram
@@ -513,14 +594,28 @@
case 2: // final phase:
// Module should crash if:
if (!opt.output->answer) { // output name not available:
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up name of output layer..."));
}
+ if (opt.trend_a->answer) {
+ out.trend[0] = atof(opt.trend_a->answer);
+ out.trend[1] = atof(opt.trend_b->answer);
+ out.trend[2] = atof(opt.trend_c->answer);
+ out.trend[3] = atof(opt.trend_d->answer);
+ }
+
+ /*// read trend
+ out.add_trend = flg.detrend->answer;
+ if (out.add_trend == TRUE) {
+ out.trend = G_matrix_init();
+ out.trend->vals = get_col_values(&map, &xD, &pnts, field, opt.trend->answer); // trend
+ } */
+
// read properties of final variogram from temp file
if (xD.i3 == TRUE) { // 3D kriging:
read_tmp_vals("variogram_final_tmp.txt", &var_pars.fin, &xD);
@@ -537,10 +632,10 @@
opt.function_var_final_vert->answer) &&
(opt.function_var_final->answer ||
opt.function_var_final_vert->answer)) {
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("If you wish to specify components of bivariate variogram please set up function type for both of them..."));
}
@@ -548,76 +643,86 @@
// if both of the function types are set up:
if (opt.function_var_final->answer &&
opt.function_var_final_vert->answer) {
- if (!opt.nugget_final->answer) { // horizontal nugget effect should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ // horizontal nugget effect should not be missing
+ if (!opt.nugget_final->answer) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up nugget effect of horizontal component of bivariate variogram..."));
}
- if (!opt.nugget_final_vert->answer) { // vertical nugget effect should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ // vertical nugget effect should not be missing
+ if (!opt.nugget_final_vert->answer) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up nugget effect of vertical component of bivariate variogram..."));
}
- if (!opt.range_final->answer) { // horizontal range should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ // horizontal range should not be missing
+ if (!opt.range_final->answer) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up range of horizontal component of bivariate variogram..."));
}
- if (!opt.range_final_vert->answer) { // vertical range should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ // vertical range should not be missing
+ if (!opt.range_final_vert->answer) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up range of vertical component of bivariate variogram..."));
}
}
}
- else { // univariate variogram
+ // univariate variogram
+ else {
if (xD.i3 == TRUE && (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0)) { // anisotropic 3D:
if (opt.function_var_final_vert->answer || atof(opt.nugget_final_vert->answer) != 0. || opt.range_final_vert->answer) { // vertical settings available:
G_fatal_error(_("Not necessary to set up vertical components properties. Anisotropic variogram will be used..."));
} // end if vert settings available
} // end if 3D
- if (!opt.function_var_final->answer) { // missing function
+ // missing function
+ if (!opt.function_var_final->answer) {
opt.function_var_final->answer = "linear";
G_message(_("Linear variogram will be computed..."));
}
- if (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0) { // nonlinear and not parabolic:
+ // nonlinear and not parabolic:
+ if (strcmp(opt.function_var_final->answer, "linear") != 0 &&
+ strcmp(opt.function_var_final->answer, "parabolic") != 0) {
if (!opt.nugget_final->answer) { // missing horizontal nugget effect
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up nugget effect..."));
} // end if nugget missing
- if (!opt.range_final->answer) { // missing horizontal range:
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp,
+ // missing horizontal range:
+ if (!opt.range_final->answer) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
"Error (see standard output). Process killed...");
- fclose(xD.report.fp);
+ fclose(report->fp);
}
G_fatal_error(_("Please set up range..."));
} // end if range missing
- if (opt.sill_final->answer) { // missing horizontal range:
+ // missing horizontal range:
+ if (opt.sill_final->answer) {
var_pars.fin.sill = atof(opt.sill_final->answer);
} // end if sill has been changed by the user
} // end nonlinear and not parabolic variogram
@@ -625,11 +730,13 @@
out.name = opt.output->answer; // Output layer name
- if (var_pars.fin.type == 3) { // if variogram is anisotropic:
+ // if variogram is anisotropic:
+ if (var_pars.fin.type == 3) {
geometric_anisotropy(&xD, &pnts); // exaggerate z coord and rebuild the spatial index
}
- T_variogram(var_pars.fin.type, xD.i3, opt, &var_pars.fin, &xD.report); // compute theoretical variogram
+ // compute theoretical variogram
+ T_variogram(var_pars.fin.type, opt, &var_pars.fin, &xD);
break;
}
Modified: grass-addons/grass7/vector/v.kriging/quantile.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/quantile.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/quantile.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -178,7 +178,7 @@
*quantile = v;
prev_v = v;
- if (report != NULL && report->write2file == TRUE)
+ if (report != NULL && report->name)
fprintf(report->fp, "%f:%f:%f\n", prev_v, v, quants[quant]);
}
}
Modified: grass-addons/grass7/vector/v.kriging/spatial_index.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/spatial_index.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/spatial_index.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -75,8 +75,8 @@
NN_sum = list->n_values;
if (NN_sum < 2) {
- G_warning(_("Point \"x=%f y=%f z=%f\" has less than 2 neighbours in its closest surrounding. The perimeter of the surrounding will be increased to include more neighbouring points"),
- *r, *(r + 1), *(r + 2));
+ /*G_fatal_error(_("Point \"x=%f y=%f z=%f\" has less than 2 neighbours in its closest surrounding. The perimeter of the surrounding will be increased to include more neighbouring points. If this warning appears regularly, please check your region settings."),
+ *r, *(r + 1), *(r + 2));*/
max_dist += dist_step;
max_dist_vert += dist_step_vert;
}
@@ -112,7 +112,7 @@
}
// check spherical (circular) surrounding
- list = find_NNs_within(dim, search, pnts, sqrt(2.) * max_dist, -1);
+ list = find_NNs_within(dim, search, pnts, max_dist, -1);
return list;
}
Modified: grass-addons/grass7/vector/v.kriging/utils.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/utils.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -16,51 +16,6 @@
return 0;
}
-void correct_indices(struct ilist *list, double *r0, struct points *pnts,
- struct parameters *var_pars)
-{
- int i;
- int n = list->n_values;
- int *vals = list->value;
- double dx, dy, dz;
- int type = var_pars->type;
-
- int n_new = 0; // # of effective indices (not identical points, nor too far)
- int *newvals, *save; // effective indices
-
- newvals = (int *)G_malloc(n * sizeof(int));
- save = &newvals[0];
-
- double *r, sqDist_i;
-
- for (i = 0; i < n; i++) {
- *vals -= 1;
- r = &pnts->r[3 * *vals];
-
- dx = *r - *r0;
- dy = *(r + 1) - *(r0 + 1);
- dz = *(r + 2) - *(r0 + 2);
- sqDist_i = SQUARE(dx) + SQUARE(dy) + SQUARE(dz);
-
- if (type != 2 && (n == 1 || (sqDist_i != 0.))) {
- *save = *vals;
- n_new++;
- save++;
- }
- vals++;
- }
-
- if (type != 2 && n_new < n) {
- list->n_values = n_new;
- list->value = (int *)G_realloc(list->value, n_new * sizeof(int));
- memcpy(list->value, newvals, n_new * sizeof(int));
- }
-
- G_free(newvals);
-
- return;
-}
-
// make coordinate triples xyz
void triple(double x, double y, double z, double *triple)
{
@@ -97,7 +52,7 @@
}
// compute horizontal radius from coordinate differences
-double radius_hz_diff(double *dr)
+double radius_hz_diff(double dr[3])
{
double rds;
@@ -106,6 +61,25 @@
return rds;
}
+double squared_distance(int direction, double *dr)
+{
+ double d;
+
+ switch (direction) {
+ case 12: // horizontal
+ d = SQUARE(dr[0]) + SQUARE(dr[1]);
+ break;
+ case 3: // vertical
+ d = SQUARE(dr[2]);
+ break;
+ case 0: // all
+ d = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
+ break;
+ }
+
+ return d;
+}
+
// compute zenith angle
double zenith_angle(double dr[3])
{
@@ -116,8 +90,66 @@
return zenith;
}
+void correct_indices(int direction, struct ilist *list, double *r0,
+ struct points *pnts, struct parameters *var_pars)
+{
+ int i;
+ int n = list->n_values;
+ int *vals = list->value;
+ double dr[3];
+ int type = var_pars->type;
+ double max_dist = var_pars->max_dist;
+
+ int skip = FALSE; // do not skip relevant points
+ int n_new = 0; // # of effective indices (not identical points, nor too far)
+ int *newvals, *save; // effective indices
+
+ // values of the new list (if necesary)
+ newvals = (int *)G_malloc(n * sizeof(int));
+ save = &newvals[0];
+
+ double *r, sqDist_i;
+
+ for (i = 0; i < n; i++) {
+ *vals -= 1; // decrease value by 1
+ r = &pnts->r[3 * *vals]; // find point coordinates
+
+ // coordinate differences
+ *dr = *r - *r0;
+ *(dr + 1) = *(r + 1) - *(r0 + 1);
+ *(dr + 2) = *(r + 2) - *(r0 + 2);
+
+ // compute squared distance
+ sqDist_i = squared_distance(direction, dr);
+
+ // compare with maximum
+ skip = sqDist_i <= max_dist ? FALSE : TRUE;
+
+ // if:
+ // - univariate variogram with only one non-zero value or
+ // - the distance (circular) does not exceed max_dist
+ if (type != 2 && (n == 1 || (sqDist_i != 0.)) || skip == FALSE) {
+ *save = *vals;
+ n_new++;
+ save++;
+ skip = FALSE;
+ }
+ vals++;
+ }
+
+ if (type != 2 && n_new < n) {
+ list->n_values = n_new;
+ list->value = (int *)G_realloc(list->value, n_new * sizeof(int));
+ memcpy(list->value, newvals, n_new * sizeof(int));
+ }
+
+ G_free(newvals);
+
+ return;
+}
+
// compute size of the lag
-double lag_size(int direction, struct points *pnts,
+double lag_size(int i3, int direction, struct points *pnts,
struct parameters *var_pars, struct write *report)
{
// local variables
@@ -168,7 +200,7 @@
coord_diff(i, (*j_vals - 1), pnts->r, dr); // compute coordinate differences
switch (direction) { // and count those which equal to zero
case 12: // horizontal
- if (dr[0] == 0. && dr[1] == 0.) {
+ if (i3 == FALSE && (dr[0] == 0. && dr[1] == 0.)) {
add_ident++;
}
break;
@@ -188,16 +220,12 @@
} //end if: n_vals > 0
else {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Error: There are no nearest neighbours of point %d..."),
i);
} // end error
- perc5 = (int)round(0.05 * n) + add_ident; // set up 5% increased by number of identical points
+ perc5 = (int)ceil(0.05 * n) + add_ident; // set up 5% increased by number of identical points
coord_diff(i, perc5, pnts->r, dr); // coordinate differences of search point and 5%th NN
switch (direction) { // compute squared distance for particular direction
@@ -216,13 +244,12 @@
d_min = MIN(d_min, dist2);
}
+ G_free_ilist(list); // free list memory
r += 3; // go to next search point
} // end i for loop: distance of i-th search pt and 5%-th NN
lagNN = sqrt(d_min);
- G_free_ilist(list); // free list memory
-
return lagNN;
}
@@ -296,7 +323,7 @@
void variogram_restricts(struct int_par *xD, struct points *pnts,
struct parameters *var_pars)
{
- struct write *report = &xD->report;
+ struct write *report = xD->report;
double *min, *max; // extend
double dr[3]; // coordinate differences
@@ -308,6 +335,7 @@
// Find extent
G_message(_("Computing %s variogram properties..."), type);
+
min = &pnts->r_min[0]; // set up pointer to minimum xyz
max = &pnts->r_max[0]; // set up pointer to maximum xyz
@@ -317,16 +345,26 @@
switch (var_pars->type) {
case 1: // vertical variogram
- var_pars->max_dist = dr[2]; // zmax - zmin
+ if (var_pars->max_dist == -1.) {
+ var_pars->max_dist = dr[2]; // zmax - zmin
+ }
var_pars->radius = SQUARE(var_pars->max_dist); // anisotropic distance (todo: try also dz)
break;
default:
- var_pars->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
- var_pars->max_dist = sqrt(var_pars->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
+ if (var_pars->max_dist == -1.) {
+ var_pars->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
+ var_pars->max_dist = sqrt(var_pars->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
+ }
+ else {
+ var_pars->radius = SQUARE(var_pars->max_dist);
+ }
break;
}
+ if (var_pars->max_dist <= 0.) {
+ G_fatal_error(_("Maximum distance must be greater than 0."));
+ }
- if (report->write2file == TRUE) { // report name available:
+ if (report->name) { // report name available:
write2file_varSetsIntro(var_pars->type, report); // describe properties
}
@@ -345,27 +383,42 @@
if (var_pars->type == 2) { // bivariate variogram
// horizontal direction:
- var_pars->lag = lag_size(12, pnts, var_pars, report); // lag distance
- var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // number of lags
- optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
- var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ if (var_pars->nLag == -1) {
+ var_pars->lag = lag_size(xD->i3, 12, pnts, var_pars, report); // lag distance
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // number of lags
+ optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ }
+ else {
+ var_pars->lag = var_pars->max_dist / var_pars->nLag; // lag size
+ }
// vertical direction
- var_pars->lag_vert = lag_size(3, pnts, var_pars, report); // lag distance
- var_pars->nLag_vert = lag_number(var_pars->lag_vert, &var_pars->max_dist_vert); // # of lags
- optimize(&var_pars->lag_vert, &var_pars->nLag_vert,
- var_pars->max_dist_vert);
- var_pars->max_dist_vert = var_pars->nLag_vert * var_pars->lag_vert; // max distance
+ if (var_pars->nLag_vert == -1) {
+ var_pars->lag_vert = lag_size(xD->i3, 3, pnts, var_pars, report); // lag distance
+ var_pars->nLag_vert = lag_number(var_pars->lag_vert, &var_pars->max_dist_vert); // # of lags
+ optimize(&var_pars->lag_vert, &var_pars->nLag_vert,
+ var_pars->max_dist_vert);
+ var_pars->max_dist_vert = var_pars->nLag_vert * var_pars->lag_vert; // max distance
+ }
+ else {
+ var_pars->lag_vert = var_pars->max_dist_vert / var_pars->nLag_vert; // lag size
+ }
}
else { // univariate variograms (hz / vert / aniso)
- var_pars->lag = lag_size(dimension, pnts, var_pars, report); // lag size
- var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // # of lags
- optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
- var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ if (var_pars->nLag == -1) {
+ var_pars->lag = lag_size(xD->i3, dimension, pnts, var_pars, report); // lag size
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // # of lags
+ optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ }
+ else {
+ var_pars->lag = var_pars->max_dist / var_pars->nLag; // lag size
+ }
}
- if (report->write2file == TRUE) { // report name available:
+ if (report->name) { // report name available:
write2file_varSets(report, var_pars); // describe properties
}
}
@@ -376,10 +429,10 @@
double ratio = xD->aniso_ratio;
if (ratio <= 0.) { // ratio is negative or zero:
- if (xD->report.name) { // close report file
- fprintf(xD->report.fp,
+ if (xD->report->name) { // close report file
+ fprintf(xD->report->fp,
"Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ fclose(xD->report->fp);
}
G_fatal_error(_("Anisotropy ratio must be greater than zero..."));
}
@@ -446,11 +499,7 @@
function = 5;
}
else {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Set up correct name of variogram function..."));
}
@@ -508,6 +557,7 @@
gp = fopen("dataE.dat", "w"); // open file to write experimental variogram (initial phase)
if (access("dataE.dat", W_OK) < 0) {
+ report_error(xD->report);
G_fatal_error(_("Something went wrong opening tmp file..."));
}
@@ -620,14 +670,16 @@
}
fclose(gp);
- if (xD->report.write2file == FALSE) {
+ if (!xD->report->name) {
remove("dataE.dat");
}
}
// plot experimental and theoretical variogram
-void plot_var(int i3, int bivar, struct parameters *var_pars)
+void plot_var(struct int_par *xD, int bivar, struct parameters *var_pars)
{
+ int i3 = xD->i3;
+ struct write *report = xD->report;
int function;
double nugget, nugget_h, nugget_v;
double sill;
@@ -666,9 +718,11 @@
if (type != 2) { // univariate variogram
gp = fopen("dataE.dat", "w"); // open file to write experimental variogram
if (access("dataE.dat", W_OK) < 0) {
+ report_error(report);
G_fatal_error(_("Something went wrong opening tmp file..."));
}
if (gp == NULL) {
+ report_error(report);
G_fatal_error(_("Error writing file"));
}
@@ -688,17 +742,20 @@
else { // bivariate variogram
gp = fopen("dataE.dat", "r"); // open file to read experimental variogram
if (access("dataE.dat", F_OK) < 0) {
+ report_error(report);
G_fatal_error(_("You have probably deleted dataE.dat - process middle phase again, please."));
}
}
if (fclose(gp) != 0) {
+ report_error(report);
G_fatal_error(_("Error closing file..."));
}
/* Theoretical variogram */
gp = fopen("dataT.dat", "w"); // open file to write theoretical variogram
if (access("dataT.dat", W_OK) < 0) {
+ report_error(report);
G_fatal_error(_("Something went wrong opening tmp file..."));
}
@@ -789,6 +846,7 @@
}
if (fclose(gp) == EOF) {
+ report_error(report);
G_fatal_error(_("Error closing file..."));
}
@@ -925,3 +983,16 @@
remove("dataE.dat");
remove("dataT.dat");
}
+
+void new_vertical(int *row0, int n)
+{
+ int i;
+ int *value;
+
+ value = &row0[0];
+
+ for (i = 0; i < n; i++) {
+ *value = 0;
+ value++;
+ }
+}
Modified: grass-addons/grass7/vector/v.kriging/utils_kriging.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_kriging.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/utils_kriging.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -45,7 +45,6 @@
switch (var_par->function) { // function of theoretical variogram
case 0: // linear variogram
G_matrix_set_element(var_par->A, nr, 0, *h);
- //G_matrix_set_element(var_par->A, nr, 1, 1.);
break;
case 1: // parabolic variogram
G_matrix_set_element(var_par->A, nr, 0, SQUARE(*h));
@@ -73,7 +72,7 @@
// Test of theoretical variogram estimation
if (var_par->T->vals == NULL) { // NULL values of theoretical variogram
- if (report->write2file == TRUE) { // close report file
+ if (report->name) { // close report file
fprintf(report->fp,
"Error (see standard output). Process killed...");
fclose(report->fp);
@@ -89,11 +88,17 @@
// coefficients of theoretical variogram (linear)
if (report->name) {
- fprintf(report->fp, "Parameters of bivariate variogram:\n");
- fprintf(report->fp, "Function: linear\n\n");
- fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + %f\n",
- var_par->T->vals[0], var_par->T->vals[1],
- var_par->T->vals[2]);
+ if (var_par->function == 5) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + %f\n",
+ var_par->T->vals[0], var_par->T->vals[1],
+ var_par->T->vals[2]);
+ }
+ else {
+ fprintf(report->fp, "Parameters of univariate variogram:\n");
+ fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + %f\n",
+ var_par->T->vals[0], var_par->T->vals[1]);
+ }
} // end if: report
}
@@ -192,6 +197,7 @@
else if (xD->univar == TRUE || (!flg->bivariate->answer && diff_sill <= diff_sill_05)) { // geometric anisotropy
var_par->fin.type = 3; // variogram type: anisotropic
var_par->fin.max_dist = var_par->hz.max_dist; // maximum distance: hz
+ var_par->fin.nLag = var_par->hz.nLag;
var_par->fin.td = var_par->hz.td; // azimuth: hz
xD->aniso_ratio = var_par->hz.h_range / var_par->vert.h_range; // anisotropic ratio
geometric_anisotropy(xD, pnts); // exaggerate z coords and build a new spatial index
@@ -315,8 +321,8 @@
}
// set up elements of G matrix
-mat_struct *set_up_G(struct points *pnts, struct parameters *var_par,
- struct write *report)
+void set_up_G(struct points *pnts, struct parameters *var_par,
+ struct write *report, struct krig_pars *krig)
{
// Local variables
int n = pnts->n; // # of input points
@@ -385,12 +391,13 @@
} // end i loop
free(dr);
- return GM;
+
+ krig->GM = G_matrix_copy(GM);
}
// make G submatrix for rellevant points
-mat_struct *submatrix(struct ilist * index, mat_struct * GM_all,
- struct write * report)
+mat_struct *submatrix(struct ilist *index, mat_struct * GM_all,
+ struct write *report)
{
// Local variables
int n = index->n_values; // # of selected points
@@ -506,7 +513,7 @@
switch (i3) {
case FALSE:
// Cell value diffs estimation using linear variogram
- dr[2] = 0.0; // !!! optimalize - not to be necessary to use this line
+ dr[2] = 0.0; // !!! optimize - not to be necessary to use this line
break;
case TRUE:
dr[2] = *(r + 2) - *(r0 + 2); // dz
@@ -571,7 +578,6 @@
for (i = 0; i < n; i++) { // for each input point:
*vt = *vo; // subval = selected original
*wt = *wo; // sub weight = selected original
- //G_debug(0,"%f %f", *vt, *wt);
// go to the next:
if (n_ind == 0) { // use all points:
@@ -594,9 +600,29 @@
return rslt_OK->vals[0];
}
+// find center
+double find_center(double r, double res)
+{
+ int n;
+ double center;
+
+ n = floor(r / res);
+ center = (n + 0.5) * res;
+
+ return center;
+}
+
+// calculate a cell for validation
+void adjacent_cell(int i3, double *r, struct reg_par *reg, double *cell)
+{
+ *cell = find_center(*r, reg->ew_res);
+ *(cell + 1) = find_center(*(r + 1), reg->ns_res);
+ *(cell + 2) = i3 == TRUE ? find_center(*(r + 2), reg->bt_res) : 0.;
+}
+
// validation
void crossvalidation(struct int_par *xD, struct points *pnts,
- struct parameters *var_par)
+ struct parameters *var_par, struct reg_par *reg)
{
int n = pnts->n; // # of input points
double *r = pnts->r; // coordinates of points
@@ -607,7 +633,7 @@
double ratio = type == 3 ? xD->aniso_ratio : 1.; // anisotropic ratio
mat_struct *GM = var_par->GM; // GM = theor_var(distances)
- int i;
+ int i, direction;
int n_vals;
double max_dist =
type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
@@ -616,13 +642,23 @@
double *search; // coordinates of the search point
struct ilist *list;
+
+ // coordinates of the adjacent cell
+ double *cell;
+
+ cell = (double *)G_malloc(3 * sizeof(double));
+
mat_struct *GM_sub;
mat_struct *GM_Inv, *g0, *w0;
double rslt_OK;
- struct write *crossvalid = &xD->crossvalid;
- struct write *report = &xD->report;
+ mat_struct *g0_cell, *w0_cell;
+ double rslt_OK_cell;
+
+ struct write *crossvalid = xD->crossvalid;
+ struct write *report = xD->report;
double *normal, *absval, *norm, *av;
+ double *normal_cell, *absval_cell, *norm_cell, *av_cell;
search = (double *)G_malloc(3 * sizeof(double));
normal = (double *)G_malloc(n * sizeof(double));
@@ -630,50 +666,74 @@
norm = &normal[0];
av = &absval[0];
+ normal_cell = (double *)G_malloc(n * sizeof(double));
+ absval_cell = (double *)G_malloc(n * sizeof(double));
+ norm_cell = &normal_cell[0];
+ av_cell = &absval_cell[0];
+
FILE *fp;
fp = fopen(crossvalid->name, "w");
+ G_message(_("Crossvalidation..."));
for (i = 0; i < n; i++) { // for each input point [r0]:
list = G_new_ilist(); // create list of overlapping rectangles
search = &pnts->r[3 * i];
if (i3 == TRUE) {
list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
+ direction = 0;
}
else {
list = find_NNs_within(2, search, pnts, max_dist, max_dist_vert);
+ direction = 12;
}
n_vals = list->n_values; // # of overlapping rectangles
if (n_vals > 0) { // if positive:
- correct_indices(list, r, pnts, var_par);
+ correct_indices(direction, list, r, pnts, var_par);
GM_sub = submatrix(list, GM, &xD->report); // create submatrix using indices
GM_Inv = G_matrix_inverse(GM_sub); // inverse matrix
G_matrix_free(GM_sub);
+ // calculate cell
+ adjacent_cell(i3, r, reg, cell);
+
g0 = set_up_g0(xD, pnts, list, r, var_par); // Diffs inputs - unknowns (incl. cond. 1))
w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
+ g0_cell = set_up_g0(xD, pnts, list, cell, var_par); // Diffs inputs - unknowns (incl. cond. 1))
+ w0_cell = G_matrix_product(GM_Inv, g0_cell); // Vector of weights, condition SUM(w) = 1 in last row
+
+
G_matrix_free(g0);
G_matrix_free(GM_Inv);
+ G_matrix_free(g0_cell);
+
rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
G_matrix_free(w0);
- //Create output
+ rslt_OK_cell = result(pnts, list, w0_cell); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0_cell);
+
+ // Create output
*norm = rslt_OK - *vals; // differences between input and interpolated values
*av = fabs(*norm); // absolute values of the differences (quantile computation)
+ *norm_cell = rslt_OK_cell - *vals; // differences between input and interpolated values
+ *av_cell = fabs(*norm_cell); // absolute values of the differences (quantile computation)
+
if (xD->i3 == TRUE) { // 3D interpolation:
- fprintf(fp, "%d %.3f %.3f %.2f %f %f %f\n", i, *r, *(r + 1),
- *(r + 2) / ratio, pnts->invals[i], rslt_OK, *norm);
+ fprintf(fp, "%d %.3f %.3f %.2f %f %f %f %f %f\n", i, *r,
+ *(r + 1), *(r + 2) / ratio, pnts->invals[i], rslt_OK,
+ *norm, rslt_OK_cell, *norm_cell);
}
else { // 2D interpolation:
- fprintf(fp, "%d %.3f %.3f %f %f %f\n", i, *r, *(r + 1), *vals,
- rslt_OK, *norm);
+ fprintf(fp, "%d %.3f %.3f %f %f %f %f %f\n", i, *r, *(r + 1),
+ *vals, rslt_OK, *norm, rslt_OK_cell, *norm_cell);
}
} // end if: n_vals > 0
@@ -687,6 +747,8 @@
vals++;
norm++;
av++;
+ norm_cell++;
+ av_cell++;
G_free_ilist(list); // free list memory
} // end i for loop
@@ -696,15 +758,160 @@
crossvalid->name);
if (report->name) {
- double quant95;
+ double quant95, quant95_cell;
+ report->fp = fopen(report->name, "a");
+
fprintf(report->fp,
"\n************************************************\n");
fprintf(report->fp, "*** Cross validation results ***\n");
test_normality(n, normal, report);
+ test_normality(n, normal_cell, report);
- fprintf(report->fp, "Quantile of absolute values\n");
+ fprintf(report->fp, "Quantile of absolute values (points)\n");
quant95 = quantile(0.95, n, absval, report);
+
+ fprintf(report->fp, "Quantile of absolute values (cells)\n");
+ quant95_cell = quantile(0.95, n, absval_cell, report);
}
}
+
+struct ilist *list_NN(struct int_par *xD, double *r0, struct points *pnts,
+ double max_dist, double max_dist_vert)
+{
+ // local variables
+ int i3 = xD->i3;
+
+ struct ilist *list;
+
+ list = G_new_ilist(); // create list of overlapping rectangles
+
+ if (i3 == TRUE) { // 3D kriging:
+ list = find_NNs_within(3, r0, pnts, max_dist, max_dist_vert);
+ }
+ else { // 2D kriging:
+ list = find_NNs_within(2, r0, pnts, max_dist, max_dist_vert);
+ }
+
+ return list;
+}
+
+int compare_NN(struct ilist *list, struct ilist *list_new, int modified)
+{
+ // local variables
+ int n = list->n_values, n_new = list_new->n_values;
+ double *list_value = list->value;
+ double *list_new_value = list_new->value;
+
+ int i, next = 0; // the samples are different
+
+ if (n == n_new) {
+ for (i = 0; i < n; i++) {
+ if (list_value[i] != list_new_value[i] - modified) {
+ goto change; // they are indeed
+ }
+ }
+ next = 1; // they are identical
+ }
+
+ change:
+ return next;
+}
+
+void make_subsamples(struct int_par *xD, struct ilist *list, double *r0,
+ int row, int col, struct points *pnts,
+ struct parameters *var_par, struct krig_pars *krig)
+{
+ // Local variables
+ int i3 = xD->i3;
+ double *vals = pnts->invals;
+ struct write *report = &xD->report;
+
+ int direction;
+ mat_struct *GM_sub;
+
+ direction = i3 == TRUE ? 0 : 12;
+
+ if (list->n_values > 1) { // positive # of selected points:
+ correct_indices(direction, list, r0, pnts, var_par);
+
+ GM_sub = submatrix(list, krig->GM, report); // make submatrix for selected points
+ krig->GM_Inv = G_matrix_inverse(GM_sub); // invert submatrix
+ G_matrix_free(GM_sub);
+ }
+ else if (list->n_values == 1) {
+ G_matrix_set_element(krig->rslt, row, col, vals[list->value[0] - 1]); // Estimated cell/voxel value rslt_OK = w x inputs
+ }
+ else if (list->n_values == 0) {
+ report_error(report);
+ G_fatal_error(_("This point does not have neighbours in given radius..."));
+ } // end else: error
+
+ //G_free_ilist(list); // free list memory
+}
+
+double interpolate(struct int_par *xD, struct ilist *list, double *r0,
+ struct points *pnts, struct parameters *var_par,
+ struct krig_pars *krig)
+{
+ double rslt;
+ mat_struct *g0, *w0;
+
+ g0 = set_up_g0(xD, pnts, list, r0, var_par); // Diffs inputs - unknowns (incl. cond. 1))
+ w0 = G_matrix_product(krig->GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
+
+ G_matrix_free(g0);
+
+ rslt = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
+
+ return rslt;
+}
+
+double trend(double *r0, struct output *out, int function, struct int_par *xD)
+{
+ int i3 = xD->i3;
+ double a, b, c, d, value;
+ double ratio;
+
+ ratio = i3 == FALSE || function == 5 ? 1. : xD->aniso_ratio;
+
+ a = out->trend[0];
+ b = out->trend[1];
+ c = out->trend[2];
+ d = out->trend[3];
+ value = a * *r0 + b * *(r0 + 1) + c * *(r0 + 2) / ratio + d;
+
+ return value;
+}
+
+int new_sample(struct int_par *xD, struct ilist *list, struct ilist *list_new,
+ struct points *pnts, int dep, int row, int col, double *r0,
+ double max_dist, double max_dist_vert, struct reg_par *reg,
+ struct parameters *var_par, struct krig_pars *krig,
+ int *new_matrix)
+{
+ int next;
+
+ // free memory
+ if (krig->modified == 0) {
+ G_free_ilist(list_new); // new list
+ }
+ G_matrix_free(krig->GM_Inv); // old inverse matrix
+
+ // make new...
+ list = list_NN(xD, r0, pnts, max_dist, max_dist_vert); //... list
+ cell_centre(col, row, dep, xD, reg, r0, var_par);
+ make_subsamples(xD, list, r0, dep * reg->nrows + row, col, pnts, var_par,
+ krig);
+
+ // options:
+ krig->modified = 0; // indices are not corrected
+ krig->new = TRUE; //
+ krig->first = TRUE; // new sample is being processed
+ *new_matrix++; // counter of skipped matrices
+ next = 2; // interpolate using new subsample
+
+ return next;
+}
Modified: grass-addons/grass7/vector/v.kriging/utils_raster.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_raster.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/utils_raster.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -10,58 +10,82 @@
out->dcell = (DCELL *) Rast_allocate_buf(DCELL_TYPE);
out->fd_2d = Rast_open_new(out->name, DCELL_TYPE); /* open output raster */
if (out->fd_2d < 0) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Unable to create 2D raster <%s>"), out->name);
}
}
/* 3D Raster layer */
if (xD->i3 == TRUE) {
out->fd_3d = (RASTER3D_Map *) Rast3d_open_cell_new(out->name, DCELL_TYPE, RASTER3D_USE_CACHE_XYZ, ®->reg_3d); // initialize pointer to 3D region
+
if (out->fd_3d == NULL) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
+ report_error(report);
G_fatal_error(_("Unable to create 3D raster <%s>"), out->name);
}
}
}
-int write2layer(struct int_par *xD, struct reg_par *reg, struct output *out,
- unsigned int col, unsigned int row, unsigned int dep,
- double rslt_OK)
+void write2layer(struct int_par *xD, struct reg_par *reg, struct output *out,
+ mat_struct * rslt)
{
// Local variables
int i3 = xD->i3;
+ int ndeps = reg->ndeps, nrows = reg->nrows, ncols = reg->ncols;
struct write *report = &xD->report;
+ int col, row, dep;
int pass = 0; /* Number of processed cells */
+ double trend, value, r0[3];
+ int nulval = 0;
- switch (i3) {
- case FALSE: /* set value to cell (2D) */
- out->dcell[col] = (DCELL) (rslt_OK);
- pass++;
- if (col == reg->ncols - 1)
- Rast_put_row(out->fd_2d, out->dcell, DCELL_TYPE);
- break;
+ for (dep = 0; dep < ndeps; dep++) {
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ value = G_matrix_get_element(rslt, dep * nrows + row, col);
+ if (value == 0.) {
+ nulval++;
+ }
- case TRUE: /* set value to voxel (based on part of r3.gwflow (Soeren Gebbert)) */
- if (Rast3d_put_double(out->fd_3d, col, row, dep, rslt_OK) == 0) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp,
- "Error (see standard output). Process killed...");
- fclose(report->fp);
+ switch (i3) {
+ case FALSE: /* set value to cell (2D) */
+ out->dcell[col] = (DCELL) (value);
+ break;
+
+ case TRUE: /* set value to voxel (based on part of r3.gwflow (Soeren Gebbert)) */
+ if (Rast3d_put_double(out->fd_3d, col, row, dep, value) ==
+ 0) {
+ report_error(report);
+ G_fatal_error(_("Error writing cell (%d,%d,%d) with value %f."),
+ row, col, dep, value);
+ }
+
+ break;
+ }
+ pass++;
+ } // end col
+ if (i3 == FALSE) {
+ Rast_put_row(out->fd_2d, out->dcell, DCELL_TYPE);
}
- G_fatal_error(_("Error writing cell (%d,%d,%d) with value %f, nrows = %d"),
- row, col, dep, rslt_OK, reg->nrows);
+ } // end row
+ } // end dep
+
+ switch (i3) {
+ case TRUE:
+ if (!Rast3d_close(out->fd_3d)) { // Close 3D raster map
+ G_fatal_error(_("Something went wrong with 3D raster: the pointer disappeared before closing..."));
}
- pass++;
break;
+ case FALSE:
+ Rast_close(out->fd_2d); // Close 2D raster map
+ break;
}
- return pass;
+
+ if (nulval > 0) {
+ G_message(_("0 values: %d"), nulval);
+ }
+
+ if (pass != ndeps * ncols * nrows) {
+ G_fatal_error(_("The number of processed cells (%d) is smaller than total number of cells (%d)..."),
+ pass, ndeps * ncols * nrows);
+ }
}
Modified: grass-addons/grass7/vector/v.kriging/utils_write.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_write.c 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/utils_write.c 2017-06-08 10:32:49 UTC (rev 71171)
@@ -20,8 +20,10 @@
void write2file_basics(struct int_par *xD, struct opts *opt)
{
- struct write *report = &xD->report;
+ struct write *report = xD->report;
+ report->fp = fopen(report->name, "a");
+
fprintf(report->fp, "************************************************\n");
fprintf(report->fp, "*** Input ***\n\n");
fprintf(report->fp, "- vector layer: %s\n", opt->input->answer);
@@ -51,8 +53,10 @@
void write2file_vector(struct int_par *xD, struct points *pnts)
{
struct select in_reg = pnts->in_reg;
- struct write *report = &xD->report;
+ struct write *report = xD->report;
+ report->fp = fopen(report->name, "a");
+
fprintf(report->fp, "\n");
fprintf(report->fp, "************************************************\n");
fprintf(report->fp, "*** Input vector layer properties ***\n\n");
@@ -135,7 +139,7 @@
double *h, *vert;
double *c;
double *gamma;
- struct write *report = &xD->report;
+ struct write *report = xD->report;
char type[12];
@@ -234,24 +238,24 @@
switch (type) {
case 0: // horizontal variogram
fp = fopen("variogram_hz_tmp.txt", "w");
- if (xD->report.write2file) { // write name of report file
- file_length = strlen(xD->report.name);
+ if (xD->report->name) { // write name of report file
+ file_length = strlen(xD->report->name);
if (file_length < 4) { // 4 types of variogram
G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
}
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report->name);
}
fprintf(fp, "%d\n", var_pars->type); // write # of lags
break;
case 1: // vertical variogram
fp = fopen("variogram_vert_tmp.txt", "w");
- if (xD->report.write2file == TRUE) { // write name of report file
- file_length = strlen(xD->report.name);
+ if (xD->report->name) { // write name of report file
+ file_length = strlen(xD->report->name);
if (file_length < 3) {
G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
}
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report->name);
}
fprintf(fp, "%d\n", var_pars->type); // write type
break;
@@ -259,12 +263,12 @@
case 2: // bivariate variogram
fp = fopen("variogram_final_tmp.txt", "w");
- if (xD->report.write2file == TRUE) { // write name of report file
- file_length = strlen(xD->report.name);
+ if (xD->report->name) { // write name of report file
+ file_length = strlen(xD->report->name);
if (file_length < 4) { // 4 types of variogram
G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
}
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report->name);
}
fprintf(fp, "%d\n", var_pars->type); // write type
@@ -279,12 +283,12 @@
case 3: // anisotropic variogram
fp = fopen("variogram_final_tmp.txt", "w");
- if (xD->report.name) { // write name of report file
- file_length = strlen(xD->report.name);
+ if (xD->report->name) { // write name of report file
+ file_length = strlen(xD->report->name);
if (file_length < 4) { // 4 types of variogram
G_fatal_error(_("File name must contain more than 4 characters...")); // todo: error
}
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report->name);
}
fprintf(fp, "%d\n", var_pars->type); // write type
fprintf(fp, "%f\n", xD->aniso_ratio); // write ratio of anisotropy
@@ -334,3 +338,11 @@
fclose(fp);
}
+
+void report_error(struct write *report)
+{
+ if (report->name) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+}
Modified: grass-addons/grass7/vector/v.kriging/v.kriging.html
===================================================================
--- grass-addons/grass7/vector/v.kriging/v.kriging.html 2017-06-08 07:59:33 UTC (rev 71170)
+++ grass-addons/grass7/vector/v.kriging/v.kriging.html 2017-06-08 10:32:49 UTC (rev 71171)
@@ -45,7 +45,9 @@
</pre></div>
<p>
-Commands based on the <a href="http://grass.osgeo.org/download/sample-data/" target="_blank">dataset</a> of <b>Slovakia 3D precipitation</b> (<i>Mitasova and Hofierka, 2004</i>). Another examples of 3D interpolation are available in (<i>Stopkova, 2014</i>).
+Commands based on the <a href="http://grass.osgeo.org/download/sample-data/" target="_blank">dataset</a> of <b>Slovakia
+3D precipitation</b> (<i>Mitasova and Hofierka, 2004</i>). For more detailed information check <a href="v.kriging.pdf"
+target="_blank">case studies</a>. Another examples of 3D interpolation are available in (<i>Stopkova, 2014</i>).
<div class="code"><pre>
v.kriging phase=initial in=precip3d at PERMANENT ic=precip report=precip3d.txt file=png --o
@@ -71,7 +73,10 @@
out=name crossval=crossval_file.txt -2
</pre></div>
-<p> Commands based on 500 random points extracted from input points of Digital Elevation Model (DEM) <i>elev_lid792_randpts</i> from the <b>North Carolina <a href="http://grass.osgeo.org/download/sample-data/" target="_blank">dataset</a></b> (<i>Neteler and Mitasova, 2004</i>).
+<p> Commands based on 500 random points extracted from input points of Digital Elevation Model (DEM)
+<i>elev_lid792_randpts</i> from the <b>North Carolina
+<a href="http://grass.osgeo.org/download/sample-data/" target="_blank">dataset</a></b> (<i>Neteler and Mitasova, 2004</i>).
+See the <a href="v.kriging.pdf" target="_blank">case studies</a>.
<div class="code"><pre>
v.kriging phase=initial in=elev_lid792_selected ic=value azimuth=45. td=45. \
@@ -84,7 +89,6 @@
<h2>TODO</h2>
<ul>
-<li>final phase still needs <b>optimization</b></li>
<li><b>anisotropy</b> in horizontal direction missing
<li>current version is suitable just for <b>metric coordinate systems</b>
<li>enable <b>mask usage</b>
@@ -92,6 +96,12 @@
<li><b>2D interpolation from 3D input layer</b> needs to be rebuilt (especially in case that there are too many points located on identical horizontal coordinates with different elevation)
</ul>
+<h2>Recommendations</h2>
+<ul>
+<li> In case of too much <i>warnings</i> about input points that have "<b>less than 2 neighbours in its closest surrounding</b>. The perimeter of the surrounding will be increased...", please consider shorter variogram range.
+<li> Save just figures with theoretical variogram (using <i>file=extension</i> in the middle and final phase). Experimental variograms are included in the theoretical variogram plot and separate "experimental" plots can be just temporal.
+</ul>
+
<h2>REFERENCES</h2>
<table>
More information about the grass-commit
mailing list