[GRASS-SVN] r65504 - grass-addons/grass7/vector/v.kriging
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 20 04:32:18 PDT 2015
Author: evas
Date: 2015-06-20 04:32:18 -0700 (Sat, 20 Jun 2015)
New Revision: 65504
Added:
grass-addons/grass7/vector/v.kriging/geostat.c
grass-addons/grass7/vector/v.kriging/getval.c
grass-addons/grass7/vector/v.kriging/quantile.c
grass-addons/grass7/vector/v.kriging/stat.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
Log:
v.kriging: kd-tree (PCL library) replaced by R-tree (GRASS GIS) and 2D kriging enabled (currently: inaccurate results for sparse and spatially heterogeneous data)
Added: grass-addons/grass7/vector/v.kriging/geostat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/geostat.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/geostat.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,615 @@
+#include "local_proto.h"
+
+/* experimental variogram
+ * based on 2D variogram (alghalandis.com/?page_id=463)) */
+ void E_variogram(int type, struct int_par *xD, struct points *pnts, struct reg_par *reg, struct var_par *pars)
+{
+ // Variogram properties
+ struct parameters *var_pars, *var_par_hz, *var_par_vert;
+ 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
+ break;
+ case 1: // vertica variogram
+ var_pars = &pars->vert; // 1: vertical variogram
+ break;
+ case 2: // bivariate variogram
+ var_pars = &pars->fin;
+
+ max_dist = var_pars->horizontal.max_dist;
+ max_dist_vert = var_pars->vertical.max_dist;
+ radius_vert = SQUARE(max_dist_vert);
+ break;
+ case 3: //anisotropic variogram
+ var_pars = &pars->fin; // default: final variogram
+ max_dist = max_dist_vert = var_pars->max_dist;
+ radius = radius_vert = SQUARE(max_dist);
+ break;
+ }
+
+ if (type < 2) {
+ max_dist = var_pars->max_dist;
+ }
+ radius = SQUARE(max_dist);
+
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *pnts_r = pnts->r; // xyz coordinates of input points
+ double *r; // pointer to xyz coordinates
+ double *vals = pnts->invals; // values to be used for interpolation
+ int i3 = xD->i3;
+ int phase = xD->phase; // phase: initial / middle / final
+ int function = var_pars->function; // type of theoretical variogram
+
+ int nLag; // # of horizontal bins
+ int nLag_vert; // # of vertical bins
+ double *vert; // pointer to vertical bins
+ double dir = var_pars->dir; // azimuth
+ double td = var_pars->td; // angle tolerance
+ double lag = var_pars->lag; // size of horizontal bin
+ double lag_vert;
+
+ if (type == 2) { // bivariate variogram
+ nLag = var_pars->horizontal.nLag;
+ nLag_vert = var_pars->vertical.nLag;
+
+ dir = var_pars->dir;
+ td = var_pars->td;
+ lag = var_pars->horizontal.lag;
+ lag_vert = var_pars->vertical.lag;
+ }
+
+ else { // univariate variogram
+ nLag = var_pars->nLag;
+ dir = var_pars->dir;
+ td = var_pars->td;
+ lag = var_pars->lag;
+ }
+
+ // depend on variogram type:
+ int k; // index of variogram to be computed
+ int nrows = nLag;
+ int ncols = type == 2 ? nLag_vert : 1;
+
+ struct write *report = &xD->report;
+
+ // Variogram processing variables
+ int s; // index of horizontal segment
+ int b; // index of verical segment (bivariate only)
+ int i, j; // indices of the input points
+ int err0=0; // number of invalid values
+
+ double *dr; // coordinate differences of each couple
+ double *h; // distance of boundary of the horizontal segment
+ double tv; // bearing of the line between the couple of the input points
+ double ddir; // the azimuth of computing variogram
+ double rv; // radius of the couple of points
+ double drv; // distance between the points
+ double rvh; // difference between point distance and the horizontal segment boundary
+ double dv; // difference of the values to be interpolated that are located on the couple of points
+
+ struct ilist *list; // list of selected nearest neighbours
+ int n_vals; // # of selected NN
+
+ double *i_vals, *j_vals; // values located on the point couples
+
+ int n1 = n+1; // number of points + 1
+ int *ii; // difference of indices between rellevant input points
+
+ double gamma_lag; // sum of dissimilarities in one bin
+ double cpls; // # of dissimilarities in one bin
+ double gamma_E; // average of dissimilarities in one bin (element of gamma matrix)
+ mat_struct *gamma_M; // gamma matrix (hz, vert or bivar)
+ mat_struct *c_M; // matrix of # of dissimilarities
+ double *gamma; // pointer to gamma matrix
+ double *c; // pointer to c matrix
+
+ unsigned int percents = 50;
+
+ double gamma_sum; // sum of gamma elements (non-nan)
+ int gamma_n; // # of gamma elements (non-nan)
+ double gamma_sill; // sill
+
+ /* Allocated vertices and matrices:
+ * --------------------------------
+ * dr - triple of coordinate differences (e.g. dx, dy, dz)
+ * h - vector of length pieces [nL x 1]
+ * vert - vector of vertical pieces [nZ x 1] (3D only)
+ * c - matrix containing number of points in each sector defined by h (and vert) [nL x nZ; in 2D nZ = 1]
+ * gama - matrix of values of experimental variogram
+ */
+
+ // allocate:
+ dr = (double *) G_malloc(3 * sizeof(double)); // vector of coordinate differences
+ if (type != 2) {
+ var_pars->h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
+ }
+
+ if (type == 2) {
+ var_pars->horizontal.h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
+ var_pars->vertical.h = (double *) G_malloc(nLag_vert * sizeof(double)); // vector of bins
+ }
+
+ // control initialization:
+ if (dr == 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);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
+ // initialize:
+ c_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix of counts
+ 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);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
+
+ // set up pointers
+ c = &c_M->vals[0];
+ gamma = &gamma_M->vals[0];
+
+ // set up starting values
+ gamma_sum = 0.;
+ gamma_n = 0;
+
+ if (percents) {
+ G_percent_reset();
+ }
+
+ if (type == 2) {
+ vert = &var_pars->vertical.h[0];
+ }
+
+ /* *** Experimental variogram computation *** */
+ for (b = 0; b < ncols; b++) { // for each vertical lag...
+ if (type == 2) { // just in case that the variogram is vertical
+ *vert = (b+1) * lag_vert; // lag distance
+ }
+
+ h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0]; // ... horizontal lags
+
+ 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 ...
+ 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
+
+ /* Compute variogram for points in relevant neighbourhood */
+ for (i = 0; i < n-1; i++) { // for each input point...
+ switch (type) { // find NNs according to variogram type
+ case 0: // horizontal variogram
+ list = find_NNs_within(2, i, pnts, max_dist, -1);
+ break;
+ case 1: // vertical variogram
+ list = find_NNs_within(1, i, pnts, -1, max_dist);
+ break;
+ default: // anisotropic or bivariate variogram
+ list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+ break;
+ }
+
+ n_vals = list->n_values; // # of input values located on NN
+ if (n_vals > 0) {
+ correct_indices(i3, 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
+
+ // Variogram processing
+ if (type == 1) { // vertical variogram:
+ tv = 0.5 * PI - zenith_angle(dr); // compute zenith angle instead of bearing
+ td = 0.5 * PI; // todo: repair, 0.25 usually
+ }
+ else { // hz / aniso / bivar variogram:
+ tv = atan2(*(dr+1), *dr); // bearing
+ }
+
+ ddir = tv - dir; // difference between bearing and azimuth
+ if (fabsf(ddir) <= 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
+
+ if (type == 1 || type == 3) { // vertical or anisotropic variogram:
+ rv += SQUARE(*(dr+2)); // consider also vertical direction
+ }
+
+ rvh = sqrt(rv) - *h; // the difference between distance and lag
+ if (rv <= radius && fabsf(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 = *(dr+2) - *vert; // compare vertical
+
+ if (fabsf(rvh) <= lag_vert) { // elevation test: vertical lag
+ goto delta_V;
+ }
+ else {
+ goto end;
+ }
+ }
+ 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
+ cpls += 1.; // count of elements
+ } // end distance test: rv <= radius ^ |rvh| <= lag
+ } // end angle test: |ddir| <= td
+ } // end i test: *ii > i
+ end:
+ ii++; // go to the next index
+ j_vals += *ii - *(ii-1); // go to the next value
+ } // 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);
+ }
+ 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
+ } // end for loop i: variogram computation for each i-th search point
+
+ if (isnan(gamma_lag) || cpls == 0.0) { // empty lags:
+ err0++; // error indicator
+ goto gamma_isnan;
+ }
+
+ *gamma = 0.5 * gamma_lag / cpls; // element of gamma matrix
+ *c = cpls; // # of values that have been used to compute gamma_e
+
+ gamma_sum += *gamma; // sum of gamma elements (non-nan)
+ gamma_n++; // # of gamma elements (non-nan)
+
+ gamma_isnan: // there are no available point couples to compute the dissimilarities in the lag:
+ h++;
+ c++;
+ gamma++;
+ } // end for loop s
+
+ if (type == 2) { // vertical variogram:
+ vert++;
+ }
+ } // 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);
+ }
+ G_fatal_error(_("Unable to compute experimental variogram..."));
+ } // end error
+
+ var_pars->gamma = G_matrix_copy(gamma_M);
+ var_pars->gamma_sum = gamma_sum;
+ var_pars->gamma_n = gamma_n;
+
+ plot_experimental_variogram(xD, var_pars);
+
+ if (phase < 2) { // initial and middle phase:
+ sill(var_pars); // compute sill
+ }
+
+ if (report->write2file == TRUE) { // 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)
+{
+ char *variogram;
+
+ // report
+ if (report->write2file = TRUE) { // report file available:
+ time(&report->now); // write down time of start
+ if (type != 1) {
+ fprintf(report->fp, "\nComputation of theoretical variogram started on %s\n", ctime(&report->now));
+ }
+ }
+
+ // set up:
+ var_pars->type = type; // hz / vert / bivar / aniso
+ var_pars->const_val = 0; // input values are not constants
+ switch (type) {
+ case 0: // horizontal variogram
+ if (i3 == TRUE) { // 3D interpolation (middle phase)
+ variogram = opt.function_var_hz->answer; // function type available:
+ if (strcmp(variogram, "linear") != 0) { // nonlinear
+ 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:
+ linear_variogram(var_pars, report);
+ }
+ }
+ else { // 2D interpolation (final phase)
+ variogram = opt.function_var_final->answer; // function type available:
+ if (strcmp(variogram, "linear") != 0) { // nonlinear
+ 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);
+ }
+ }
+ else { // function type not available:
+ linear_variogram(var_pars, report);
+ }
+ }
+
+ if (report->name) {
+ 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);
+ fprintf(report->fp, "Range: %f\n", var_pars->h_range);
+ }
+ break;
+
+ case 1: // vertical variogram
+ 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;
+ if (report->name) {
+ 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)
+ 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)
+ linear_variogram(var_pars, report);
+ }
+
+ else { // function type was set up by the user:
+ var_pars->horizontal.nugget = atof(opt.nugget_final->answer);
+ var_pars->horizontal.h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_pars->horizontal.sill = atof(opt.sill_final->answer);
+ }
+
+ var_pars->vertical.nugget = atof(opt.nugget_final_vert->answer);
+ var_pars->vertical.h_range = atof(opt.range_final_vert->answer);
+ if (opt.sill_final_vert->answer) {
+ var_pars->vertical.sill = atof(opt.sill_final_vert->answer);
+ }
+
+ var_pars->horizontal.function = set_function(opt.function_var_final->answer, report);
+ var_pars->vertical.function = set_function(opt.function_var_final_vert->answer, report);
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Nugget effect (hz): %f\n", var_pars->horizontal.nugget);
+ fprintf(report->fp, "Sill (hz): %f\n", var_pars->horizontal.sill);
+ fprintf(report->fp, "Range (hz): %f\n", var_pars->horizontal.h_range);
+ fprintf(report->fp, "Function: %s\n\n", opt.function_var_hz->answer);
+ fprintf(report->fp, "Nugget effect (vert): %f\n", var_pars->vertical.nugget);
+ fprintf(report->fp, "Sill (vert): %f\n", var_pars->vertical.sill);
+ fprintf(report->fp, "Range (vert): %f\n", var_pars->vertical.h_range);
+ fprintf(report->fp, "Function: %s\n", opt.function_var_vert->answer);
+ }
+ }
+
+ plot_var(i3, TRUE, var_pars); // Plot variogram using gnuplot
+ break;
+
+ case 3: // univariate (just final phase)
+ variogram = opt.function_var_final->answer;
+ if (strcmp(variogram, "linear") != 0) { // nonlinear variogram:
+ 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 (i3 == TRUE) { // 3D interpolation:
+ fprintf(report->fp, "Parameters of anisotropic variogram:\n");
+ }
+ else {
+ fprintf(report->fp, "Parameters of 2D 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);
+ }
+ }
+ else {
+ linear_variogram(var_pars, report);
+ }
+ break;
+ }
+
+ if (type != 2) {
+ var_pars->function = set_function(variogram, report);
+ plot_var(i3, FALSE, var_pars); // Plot variogram using gnuplot
+ }
+}
+
+void ordinary_kriging(struct int_par *xD, struct reg_par *reg, struct points *pnts, struct var_par *pars, struct output *out)
+{
+ // Local variables
+ int i3 = xD->i3;
+ int n = pnts->n; // number of input points
+ double *r = pnts->r; // xyz coordinates of input points
+ double *vals = pnts->invals; // values to be used for interpolation
+ struct write *report = &xD->report;
+ struct write *crossvalid = &xD->crossvalid;
+ struct parameters *var_par = &pars->fin;
+
+ int type = var_par->type;
+ double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
+ double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
+ double ratio = var_par->type == 3 ? xD->aniso_ratio : 1.;
+
+ unsigned int passed=0; // number of successfully interpolated valuesy
+ unsigned int percents=50; // counter
+ unsigned int nrcd; // number of cells/voxels
+ unsigned int row, col, dep; // indices of cells/voxels
+ double rslt_OK; // interpolated value located on r0
+
+ int i, j;
+ unsigned int n1;
+
+ pnts->max_dist = var_par->lag;
+ struct ilist *list;
+ int n_vals;
+
+ 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
+
+ // Cell/voxel center coords (location of interpolated value)
+ r0 = (double *) G_malloc(3 * sizeof(double));
+
+ FILE *fp;
+
+ if (report->name) { // report file available:
+ report->write2file = TRUE;
+ report->fp = fopen(report->name, "a");
+ time(&report->now);
+ fprintf(report->fp, "Interpolating values started on %s\n\n", ctime(&report->now));
+ }
+
+ G_message(_("Interpolating unknown values..."));
+ if (percents) {
+ G_percent_reset();
+ }
+
+ open_layer(xD, reg, out); // open 2D/3D raster
+
+ if (var_par->const_val == 1) { // input values are constant:
+ 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
+
+ // perform cross validation...
+ if (crossvalid->name) { // ... if desired
+ crossvalidation(xD, pnts, var_par);
+ }
+
+ constant_voxel_centre:
+ for (dep=0; dep < reg->ndeps; dep++) {
+ if (xD->i3 == TRUE) {
+ if (percents) {
+ G_percent(dep, reg->ndeps, 1);
+ }
+ }
+ for (row=0; row < reg->nrows; row++) {
+ if (xD->i3 == FALSE) {
+ if (percents) {
+ G_percent(row, reg->nrows, 1);
+ }
+ }
+ //#pragma omp parallel for private(col, r0, cellCent, ind, sqDist, n1, GM, GM_Inv, g0, w0, rslt_OK)
+ for (col=0; col < reg->ncols; col++) {
+
+ if (var_par->const_val == 1) { // constant input values
+ goto constant_voxel_val;
+ }
+
+ cell_centre(col, row, dep, xD, reg, r0, var_par); // coords of output point
+
+ // add cell centre to the R-tree
+ list = G_new_ilist(); // create list of overlapping rectangles
+
+ if (i3 == TRUE) { // 3D kriging:
+ list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+ }
+ else { // 2D kriging:
+ list = find_NNs_within(2, i, pnts, max_dist, max_dist_vert);
+ }
+
+ n_vals = list->n_values; // # of selected points
+ if (n_vals > 0) { // positive # of selected points:
+ correct_indices(i3, list, r0, pnts, var_par);
+
+ 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);
+
+ 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
+
+ G_matrix_free(GM_Inv);
+ G_matrix_free(g0);
+
+ rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+
+ // Create output
+ constant_voxel_val:
+ if (var_par->const_val == 1) { // constant input values:
+ rslt_OK = (double) *vals; // setup input as output
+ }
+
+ // write output to the (3D) raster layer
+ if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
+ if (report->name) { // 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..."));
+ }
+ G_free_ilist(list); // free list memory
+
+ } // end if searchRadius
+ else { // no selected points:
+ if (report->name) { // 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
+ } // end col
+ } // end row
+ } // end dep
+
+ 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;
+ }
+}
Added: grass-addons/grass7/vector/v.kriging/getval.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/getval.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/getval.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,548 @@
+#include "local_proto.h"
+
+#ifndef MAX
+#define MIN(a,b) ((a<b) ? a : b)
+#define MAX(a,b) ((a>b) ? a : b)
+#endif
+
+/* 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 select *in_reg = &pnts->in_reg;
+
+ int i, n, nrec, type, ctype;
+ struct field_info *Fi;
+
+ dbCatValArray cvarr;
+ dbDriver *Driver;
+
+ int *index;
+ double *values, *vals;
+
+ int save;
+ if (xD == NULL) {
+ save = 0;
+ }
+ else {
+ save = 1;
+ }
+
+ G_message(_("Reading values from attribute column <%s>..."), column);
+ 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, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Database connection not defined for layer %d"), field);
+ }
+
+ 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, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+ }
+
+ /* Note do not check if the column exists in the table because it may be expression */
+
+ /* TODO: only select values we need instead of all in column */
+
+ /* Select pairs key/value to array, values are sorted by key (must be integer) */
+ nrec = 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, "Error (see standard output). Process killed...");
+ 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, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Column must be numeric"));
+ }
+
+ db_close_database_shutdown_driver(Driver);
+ G_free(Fi);
+
+ /* Output cats/values list */
+ switch (in_reg->out) {
+ case 0:
+ index = NULL;
+ n = cvarr.n_values;
+ break;
+ default:
+ index = &in_reg->indices[0];
+ n = in_reg->n;
+ break;
+ }
+
+ values = (double *) G_malloc(n * sizeof(double));
+ vals = &values[0];
+
+ for (i = 0; i < n; i++) {
+ /* TODO: Only for point data:
+ * - store indexes of skipped points in read_points and compare them with database
+ */
+ if (index == NULL || (index != NULL && *index == i)) {
+ if (ctype == DB_C_TYPE_INT)
+ *vals = (double) cvarr.value[i].val.i;
+ else if (ctype == DB_C_TYPE_DOUBLE) {
+ *vals = cvarr.value[i].val.d;
+ }
+ vals++;
+ if (index != NULL)
+ index++;
+ }
+ }
+
+ 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);
+ }
+
+ db_CatValArray_free(&cvarr); // free memory of the array of categories and values
+
+ return values;
+}
+
+/* get coordinates of input points */
+void read_points(struct Map_info *map, struct reg_par *reg, struct points *point, struct int_par *xD, const char *zcol, int field, struct write *report)
+{
+ int type; // check elements of vector layer
+ int i3 = xD->i3; // 2D / 3D interpolation
+ int n; // # of all elements in vector layer
+ int nskipped; // # of skipped elements (all except points)
+ int n_in_reg = 0; // # of points within a region
+ int out_reg = 0; // # of points out of region
+ double *z_attr; // vertical coordinate from attribute value
+
+ struct line_pnts *Points; // structures to hold line *Points (map)
+
+ // create spatial index R-tree
+ point->Rtree_hz = RTreeCreateTree(-1, 0, 2); // create 2D spatial index
+ if (i3 == TRUE) { // 3D:
+ point->R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
+ point->Rtree_vert = RTreeCreateTree(-1, 0, 1); // create 1D spatial index
+ }
+
+ int i;
+ int sid=0; // id of the rectangle (will be increased by 1)
+
+ double *rx, *ry, *rz, *r; // pointers to coordinates
+ int ind = 0; // index of the point
+ int *indices; // indices of the points within the region
+ int *index; // pointer to the vector of indices
+
+ G_message(_("Reading coordinates..."));
+
+ Points = Vect_new_line_struct(); // Make point structure
+ n = Vect_get_num_primitives(map, GV_POINT); // Number of points - topology required
+
+ point->r = (double *) G_malloc(n * 3 * sizeof(double)); // coords: x0 y0 z0 ... xn yn zn
+ indices = (int *) G_malloc(n * sizeof(int)); // indices of pts within region
+
+ point->r_min = (double *) G_malloc(3 * sizeof(double)); // minimum
+ point->r_max = (double *) G_malloc(3 * sizeof(double)); // maximum
+
+ rx = &point->r[0]; // pointer to x coordinate
+ ry = &point->r[1]; // pointer to y coordinate
+ rz = &point->r[2]; // pointer to z coordinate
+ index = &indices[0]; // pointer to indices
+
+ // 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);
+ }
+
+ nskipped = 0; // # of skipped elements (lines, areas)
+
+ // Read points and set them into the structures
+ while (TRUE) {
+ type = Vect_read_next_line(map, Points, NULL);
+ if (type == -1) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to read vector map"));
+ }
+ if (type == -2) {
+ break;
+ }
+
+ // skip everything except points
+ if (type != GV_POINT) {
+ nskipped++;
+ continue;
+ }
+
+ if (isnan(Points->x[0]) || isnan(Points->y[0]) || isnan(Points->z[0])) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Input coordinates must be a number. Check point no. %d..."), ind);
+ }
+
+ // take every point in region...
+ if ((reg->west <= Points->x[0] && Points->x[0] <= reg->east) && (reg->south <= Points->y[0] && Points->y[0] <= reg->north)) {
+ *rx = Points->x[0]; // set up x coordinate
+ *ry = Points->y[0]; // set up y coordinate
+
+ insert_rectangle(2, sid, point); // R-tree node 2D (2D and 3D kriging)
+
+ if (i3 == TRUE) { // 3D interpolation:
+ if (reg->bot <= Points->z[0] && Points->z[0] <= reg->top) {
+ if (zcol == NULL) { // 3D layer:
+ *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);
+ }
+ G_fatal_error(_("Error reading input coordinates z..."));
+ }
+ }
+ else { // 2D layer with z-column:
+ *rz = *z_attr; // set up x coordinate
+ } // end else (2,5D layer)
+
+ insert_rectangle(3, sid, point); // R-tree node 3D (just 3D kriging)
+ insert_rectangle(1, sid, point); // R-tree node 1D (just 3D kriging)
+ } // end if rb, rt
+
+ else { // point is out of region:
+ goto out_of_region;
+ }
+ } // end if 3D interpolation
+
+ else { // 2D interpolation:
+ *rz = 0.; // set up z coordinate
+ }
+
+ // Find extends
+ if (n_in_reg == 0) {
+ triple(*rx, *ry, *rz, point->r_min);
+ triple(*rx, *ry, *rz, point->r_max);
+ }
+
+ else {
+ triple(MIN(*point->r_min, *rx), MIN(*(point->r_min+1), *ry), MIN(*(point->r_min+2), *rz), point->r_min);
+ triple(MAX(*point->r_max, *rx), MAX(*(point->r_max+1), *ry), MAX(*(point->r_max+2), *rz), point->r_max);
+ } // end else (extent test)
+
+ n_in_reg++; // # of points in region
+ *index = ind; // store index of point within the region
+ index++; // go to next index
+
+ rx += 3; // go to new x coordinate
+ ry += 3; // go to new y coordinate
+ rz += 3; // go to new z coordinate
+
+ sid++; // go to new ID for spatial indexing
+
+ goto finish;
+ } // end in region
+
+ out_of_region:
+ out_reg++;
+ continue;
+
+ finish:
+ if (zcol != NULL) {
+ z_attr++;
+ }
+ ind++;
+ } // end while (TRUE)
+
+ point->in_reg.total = n;
+ point->in_reg.n = n_in_reg;
+ point->in_reg.out = out_reg;
+ point->n = n_in_reg;
+
+ point->in_reg.indices = (int *) G_malloc(n_in_reg * sizeof(int));
+ memcpy(point->in_reg.indices, indices, n_in_reg * sizeof(int));
+ G_free(indices);
+
+ if (nskipped > 0) {
+ G_warning(_("%d features skipped, only points are accepted"), nskipped);
+ }
+
+ Vect_destroy_line_struct(Points);
+
+ if (n_in_reg == 0) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to read coordinates of point layer (all points are out of region)"));
+ }
+
+ if (xD->phase == 0) { // initial phase:
+ write2file_vector(xD, point); // describe properties
+ }
+}
+
+void get_region_pars(struct int_par *xD, struct reg_par *reg)
+{
+ /* Get extent and resolution of output raster 2D/3D */
+ if (xD->i3 == FALSE) { /* 2D region */
+ /* Get database region parameters */
+ Rast_get_window(®->reg_2d); /* stores the current default window in region */
+ reg->west = reg->reg_2d.west;
+ reg->east = reg->reg_2d.east;
+ reg->north = reg->reg_2d.north;
+ reg->south = reg->reg_2d.south;
+ reg->ew_res = reg->reg_2d.ew_res;
+ reg->ns_res = reg->reg_2d.ns_res;
+ reg->nrows = reg->reg_2d.rows;
+ reg->ncols = reg->reg_2d.cols;
+ reg->ndeps = 1;
+ }
+ if (xD->i3 == TRUE) { /* 3D region */
+ Rast3d_init_defaults();
+ /* Get database region parameters */
+ 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, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error("To process 3D interpolation, please set 3D region settings.");
+ }
+ reg->west = reg->reg_3d.west;
+ reg->east = reg->reg_3d.east;
+ reg->north = reg->reg_3d.north;
+ reg->south = reg->reg_3d.south;
+ reg->bot = reg->reg_3d.bottom; /* output 3D raster */
+ reg->top = reg->reg_3d.top; /* output 3D raster */
+ reg->ew_res = reg->reg_3d.ew_res;
+ reg->ns_res = reg->reg_3d.ns_res;
+ reg->bt_res = reg->reg_3d.tb_res;
+ reg->nrows = reg->reg_3d.rows;
+ reg->ncols = reg->reg_3d.cols;
+ reg->ndeps = xD->i3 == TRUE ? reg->reg_3d.depths : 1; /* 2D/3D interpolation */
+ }
+ reg->nrcd = reg->nrows * reg->ncols * reg->ndeps; /* number of cells */
+}
+
+// read values from temporary files
+void read_tmp_vals(const char *file_name, struct parameters *var_par, struct int_par *xD)
+{
+ FILE *fp;
+
+ int j, nLag_vert;
+ double lag_vert, max_dist_vert;
+ double *v_elm, *g_elm, *c_elm;
+ double sill_hz, sill_vert;
+
+
+ fp = fopen(file_name, "r");
+ if (fp == NULL) {
+ G_fatal_error(_("File is missing..."));
+ }
+
+ else { // file exists:
+ int i, type;
+ int nLag;
+ double lag, max_dist, td_hz, sill;
+ double *h, *h_elm, *c, *gamma;
+ int function;
+ int file, file_length;
+
+ for (i=0; i<2; i++) {
+ if (fscanf(fp, "%d", &file_length) == 0) { // filename length
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ if (file_length > 3) {
+ if (fscanf(fp, "%d", &file) == 0) {
+ 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;
+ }
+ 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;
+ }
+ }
+ else {
+ type = file_length;
+ goto no_file;
+ }
+ } // todo: test without report and crossval
+
+ if (fscanf(fp, "%d", &type) == 0) { // read type
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ no_file:
+
+ switch (type) {
+ case 2: // bivariate variogram
+ var_par->type = type;
+ if (fscanf(fp, "%d %lf %lf %d %lf %lf %lf", &nLag_vert, &lag_vert, &max_dist_vert, &nLag, &lag, &max_dist, &td_hz) < 7) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ var_par->horizontal.h = (double *) G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->horizontal.h[0];
+ for (i=0; i<nLag; i++) {
+ if (fscanf(fp, "%lf", h_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ h_elm++;
+ }
+
+ var_par->vertical.h = (double *) G_malloc(nLag_vert * sizeof(double));
+ v_elm = &var_par->vertical.h[0];
+ for (i=0; i<nLag_vert; i++) {
+ if (fscanf(fp, "%lf", v_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ v_elm++;
+ }
+
+ var_par->gamma = G_matrix_init(nLag, nLag_vert, nLag);
+ g_elm = &var_par->gamma->vals[0];
+ for (i=0; i<nLag_vert; i++) {
+ for (j=0; j<nLag; j++) {
+ if (fscanf(fp, "%lf", g_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ g_elm++;
+ }
+ }
+
+ if (fscanf(fp, "%lf %lf", &sill_hz, &sill_vert) < 2) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ var_par->vertical.nLag = nLag_vert;
+ var_par->vertical.lag = lag_vert;
+ var_par->vertical.max_dist = max_dist_vert;
+ var_par->vertical.sill = sill_vert;
+
+ var_par->horizontal.nLag = nLag;
+ var_par->horizontal.lag = lag;
+ var_par->horizontal.max_dist = max_dist;
+ var_par->horizontal.sill = sill_hz;
+ break;
+
+ default:
+ if (type == 3) { // anisotropic variogram:
+ if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0) { // anisotropic ratio
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ }
+
+ if (fscanf(fp, "%d %lf %lf", &nLag, &lag, &max_dist) < 3) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ 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];
+
+ var_par->gamma = G_matrix_init(nLag, 1, nLag);
+ gamma = &var_par->gamma->vals[0];
+
+ for (i=0; i<nLag; i++) {
+ if (fscanf(fp, "%lf %lf", h_elm, gamma) < 2) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ h_elm++;
+ gamma++;
+ }
+
+ if (fscanf(fp, "%lf", &sill) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ var_par->sill = sill;
+ break;
+ }
+
+ if (type != 2) {
+ var_par->type = type;
+ var_par->nLag = nLag;
+ var_par->lag = lag;
+ var_par->max_dist = max_dist;
+ var_par->td = td_hz;
+ }
+ }
+ fclose(fp);
+}
Added: grass-addons/grass7/vector/v.kriging/quantile.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/quantile.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/quantile.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,223 @@
+#include "local_proto.h"
+
+/* These functions are taken from the module r.quantile (Clemens, G.) */
+
+struct bin
+{
+ unsigned long origin;
+ double minimum, maximum;
+ int base, count;
+};
+
+static double minimum, maximum;
+static int num_quants;
+static double *quants;
+static int num_slots;
+static unsigned int *slots;
+static double slot_size;
+static unsigned long total;
+static int num_values;
+static unsigned short *slot_bins;
+static int num_bins;
+static struct bin *bins;
+static double *values;
+
+static inline int get_slot(double c)
+{
+ int i = (int) floor((c - minimum) / slot_size);
+
+ if (i < 0)
+ i = 0;
+ if (i > num_slots - 1)
+ i = num_slots - 1;
+ return i;
+}
+
+static inline double get_quantile(int n)
+{
+ return (double)total * quants[n]; // # of lower values
+}
+
+static void get_slot_counts(int n, double *data)
+{
+ int i;
+
+ total = 0;
+ for (i = 0; i < n; i++) {
+ int j;
+
+ // todo nan
+ j = get_slot(data[i]);
+
+ slots[j]++; // rise value of j-th slot
+ total++;
+ }
+ //G_percent(i, n, 2);
+}
+
+static void initialize_bins(void)
+{
+ int slot; // index of slot
+ double next; // percentile
+ int bin = 0; // adjacent bin
+ unsigned long accum = 0;
+ int quant = 0; // index of percentile
+
+ num_values = 0;
+ next = get_quantile(quant); // # of lower values
+
+ for (slot = 0; slot < num_slots; slot++) {
+ unsigned int count = slots[slot]; // assign # of elements in each slots to the count
+ unsigned long accum2 = accum + count; // total # of elements
+
+ if (accum2 > next) { // # of values is greater than percentile
+ struct bin *b = &bins[bin]; // make bin
+ slot_bins[slot] = ++bin;
+
+ b->origin = accum; // origin of bin = total # of elements
+ b->base = num_values;
+ b->count = 0;
+ b->minimum = minimum + slot_size * slot;
+ b->maximum = minimum + slot_size * (slot + 1);
+
+ while (accum2 > next)
+ next = get_quantile(++quant);
+
+ num_values += count;
+ }
+
+ accum = accum2;
+ }
+
+ num_bins = bin;
+}
+
+static void fill_bins(int n, double *data)
+{
+ int i;
+
+ for (i = 0; i < n; i++) {
+ int j, bin;
+ struct bin *b;
+
+ // todo nan
+
+ j = get_slot(data[i]);
+
+ if (!slot_bins[j])
+ continue;
+
+ bin = slot_bins[j] - 1;
+ b = &bins[bin];
+
+ values[b->base + b->count++] = data[i];
+ }
+
+ //G_percent(i, n, 2);
+}
+
+static int compare(const void *aa, const void *bb)
+{
+ double a = *(const double *)aa;
+ double b = *(const double *)bb;
+
+ if (a < b)
+ return -1;
+ if (a > b)
+ return 1;
+ return 0;
+}
+
+static void sort_bins(void)
+{
+ int bin;
+
+ for (bin = 0; bin < num_bins; bin++) {
+ struct bin *b = &bins[bin];
+
+ qsort(&values[b->base], b->count, sizeof(double), compare);
+ //G_percent(bin, num_bins, 2);
+ }
+ //G_percent(bin, num_bins, 2);
+}
+
+static void compute_quantiles(int recode, double *quantile, struct write *report)
+{
+ int bin = 0;
+ double prev_v = minimum;
+ int quant;
+
+ for (quant = 0; quant < num_quants; quant++) {
+ struct bin *b;
+ double next = get_quantile(quant);
+ double k, v;
+ int i0, i1;
+
+ for (; bin < num_bins; bin++) {
+ b = &bins[bin];
+ if (b->origin + b->count >= next)
+ break;
+ }
+
+ if (bin < num_bins) {
+ k = next - b->origin;
+ i0 = (int)floor(k);
+ i1 = (int)ceil(k);
+
+ if (i1 > b->count - 1)
+ i1 = b->count - 1;
+
+ v = (i0 == i1)
+ ? values[b->base + i0]
+ : values[b->base + i0] * (i1 - k) +
+ values[b->base + i1] * (k - i0);
+ }
+ else
+ v = maximum;
+
+ *quantile = v;
+ prev_v = v;
+ if (report != NULL && report->write2file == TRUE)
+ fprintf(report->fp, "%f:%f:%f\n", prev_v, v, quants[quant]);
+ }
+}
+
+double quantile(double q, int n, double *data, struct write *report)
+{
+ int i, recode=FALSE;
+ double quantile;
+
+ num_slots = 1000000; // # of slots
+
+ num_quants = 1;
+ quants = (double *) G_malloc(num_quants * sizeof(double));
+ quants[0] = q; // compute quantile for 0.05
+
+ // initialize list of slots
+ slots = (unsigned int *) G_calloc(num_slots, sizeof(unsigned int));
+ slot_bins = (unsigned short *) G_calloc(num_slots, sizeof(unsigned short));
+
+ // get min and max of the data
+ minimum = data[0];
+ maximum = data[0];
+ for (i=0; i<n; i++) {
+ minimum = MIN(data[i], minimum);
+ maximum = MAX(data[i], maximum);
+ }
+
+ slot_size = (maximum - minimum) / num_slots;
+ get_slot_counts(n, data); // # of data in each slot
+
+ bins = (struct bin *) G_calloc(num_quants, sizeof(struct bin));
+ initialize_bins();
+ G_free(slots);
+
+ values = (double *) G_calloc(num_values, sizeof(double));
+ fill_bins(n, data);
+ G_free(slot_bins);
+
+ sort_bins();
+ compute_quantiles(recode, &quantile, report);
+
+ return quantile;
+}
Added: grass-addons/grass7/vector/v.kriging/stat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/stat.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/stat.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,52 @@
+#include "local_proto.h"
+
+void test_normality(int n, double *data, struct write *report)
+{
+ int i;
+ double *elm;
+ double sum_data = 0.;
+ double res;
+ double sum_res2, sum_res3, sum_res4;
+ double tmp;
+ double mean, skewness, kurtosis;
+ double limit5, limit1;
+
+ // compute mean
+ elm = &data[0];
+ for (i=0; i<n; i++) {
+ sum_data += *elm;
+ elm++;
+ }
+ mean = sum_data / n;
+
+ // compute skewness and curtosis
+ elm = &data[0];
+ sum_res2 = sum_res3 = sum_res4 = 0.;
+ for (i=0; i<n; i++) {
+ res = *elm - mean;
+ sum_res2 += SQUARE(res);
+ sum_res3 += POW3(res);
+ sum_res4 += POW4(res);
+ elm++;
+ }
+
+ tmp = sum_res2 / (n-1);
+ skewness = (sum_res3 / (n-1)) / sqrt(POW3(tmp));
+ kurtosis = (sum_res4 / (n-1)) / SQUARE(tmp) - 3.;
+
+ // expected RMSE
+ limit5 = sqrt(6./n);
+ limit1 = sqrt(24./n);
+
+ fprintf(report->fp,"\n");
+ fprintf(report->fp,"Test of normality\n");
+ fprintf(report->fp,"-----------------\n");
+ fprintf(report->fp,"Test statistics | Value\n");
+ fprintf(report->fp,"skewness | %f\n", skewness);
+ fprintf(report->fp,"kurtosis | %f\n", kurtosis);
+ fprintf(report->fp,"\n");
+ fprintf(report->fp,"Level of confidence | Critical value\n");
+ fprintf(report->fp," 0.05 | %f\n", limit5);
+ fprintf(report->fp," 0.01 | %f\n", limit1);
+ fprintf(report->fp,"-----------------\n\n");
+}
Added: grass-addons/grass7/vector/v.kriging/utils.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,852 @@
+#include "local_proto.h"
+
+#define GNUPLOT "gnuplot -persist"
+
+int cmpVals(const void *v1, const void *v2)
+ {
+ int *p1, *p2;
+
+ p1 = (int *)v1;
+ p2 = (int *)v2;
+ if (p1[0] > p2[0])
+ return 1;
+ else if (p1[0] < p2[0])
+ return -1;
+ else
+ return 0;
+ }
+
+void correct_indices(int i3, struct ilist *list, double *r0, struct points *pnts, struct parameters *var_pars)
+{
+ int i;
+ int n = list->n_values;
+ int *vals = list->value;
+ double *pnt = pnts->r;
+ int type = var_pars->type;
+ double sqDist = type == 2 ? SQUARE(var_pars->horizontal.max_dist) : SQUARE(var_pars->max_dist);
+
+ int n_new = 0;
+ int *newvals;
+ newvals = (int *) G_malloc(n * sizeof(int));
+
+ double *r;
+
+ for (i=0; i<n; i++) {
+ *vals -= 1;
+ r = &pnt[3 * *vals];
+
+ if (type != 2 && SQUARE(*r - *r0) + SQUARE(*(r+1) - *(r0+1)) + SQUARE(*(r+2) - *(r0+2)) <= sqDist) {
+ *newvals = *vals;
+ n_new++;
+ newvals++;
+ }
+ vals++;
+ }
+
+ if (type != 2 && n_new < n) {
+ list->n_values = n_new;
+ G_realloc(list->value, n_new);
+ G_realloc(newvals, n_new);
+ memcpy(list->value, newvals, n_new * sizeof(int));
+ }
+}
+
+// make coordinate triples xyz
+void triple(double x, double y, double z, double *triple)
+{
+ double *t;
+ t = &triple[0];
+
+ *t = x;
+ *(t+1) = y;
+ *(t+2) = z;
+}
+
+// compute coordinate differences
+void coord_diff(int i, int j, double *r, double *dr)
+{
+ int k = 3*i, l = 3*j;
+ double *rk, *rl, zi, zj, *drt;
+ rk = &r[k];
+ rl = &r[l];
+ drt = &dr[0];
+
+ if ((*rk == 0.) && (*(rk+1) == 0.) && (*(rk+2) == 0.)) {
+ G_fatal_error(_("Coordinates of point no. %d are zeros."),i);
+ }
+
+ *drt = *rl - *rk; // dx
+ *(drt+1) = *(rl+1) - *(rk+1); // dy
+ *(drt+2) = *(rl+2) - *(rk+2); // dz
+}
+
+// compute horizontal radius from coordinate differences
+double radius_hz_diff(double *dr)
+{
+ double rds;
+ rds = SQUARE(dr[0]) + SQUARE(dr[1]);
+ return rds;
+}
+
+// compute zenith angle
+double zenith_angle(double dr[3])
+{
+ double zenith;
+ zenith = atan2(dr[2], sqrt(SQUARE(dr[0]) + SQUARE(dr[1])));
+
+ return zenith;
+}
+
+// compute size of the lag
+double lag_size(int direction, struct int_par *xD, struct points *pnts, struct parameters *var_pars, struct write *report)
+{
+ // local variables
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz of input points
+ int fction = var_pars->function;
+ int type = var_pars->type;
+ double max_dist = type == 2 ? var_pars->horizontal.max_dist : var_pars->max_dist; // maximum horizontal distance (or vertical for vertical variogram)
+ double max_dist_vert = type == 2 ? var_pars->vertical.max_dist : var_pars->max_dist; // maximum vertical distance (just for bivariate variogram)
+
+
+ int i, j, k; // loop indices
+ int n_vals; // number of the nearest neighbors (NN)
+ int *j_vals; // indices of NN (NOT sorted by distance)
+ struct ilist *list; // list of NN
+
+ double dr[3]; // coordinate differences
+ double d_min = SQUARE(max_dist) + SQUARE(max_dist_vert); // the shortest distance of the 5%-th NN
+ double lagNN; // square root of d_min -> lag
+ double dist2; // square of the distances between search point and NN
+
+ int perc5; // 5% from total number of NN
+ int add_ident; // number of identical points to add to perc5
+
+ for (i=0; i < n; i++) { // for each input point...
+
+ add_ident = 0; // number of points identical to search point equals to zero
+
+ switch (direction) {
+ case 12: // horizontal variogram
+ list = find_NNs_within(2, i, pnts, max_dist, -1);
+ break;
+ case 3: // vertical variogram
+ list = find_NNs_within(1, i, pnts, max_dist, max_dist_vert);
+ break;
+ default: // anisotropic and bivariate variogram
+ list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+ break;
+ }
+ n_vals = list->n_values; // number of overlapping rectangles
+
+ // find number of points identical to the search point
+ if (n_vals > 0) {
+ j_vals = &list->value[0]; // indices of overlapping rectangles (note: increased by 1)
+
+ for (j=0; j < n_vals; j++) { // for each overlapping rectangle:
+ 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.) {
+ add_ident++;
+ }
+ break;
+ case 3: // vertical
+ if (dr[2] == 0.) {
+ add_ident++;
+ }
+ break;
+ default: // aniso / bivar
+ if(dr[0] == 0. && dr[1] == 0. && dr[2] == 0.) {
+ add_ident++;
+ }
+ break;
+ } // end switch: test direction
+ j_vals++; // index of the next NN
+ } // end j for loop
+ } //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);
+ }
+ 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
+
+ coord_diff(i, perc5, pnts->r, dr); // coordinate differences of search point and 5%th NN
+ switch (direction) { // compute squared distance for particular direction
+ case 12: // horizontal
+ dist2 = SQUARE(dr[0]) + SQUARE(dr[1]);
+ break;
+ case 3: // vertical
+ dist2 = SQUARE(dr[2]);
+ break;
+ default: // aniso / bivar
+ dist2 = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
+ break;
+ }
+
+ if (dist2 != 0.) { // and compare it with minimum distance between search point and 5% NN
+ d_min = MIN(d_min, dist2);
+ }
+
+ r += 3; // go to next search point
+ } // end i for loop: distance of i-th search pt and 5%-th NN
+
+ G_free_ilist(list); // free list memory
+ lagNN = sqrt(d_min);
+
+ return lagNN;
+}
+
+double linear(double x, double a, double b)
+{
+ double y;
+ y = a * x + b;
+
+ return y;
+}
+
+double exponential(double x, double nugget, double part_sill, double h_range)
+{
+ double y;
+ y = nugget + part_sill * (1. - exp(- 3. * x/h_range)); // practical
+
+ return y;
+}
+
+double spherical(double x, double a, double b, double c)
+{
+ double y;
+ double ratio;
+ if (x < c) {
+ ratio = x / c;
+ y = a + b * (1.5 * ratio - 0.5 * POW3(ratio));
+ }
+ else
+ y = a + b;
+
+ return y;
+}
+
+double gaussian(double x, double a, double b, double c)
+{
+ double y;
+ double c2 = SQUARE(c);
+ double ratio;
+ ratio = x / c2;
+ y = a + b * (1 - exp(-3. * ratio));
+
+ return y;
+}
+
+// compute # of lags
+int lag_number(double lag, double *varpar_max)
+{
+ int n; // number of lags
+ n = (int) round(*varpar_max / lag); // compute number of lags
+ *varpar_max = n * lag; // set up modified maximum distance (to not have empty lags)
+
+ return n;
+}
+
+// maximal horizontal and vertical distance to restrict variogram computing
+void variogram_restricts(struct int_par *xD, struct points *pnts, struct parameters *var_pars)
+{
+ int n = pnts->n; // # of points
+ double *r = pnts->r; // xyz of points
+ struct write *report = &xD->report;
+
+ int i;
+ double *min, *max; // extend
+ double dr[3]; // coordinate differences
+ double h_maxG; // modify maximum distance (to not have empty lags)
+
+ int dimension; // dimension: hz / vert / aniso
+ char type[12];
+ variogram_type(var_pars->type, type); // set up: 0 - "hz", 1 - "vert", 2 - "bivar", 3 - "aniso"
+
+ // 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
+
+ dr[0] = *max - *min; // x range
+ dr[1] = *(max+1) - *(min+1); // y range
+ dr[2] = *(max+2) - *(min+2); // z range
+
+ switch (var_pars->type) {
+ case 1: // vertical variogram
+ 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)
+ break;
+ }
+
+ if (report->write2file == TRUE) { // report name available:
+ write2file_varSetsIntro(var_pars->type, report); // describe properties
+ }
+
+ // set up code according to type
+ switch (var_pars->type) {
+ case 0: // horizontal
+ dimension = 12;
+ break;
+ case 1: // vertical
+ dimension = 3;
+ break;
+ case 3: // anisotropic
+ dimension = 0;
+ break;
+ }
+
+ if (var_pars->type == 2) { // bivariate variogram
+ // horizontal direction:
+ var_pars->lag = lag_size(12, xD, pnts, var_pars, report); // lag distance
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // number of lags
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+
+ // vertical direction
+ var_pars->lag_vert = lag_size(3, xD, pnts, var_pars, report); // lag distance
+ var_pars->nLag_vert = lag_number(var_pars->lag_vert, &var_pars->max_dist_vert); // # of lags
+ var_pars->max_dist_vert = var_pars->nLag_vert * var_pars->lag_vert; // max distance
+ }
+
+ else { // univariate variograms (hz / vert / aniso)
+ var_pars->lag = lag_size(dimension, xD, pnts, var_pars, report); // lag size
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // # of lags
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ }
+
+ if (report->write2file == TRUE) { // report name available:
+ write2file_varSets(report, var_pars); // describe properties
+ }
+}
+
+void geometric_anisotropy(struct int_par *xD, struct points *pnts)
+{
+ int i;
+ double ratio = xD->aniso_ratio;
+
+ if (ratio <= 0.) { // ratio is negative or zero:
+ if (xD->report.name) { // close report file
+ fprintf(xD->report.fp, "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Anisotropy ratio must be greater than zero..."));
+ }
+
+ double *r;
+ r = &pnts->r[0];
+
+ struct RTree *R_tree; // spatial index for input
+ struct RTree_Rect *rect;
+
+ R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
+
+ for (i=0; i < pnts->n; i++) {
+ *(r+2) = ratio * *(r+2);
+
+ rect = RTreeAllocRect(R_tree);
+ RTreeSetRect3D(rect, R_tree, *r, *r, *(r+1), *(r+1), *(r+2), *(r+2));
+ RTreeInsertRect(rect, i+1, R_tree);
+
+ r += 3;
+ }
+
+ pnts->R_tree = R_tree;
+}
+
+// Least Squares Method
+mat_struct *LSM(mat_struct *A, mat_struct *x)
+{
+ int i, nr, nc;
+ mat_struct *AT, *ATA, *ATA_Inv, *ATx, *T;
+
+ /* A[nr x nc] */
+ nr = A->rows;
+ nc = A->cols;
+
+ /* LMS */
+ AT = G_matrix_transpose(A); /* Transposed design matrix */
+ ATA = G_matrix_product(AT, A); /* AT*A */
+ ATA_Inv = G_matrix_inverse(ATA); /* (AT*A)^-1 */
+ ATx = G_matrix_product(AT, x); /* AT*x */
+ T = G_matrix_product(ATA_Inv, ATx); /* theta = (AT*A)^(-1)*AT*x */
+
+ return T;
+}
+
+// set up type of function according to GUI (user)
+int set_function(char *variogram, struct write *report)
+{
+ int function;
+ if (strcmp(variogram, "linear") == 0) {
+ function = 0;
+ }
+ else if (strcmp(variogram, "parabolic") == 0) {
+ function = 1;
+ }
+ else if (strcmp(variogram, "exponential") == 0) {
+ function = 2;
+ }
+ else if (strcmp(variogram, "spherical") == 0) {
+ function = 3;
+ }
+ else if (strcmp(variogram, "gaussian") == 0) {
+ function = 4;
+ }
+ else if (strcmp(variogram, "bivariate") == 0) {
+ function = 5;
+ }
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Set up correct name of variogram function..."));
+ }
+
+ return function;
+}
+
+// set up terminal and format of output plot
+void set_gnuplot(char *fileformat, struct parameters *var_pars)
+{
+ if (strcmp(fileformat,"cdr") == 0) {
+ strcpy(var_pars->term, "corel");
+ strcpy(var_pars->ext, "cdr");
+ }
+ if (strcmp(fileformat,"dxf") == 0) {
+ strcpy(var_pars->term, "dxf");
+ strcpy(var_pars->ext, "dxf");
+ }
+ if (strcmp(fileformat,"eps") == 0) {
+ strcpy(var_pars->term, "postscript");
+ strcpy(var_pars->ext, "eps");
+ }
+ if (strcmp(fileformat,"pdf") == 0) {
+ strcpy(var_pars->term, "pdfcairo");
+ strcpy(var_pars->ext, "pdf");
+ }
+ if (strcmp(fileformat,"png") == 0) {
+ strcpy(var_pars->term, "png");
+ strcpy(var_pars->ext, "png");
+ }
+ if (strcmp(fileformat,"svg") == 0) {
+ strcpy(var_pars->term, "svg");
+ strcpy(var_pars->ext, "svg");
+ }
+}
+
+// plot experimental variogram
+void plot_experimental_variogram(struct int_par *xD, struct parameters *var_par)
+{
+ int bivar = var_par->type == 2 ? TRUE : FALSE;
+
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert;
+ mat_struct *gamma_M = var_par->gamma; // pointer to gamma matrix
+ FILE *gp; // pointer to file
+
+ nr = gamma_M->rows; // # of rows of gamma matrix
+ nc = gamma_M->cols; // # of cols of gamma matrix
+
+ double *gamma;
+ gamma = &gamma_M->vals[0]; // values of gamma matrix
+
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (access("dataE.dat", W_OK) < 0) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
+ }
+
+ /* 3D experimental variogram */
+ if (bivar == TRUE) {
+ vert = &var_par->vertical.h[0];
+ for (i=0; i < nc+1; i++) { // for each row (nZ)
+ h = &var_par->horizontal.h[0];
+ for (j=0; j < nr+1; j++) { // for each col (nL)
+ if (i == 0 && j == 0) {
+ fprintf(gp, "%d ", nr); // write to file count of columns
+ }
+ else if (i == 0 && j != 0) { // 1st row
+ fprintf(gp, "%f", *h); // write h to 1st row of the file
+ h++;
+ if (j < nr) {
+ fprintf(gp, " "); // spaces between elements
+ }
+ } // end 1st row setting
+
+ else {
+ if (j == 0 && i != 0) // 1st column
+ fprintf(gp, "%f ", *vert);// write vert to 1st column of the file
+ else {
+ if (isnan(*gamma)) {
+ fprintf(gp, "NaN"); // write other elements
+ }
+ else {
+ fprintf(gp, "%f", *gamma); // write other elements
+ }
+
+ if (j != nr) {
+ fprintf(gp, " "); // spaces between elements
+ }
+ gamma++;
+ }
+ } // end columns settings
+ } // end j
+ fprintf(gp, "\n");
+ if (i!=0) { // vert must not increase in the 1st loop
+ vert++;
+ }
+ } // end i
+ } // end 3D
+
+ /* 2D experimental variogram */
+ else {
+ h = &var_par->h[0];
+ for (i=0; i < nr; ++i) {
+ if (isnan(*gamma)) {
+ fprintf(gp, "%f NaN\n", *h);
+ }
+ else {
+ fprintf(gp, "%f %f\n", *h, *gamma);
+ }
+ h++;
+ gamma++;
+ }
+ }
+ fclose(gp);
+
+ gp = popen(GNUPLOT, "w"); /* open file to create plots */
+ if (gp == NULL)
+ G_message("Unable to plot variogram");
+
+ fprintf(gp, "set terminal wxt %d size 800,450\n", var_par->type);
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp, "set title \"Experimental variogram (3D) of <%s>\"\n", var_par->name);
+ fprintf(gp, "set xlabel \"h interval No\"\n");
+ fprintf(gp, "set ylabel \"dz interval No\"\n");
+ fprintf(gp, "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
+ fprintf(gp, "set key below\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
+ fprintf(gp, "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\"\n");
+ }
+
+ else { // univariate variogram
+ char dim[6], dist[2];
+ if (xD->i3 == TRUE) { // 3D
+ if (var_par->type == 0) {
+ strcpy(dim,"hz");
+ }
+ if (var_par->type == 1) {
+ strcpy(dim, "vert");
+ }
+ if (var_par->type == 3) {
+ strcpy(dim, "aniso");
+ }
+ fprintf(gp, "set title \"Experimental variogram (3D %s) of <%s>\"\n", dim, var_par->name);
+ }
+ else { // 2D
+ fprintf(gp, "set title \"Experimental variogram (2D) of <%s>\"\n", var_par->name);
+ }
+
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5\n");
+ }
+ fclose(gp);
+}
+
+// plot experimental and theoretical variogram
+void plot_var(int i3, int bivar, struct parameters *var_pars)
+{
+ int function, func_h, func_v;
+ double nugget, nugget_h, nugget_v;
+ double sill, sill_h, sill_v;
+ double h_range, h_range_h, h_range_v;
+ double *T;
+
+ // setup theoretical variogram parameters
+ switch (bivar) {
+ case FALSE: // univariate variogram
+ function = var_pars->function;
+ if (function > 0) { // nonlinear variogram
+ nugget = var_pars->nugget;
+ sill = var_pars->sill - nugget;
+ h_range = var_pars->h_range;
+ }
+ else {
+ T = &var_pars->T->vals[0];
+ }
+ break;
+
+ case TRUE: // bivariate variogram
+ if (var_pars->function != 5) {
+ func_h = var_pars->horizontal.function;
+ func_v = var_pars->vertical.function;
+
+ nugget_h = var_pars->horizontal.nugget;
+ sill_h = var_pars->horizontal.sill - nugget_h;
+ h_range_h = var_pars->horizontal.h_range;
+
+ nugget_v = var_pars->vertical.nugget;
+ sill_v = var_pars->vertical.sill - nugget_v;
+ h_range_v = var_pars->vertical.h_range;
+ }
+ else {
+ T = &var_pars->T->vals[0];
+ }
+ break;
+ }
+
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ int type = var_pars->type;
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert; // pointer to vertical bins
+ double h_ratio;
+ double g_teor;
+ double *gamma_var; // pointer to gamma matrix
+ FILE *gp; // pointer to file
+
+ nr = var_pars->gamma->rows; // # of rows of gamma matrix
+ nc = var_pars->gamma->cols; // # of cols of gamma matrix
+
+ gamma_var = &var_pars->gamma->vals[0]; // values of gamma matrix
+
+ if (type != 2) { // univariate variogram
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (access("dataE.dat", W_OK) < 0) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
+ }
+ if (gp == NULL) {
+ G_fatal_error(_("Error writing file"));
+ }
+
+ h = &var_pars->h[0];
+ for (i=0; i < nr; i++) {
+ if (isnan(*gamma_var)) {
+ fprintf(gp, "%f NaN\n", *h); // write other elements
+ }
+ else {
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ }
+ h++;
+ gamma_var++;
+ }
+ }
+
+ else { // bivariate variogram
+ gp = fopen("dataE.dat","r"); // open file to read experimental variogram
+ if (access("dataE.dat", F_OK) < 0) {
+ G_fatal_error(_("You have probably deleted dataE.dat - process middle phase again, please."));
+ }
+ }
+
+ if (fclose(gp) != 0) {
+ 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) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
+ }
+
+ double dd;
+ double hh;
+ double hv[2];
+ h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
+
+ // Compute theoretical variogram
+ switch (bivar) {
+ case 0: // Univariate variogram:
+ for (i=0; i <= nr; i++) {
+ hh = i==0 ? 0. : *h;
+
+ switch (function) {
+ case 0: // linear
+ dd = *T * hh + *(T+1);
+ break;
+ case 1: // parabolic
+ dd = *T * SQUARE(hh) + *(T+1);
+ break;
+ case 2: // exponential
+ dd = nugget + sill * (1 - exp(- 3. * hh/h_range)); // practical
+ break;
+ case 3: // spherical
+ if (hh < h_range) {
+ h_ratio = hh / h_range;
+ dd = nugget + sill * (1.5 * h_ratio - 0.5 * POW3(h_ratio));
+ }
+ else {
+ dd = sill + nugget;
+ }
+ break;
+ case 4: // Gaussian
+ h_ratio = SQUARE(hh) / SQUARE(h_range);
+ dd = nugget + sill * (1 - exp(-3. * h_ratio));
+ break;
+ }
+ fprintf(gp, "%f %f\n", hh, dd);
+
+ if (i > 0) {
+ h++;
+ }
+ } // end i for cycle
+ break;
+
+ case 1: // bivariate (3D)
+ // Coefficients of theoretical variogram
+ T = &var_pars->T->vals[0]; // values
+ vert = &var_pars->vertical.h[0];
+ for (i=0; i < nc+1; i++) { // for each row
+ h = &var_pars->horizontal.h[0];
+ for (j=0; j < nr+1; j++) { // for each col
+ if (i == 0 && j == 0) { // the 0-th element...
+ fprintf(gp, "%d ", nr); // ... is number of columns
+ }
+ else if (i == 0 && j != 0) { // for the other elements in 1st row...
+ fprintf(gp, "%f", *h); // ... write h
+ if (j < nr) {
+ fprintf(gp, " "); // ... write a space
+ }
+ h++; // go to next h
+ }
+ else { // other rows
+ if (j == 0 && i != 0) {// write vert to 1st column
+ fprintf(gp, "%f ", *vert);
+ }
+ else { // fill the other columns with g_teor
+ hv[0] = *h;
+ hv[1] = *vert;
+ g_teor = variogram_fction(var_pars, hv);
+ fprintf(gp, "%f", g_teor);
+ if (j != nr) {
+ fprintf(gp, " ");
+ }
+ h++;
+ } // end else: fill the other columns
+ } // end fill the other rows
+ } //end j
+ if (i != 0 && j != 0) {
+ vert++;
+ }
+ fprintf(gp, "\n");
+ } //end i
+ break;
+ }
+
+ if (fclose(gp) == EOF) {
+ G_fatal_error(_("Error closing file..."));
+ }
+
+ gp = popen(GNUPLOT, "w"); // open file to create plots
+ if (gp == NULL) {
+ G_message(_("Unable to plot variogram"));
+ }
+
+ if (strcmp(var_pars->term," ") != 0) {
+ fprintf(gp, "set terminal %s size 750,450\n", var_pars->term);
+
+ switch (var_pars->type) {
+ case 0: // horizontal component of bivariate variogram
+ if (var_pars->function == 5) { // component of 3D bivariate
+ fprintf(gp, "set output \"variogram_bivar.%s\"\n", var_pars->ext);
+ }
+ else { // horizontal variogram
+ fprintf(gp, "set output \"variogram_hz.%s\"\n", var_pars->ext);
+ }
+ break;
+ case 1: // vertical variogram
+ fprintf(gp, "set output \"variogram_vertical.%s\"\n", var_pars->ext);
+ break;
+ case 2: // bivariate variogram
+ if (bivar == TRUE) {
+ fprintf(gp, "set output \"variogram_bivariate.%s\"\n", var_pars->ext);
+ }
+ break;
+ case 3: // anisotropic variogram
+ fprintf(gp, "set output \"variogram_final.%s\"\n", var_pars->ext);
+ break;
+ }
+ }
+ else {
+ G_warning("\nVariogram output> If you wish some special file format, please set it at the beginning.\n");
+ fprintf(gp, "set terminal wxt %d size 800,450\n", var_pars->type);
+ }
+
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_pars->name);
+
+ fprintf(gp, "set xlabel \"h interval No\"\n");
+ fprintf(gp, "set ylabel \"dz interval No\"\n");
+ fprintf(gp, "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
+ fprintf(gp, "set key below\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
+ fprintf(gp, "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\", 'dataT.dat' every ::1:1 matrix with lines title \"theoretical variogram\" palette\n");
+ }
+
+ else { // univariate variogram
+ if (i3 == TRUE) { // 3D
+ if (var_pars->type == 0) { // horizontal variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D hz) of <%s>\"\n", var_pars->name);
+ //fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ }
+ else if (var_pars->type == 1) { // vertical variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D vert) of <%s>\"\n", var_pars->name);
+ //fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_par->T->vals[0], var_par->T->vals[1]);
+ }
+ else if (var_pars->type == 3) // anisotropic variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_pars->name);
+ }
+
+ else { // 2D
+ fprintf(gp, "set title \"Experimental and theoretical variogram (2D) of <%s>\"\n", var_pars->name);
+ }
+
+ if (var_pars->type == 2) {
+ switch (var_pars->function) {
+ case 0: // linear
+ fprintf(gp, "set label \"linear: gamma(h) = %f*h + %f\" at screen 0.30,0.90\n", var_pars->T->vals[0], var_pars->T->vals[1]);
+ break;
+ case 1: // parabolic
+ fprintf(gp, "set label \"parabolic: gamma(h) = %f*h^2 + %f\" at screen 0.25,0.90\n", var_pars->T->vals[0], var_pars->T->vals[1]);
+ break;
+ case 2: // exponential
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp, "set label \"exponential: gamma(h) = %f + %f * (1 - exp(-3*h / %f))\" at screen 0.10,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range);
+ break;
+ case 3: // spherical
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp, "set label \"spherical: gamma(h) = %f+%f*(1.5*h/%f-0.5*(h/%f)^3)\" at screen 0.05,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range, var_pars->h_range);
+ break;
+ case 4: // gaussian
+ fprintf(gp, "set label \"gaussian: gamma(h) = %f + %f * (1 - exp(-3*(h / %f)^2))\" at screen 0.10,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range);
+ break;
+ }
+ }
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5, 'dataT.dat' using 1:2 title \"theoretical variogram\" with linespoints\n");
+ }
+ fclose(gp);
+
+ remove("dataE.dat");
+ remove("dataT.dat");
+}
Added: grass-addons/grass7/vector/v.kriging/utils_kriging.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_kriging.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_kriging.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,686 @@
+#include "local_proto.h"
+
+ __inline double square(double x)
+ {
+ return x*x;
+ }
+
+void linear_variogram(struct parameters *var_par, struct write *report)
+{
+ int nZ, nL; // # of gamma matrix rows and columns
+ int nr, nc; // # of plan matrix rows and columns
+ int i, j; // indices
+ double *h, *vert, *gamma;
+ mat_struct *gR;
+
+ nL = var_par->gamma->rows; // # of cols (horizontal bins)
+ nZ = var_par->gamma->cols; // # of rows (vertical bins)
+ gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
+
+ // Test length of design matrix
+ nr = 0; // # of rows of design matrix A - depend on # of non-nan gammas
+ for (i = 0; i < nZ; i++) {
+ for (j = 0; j < nL; j++) {
+ if (!isnan(*gamma)) { // todo: upgrade: !isnan(*c) L434 too
+ nr++;
+ }
+ gamma++;
+ } // end j
+ } // end i
+
+ // Number of columns of design matrix A
+ nc = var_par->function == 5 ? 3 : 2;
+
+ var_par->A = G_matrix_init(nr, nc, nr); // initialise design matrix
+ gR = G_matrix_init(nr, 1, nr); // initialise vector of observations
+
+ // Set design matrix A
+ nr = 0; // index of matrix row
+ gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
+
+ if (var_par->function == 5) { // in case of bivariate variogram
+ vert = &var_par->vertical.h[0]; // vertical direction
+ }
+
+ for (i = 0; i < nZ; i++) {
+ h = var_par->function == 5 ? &var_par->horizontal.h[0] : &var_par->h[0];
+ for (j = 0; j < nL; j++) {
+ if (!isnan(*gamma)) { // write to matrix - each valid element of gamma
+ switch (var_par->function) { // function of theoretical variogram
+ case 5: // bivariate planar
+ G_matrix_set_element(var_par->A, nr, 0, *h);
+ G_matrix_set_element(var_par->A, nr, 1, *vert);
+ G_matrix_set_element(var_par->A, nr, 2, 1.);
+ G_matrix_set_element(gR, nr, 0, *gamma);
+ break;
+ default:
+ G_matrix_set_element(var_par->A, nr, 0, *h);
+ G_matrix_set_element(var_par->A, nr, 1, 1.);
+ G_matrix_set_element(gR, nr, 0, *gamma);
+ break;
+ } // end switch variogram fuction
+ nr++; // length of vector of valid elements (not null)
+ } // end test if !isnan(*gamma)
+ h++;
+ gamma++;
+ } // end j
+ if (var_par->function == 5) {
+ vert++;
+ }
+ } // end i
+ end_loop:
+
+ // Estimate theoretical variogram coefficients
+ var_par->T = LSM(var_par->A, gR); // Least Square Method
+
+ // Test of theoretical variogram estimation
+ if (var_par->T->vals == NULL) { // NULL values of theoretical variogram
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error("Unable to compute 3D theoretical variogram...");
+ } // error
+
+ // constant raster
+ if (var_par->T->vals[0] == 0. && var_par->T->vals[1] == 0.) { //to do: otestuj pre 2d
+ var_par->const_val = 1;
+ G_message(_("Input values to be interpolated are constant."));
+ } // todo: cnstant raster for exponential etc.
+
+ // 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]);
+ }
+}
+
+double bivar_sill(int direction, mat_struct *gamma)
+{
+ int i, j, k;
+ int n; // # of lags
+ int n_gamma = 0; // # of real gammas
+ double gamma_i; // each real gamma
+ double sum_gamma = 0.; // sum of real gammas
+ double sill; // variogram sill
+
+ switch (direction) { // direction is:
+ case 12: // horizontal
+ n = gamma->rows; // use 1st column of experimental variogram matrix
+ break;
+ case 3: // vertical
+ n = gamma->cols; // use 1st row of experimental variogram matrix
+ break;
+ }
+
+ for (i=0; i<n; i++) {
+ switch (direction) { // direction is:
+ case 12: // horizontal
+ j = i; // row index is increasing
+ k = 0; // column index is constant
+ break;
+ case 3: // vertical
+ j = 0; // row index is constant
+ k = i; // column index is increasing
+ break;
+ }
+
+ gamma_i = G_matrix_get_element(gamma, j, k);
+ if (!isnan(gamma_i)) { // gamma is real:
+ sum_gamma += gamma_i; // sum all real elements of the matrix
+ n_gamma++; // count them
+ } // end if
+ } // end for
+
+ sill = sum_gamma / n_gamma;
+
+ return sill;
+}
+
+// Compute sill
+void sill(struct parameters *var_par)
+{
+ // Local variables
+ int type = var_par->type;
+ mat_struct *gamma = var_par->gamma;
+ int nrows = gamma->rows;
+ int ncols = gamma->cols;
+
+ char var_type[12];
+
+ switch (type) {
+ case 2: // bivariate variogram
+ var_par->horizontal.sill = bivar_sill(12, gamma);
+ var_par->vertical.sill = bivar_sill(3, gamma);
+ break;
+
+ default: // hz / vert / aniso variogram
+ var_par->sill = var_par->gamma_sum / var_par->gamma_n; // mean of gamma elements
+
+ variogram_type(var_par->type, var_type); // describe code by string
+ G_message(_("Default sill of %s variogram: %f"), var_type, var_par->sill);
+ break;
+ } // end switch: variogram type
+}
+
+// compare sills
+int sill_compare(struct int_par *xD, struct flgs *flg, struct var_par *var_par, struct points *pnts)
+{
+ // local variables
+ double sill_hz = var_par->hz.sill; // sill in horizontal direction
+ double sill_vert = var_par->vert.sill; // sill in vertical direction
+
+ double diff_sill_05, diff_sill;
+
+ diff_sill_05 = sill_hz > sill_vert ? 0.05 * sill_hz : 0.05 * sill_vert; // critical value as 5% from bigger sill
+ diff_sill = fabsf(sill_hz - sill_vert); // difference between the sills
+
+ if (xD->bivar == TRUE || (!flg->univariate->answer && diff_sill > diff_sill_05)) { // zonal anisotropy
+ var_par->fin.type = 2; // code for bivariate variogram
+ var_par->fin.td = var_par->hz.td; // azimuth limit
+ var_par->fin.horizontal.max_dist = var_par->hz.max_dist; // maximum dist in hz direction
+ var_par->fin.horizontal.nLag = var_par->hz.nLag;
+ var_par->fin.horizontal.lag = var_par->hz.lag;
+
+ var_par->fin.vertical.max_dist = var_par->vert.max_dist; // maximum dist in vert direction
+ var_par->fin.vertical.nLag = var_par->vert.nLag;
+ var_par->fin.vertical.lag = var_par->vert.lag;
+ }
+
+ 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.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
+ }
+}
+
+// formulation of variogram functions
+double variogram_fction(struct parameters *var_par, double *dr)
+{
+ // Local variables
+ int i;
+ int variogram = var_par->function; // function of theoretical variogram
+ int type = var_par->type; // hz / vert / bivar / aniso
+ double nugget;
+ double part_sill;
+ double h_range;
+ double *T;
+
+ if (type == 2 || var_par->function == 0) { // bivariate variogram:
+ T = &var_par->T->vals[0]; // coefficients of theoretical variogram
+ }
+
+ double radius; // square distance of the point couple
+ double h;
+ double vert;
+ double h_ratio;
+ double teor_var, result = 0.; // estimated value of the theoretical variogram
+
+ int n_cycles = (type == 2 && var_par->function != 5) ? 2 : 1; // # of cycles (bivariate (not linear) or univariate)
+
+ for (i=0; i < n_cycles; i++) {
+ if (n_cycles == 2) { // bivariate variogram
+ variogram = i==0 ? var_par->horizontal.function : var_par->vertical.function;
+ h = i==0 ? dr[0] : dr[1];
+ radius = SQUARE(h);
+ nugget = i==0 ? var_par->horizontal.nugget : var_par->vertical.nugget;
+ part_sill = i==0 ? var_par->horizontal.sill - nugget : var_par->vertical.sill - nugget;
+ h_range = i==0 ? var_par->horizontal.h_range : var_par->vertical.h_range;
+ }
+ else { // univariate variogram or linear bivariate
+ variogram = var_par->function;
+ if (variogram == 5) {
+ h = dr[0];
+ vert = dr[1];
+ }
+ else { // univariate variogram
+ radius = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]); // weighted dz
+ h = sqrt(radius);
+
+ nugget = var_par->nugget;
+ part_sill = var_par->sill - nugget;
+ h_range = var_par->h_range;
+ }
+ }
+
+ switch (variogram) {
+ case 0: // linear function
+ teor_var = linear(h, *T, *(T+1));
+ break;
+ case 1: // parabolic function TODO: delete
+ teor_var = *T * radius + *(T+1);
+ break;
+ case 2: // exponential function
+ teor_var = exponential(h, nugget, part_sill, h_range);
+ break;
+ case 3: // spherical function
+ teor_var = spherical(h, nugget, part_sill, h_range);
+ break;
+ case 4: // Gaussian function
+ teor_var = gaussian(radius, nugget, part_sill, h_range);
+ break;
+ case 5:
+ teor_var = *T * h + *(T+1) * vert + *(T+2);
+ break;
+ } // end switch (variogram)
+
+ result += teor_var;
+ }
+ if (type == 2 && var_par->function != 5) {
+ result -= 0.5 * (var_par->horizontal.sill - nugget + var_par->vertical.sill - nugget);
+ }
+
+ return result;
+}
+
+// compute coordinates of reference point - cell centre
+void cell_centre(unsigned int col, unsigned int row, unsigned int dep, struct int_par *xD, struct reg_par *reg, double *r0, struct parameters *var_par)
+{
+ // Local variables
+ int i3 = xD->i3;
+
+ r0[0] = reg->west + (col + 0.5) * reg->ew_res; // x0
+ r0[1] = reg->north - (row + 0.5) * reg->ns_res; // y0
+
+ if (i3 == TRUE) { // 3D interpolation
+ switch (var_par->function) {
+ case 5: // bivariate variogram
+ r0[2] = reg->bot + (dep + 0.5) * reg->bt_res; // z0
+ break;
+ default: // univariate (anisotropic) function
+ r0[2] = xD->aniso_ratio * (reg->bot + (dep + 0.5) * reg->bt_res); // z0
+ break;
+ } // end switch
+ } // end if
+
+ else { // 2D interpolation
+ r0[2] = 0.;
+ }
+}
+
+// set up elements of G matrix
+mat_struct *set_up_G(struct points *pnts, struct parameters *var_par, struct write *report)
+{
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz coordinates of input points
+ int type = var_par->type; // hz / vert / aniso / bivar
+
+ int i, j; // indices of matrix rows/cols
+ int n1 = n+1; // # of matrix rows/cols
+ double theor_var; // GM element = theor_var(distance)
+ double *dr; // dx, dy, dz between point couples
+ mat_struct *GM; // GM matrix
+
+ dr = (double *) G_malloc(3 * sizeof(double));
+ GM = G_matrix_init(n1, n1, n1); // G[n1][n1] matrix
+
+ doublereal *md, *mu, *ml, *dbu, *dbl, *m1r, *m1c;
+ dbu = &GM->vals[0]; // upper matrix elements
+ dbl = &GM->vals[0]; // lower matrix elements
+ md = &GM->vals[0]; // diagonal matrix elements
+ m1r = &GM->vals[n1-1]; // GM - last matrix row
+ m1c = &GM->vals[n1*n]; // GM - last matrix col
+
+ // G[n;n]
+ for (i = 0; i < n1; i++) { // for elements in each row
+ dbu += n1; // new row of the U
+ dbl++; // new col of the L
+ mu = dbu; // 1st element in the new U row
+ ml = dbl; // 1st element in the new L col
+ for (j = i; j < n1; j++) { // elements in cols of upper matrix
+ if (i != j && i != n && j != n) { // non-diagonal elements except last row/col
+ coord_diff(i, j, r, dr); // compute coordinate differences
+ if (type == 2) { // bivariate variogram
+ dr[0] = sqrt(radius_hz_diff(dr));
+ dr[1] = dr[2];
+ }
+ theor_var = variogram_fction(var_par, dr); // compute GM element
+
+ if (isnan(theor_var)) { // not a number:
+ if (report->name) {
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Theoretical variogram is NAN..."));
+ }
+
+ *mu = *ml = (doublereal) theor_var; // set the value to the matrix
+ mu += n1; // go to next element in the U row
+ ml++; // go to next element in the L col
+ } // end non-diagonal elements condition
+ } // end j loop
+
+ // go to the diagonal element in the next row
+ dbu++; // U
+ dbl += n1; // L
+ *md = 0.0; // set diagonal
+ md += n1+1;
+
+ // Set ones
+ if (i<n1-1) { // all elements except the last one...
+ *m1r = *m1c = 1.0; // ... shall be 1
+ m1r += n1; // go to next col in last row
+ m1c++; // go to next row in last col
+ } // end "last 1" condition
+ } // end i loop
+
+ free(dr);
+ return GM;
+}
+
+// make G submatrix for rellevant points
+mat_struct *submatrix(struct ilist *index, mat_struct *GM_all, struct write *report)
+{
+ // Local variables
+ int n = index->n_values; // # of selected points
+ mat_struct *GM = GM_all; // whole G matrix
+
+ int i, j, k, N1 = GM->rows, n1 = n+1, *dinR, *dini, *dinj;
+ doublereal *dbo, *dbx, *dbu, *dbl, *md, *mu, *ml, *m1r, *m1c;
+
+ mat_struct *sub; // new submatrix
+ sub = G_matrix_init(n1, n1, n1);
+
+ if (sub == NULL) {
+ if (report->name) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to initialize G-submatrix..."));
+ }
+
+ // Initialize indexes: GM[0;1]
+ // - cols
+ dini = &index->value[0]; // i - set processing column
+
+ // initialize new col
+ dinR = &index->value[1]; // return to first processing row - lower GM
+
+ dbo = &GM->vals[*dini * N1 + *dinR]; // 1st value in GM_all
+ md = &sub->vals[0]; // GM - diagonal
+ dbu = &sub->vals[n1]; // GM - upper
+ dbl = &sub->vals[1]; // GM - lower
+ m1r = &sub->vals[n1-1]; // GM - last row
+ m1c = &sub->vals[n1*n]; // GM - last col
+
+ for (i=0; i<=n; i++) { // set up cols of submatrix
+ // original matrix
+ dbx = dbo; // new col in GM_all
+ dinj = dinR; // orig: inicialize element (row) in (col) - upper matrix
+ // submatrix
+ mu = dbu; // sub: start new column
+ ml = dbl; // sub: start new row
+ for (j=i+1; j<n; j++) { // for rows
+ //submatrix
+ *mu = *ml = *dbx; // general elements
+ mu += n1; // go to element in next column
+ ml++; // go to element in next row
+ // Original matrix
+ dinj++; // set next ind
+ dbx += *dinj - *(dinj-1);
+ } // end j
+
+ // Original matrix
+ dini++; // set next ind
+
+ dbu += n1+1; // new col in GM
+ dbl += n1+1; // new row in GM
+
+ dinR++;
+
+ dbo += (*dini - *(dini-1)) * N1 + (*dinR - *(dinR-1));
+ // Set ones
+ *m1r = *m1c = 1.0;
+ m1r += n1;
+ m1c++;
+ // Set diagonals
+ *md = 0.0;
+ md += n1+1;
+ } // end i
+
+ return sub;
+}
+
+// set up g0 vector
+mat_struct *set_up_g0(struct int_par *xD, struct points *pnts, struct ilist *index, double *r0, struct parameters *var_par)
+{
+ // Local variables
+ int i3 = xD->i3; // interpolation: 2D or 3D
+ int bivar = xD->bivar; // variogram: uni- or bivariate
+ int type = var_par->type; // variogram: hz / vert / aniso / bivar
+ double *r; // xyz coordinates of input points
+
+ int n; // # of points (all or selected)
+ int n_ind = index->n_values; // # of selected points
+ int *lind; // pointer to indices of selected points
+
+ if (n_ind == 0) { // no selected points:
+ n = pnts->n; // use all points
+ r = &pnts->r[0];
+ }
+ else { // any selected points:
+ n = n_ind; // reduce # of input points
+ lind = &index->value[0];
+ r = &pnts->r[*lind * 3];
+ }
+
+ int n1 = n + 1;
+
+ int i; // index of elements and length of g0 vector
+ double teor_var; // estimated value based on theoretical variogram and distance of the input points
+ double *dr; // coordinate differences dx, dy, dz between couples of points
+ mat_struct *g0; // vector of estimated differences between known and unknown values
+
+ dr = (double *) G_malloc(3 * sizeof(double)); // Coordinate differences
+ g0 = G_matrix_init(n1,1,n1);
+
+ doublereal *g = g0->vals;
+
+ for (i=0; i<n; i++) { // count of input points
+ // Coord diffs (input points and cell/voxel center)
+ dr[0] = *r - *r0; // dx
+ dr[1] = *(r+1) - *(r0+1); // dy
+ switch (i3) {
+ case FALSE:
+ // Cell value diffs estimation using linear variogram
+ dr[2] = 0.0; // !!! optimalize - not to be necessary to use this line
+ break;
+ case TRUE:
+ dr[2] = *(r+2) - *(r0+2); // dz
+ break;
+ }
+
+ if (type == 2) { // bivariate variogram
+ dr[0] = sqrt(radius_hz_diff(dr));
+ dr[1] = dr[2];
+ }
+
+ *g = variogram_fction(var_par, dr); // compute GM element using theoretical variogram function
+ g++;
+
+ if (n_ind == 0) { // no selected points:
+ r += 3; // use each point
+ }
+ else { // selected points:
+ lind++; // go to the next index
+ r += 3 * (*lind - *(lind-1)); // use next selected point
+ }
+ } // end i loop
+
+ // condition: sum of weights = 1
+ *g = 1.0; // g0 [n1x1]
+
+ return g0;
+}
+
+// compute result to be located on reference point
+double result(struct points *pnts, struct ilist *index, mat_struct *w0)
+{
+ // local variables
+ int n; // # of points
+ double *vo; // pointer to input values
+
+ int n_ind = index->n_values; // # of selected points
+ int *indo; // pointer to indices of selected pt
+
+ // Setup n depending on NN points selection
+ n = n_ind == 0 ? pnts->n : n_ind;
+
+ int i;
+ mat_struct *ins, *w, *rslt_OK;
+ doublereal *vt, *wo, *wt;
+
+ ins = G_matrix_init(n, 1, n); // matrix of selected vals
+ w = G_matrix_init(1, n, 1); // matrix of selected weights
+
+ if (n_ind == 0) { // there are no selected points:
+ vo = &pnts->invals[0]; // use all input values
+ }
+ else { // there are selected points:
+ indo = &index->value[0]; // original indexes
+ vo = &pnts->invals[*indo]; // original input values
+ }
+
+ wo = &w0->vals[0]; // pointer to weights of input vals
+ vt = &ins->vals[0]; // pointer to selected inputs values
+ wt = &w->vals[0]; // pointer to weights without the last one
+
+ for (i=0; i<n; i++) { // for each input point:
+ *vt = *vo; // subval = selected original
+ *wt = *wo; // sub weight = selected original
+
+ // go to the next:
+ if (n_ind == 0) { // use all points:
+ vo++; // input value
+ }
+ else { // use selected points:
+ indo++; // index of selected point
+ vo += *indo - *(indo-1); // selected input value
+ }
+ vt++; // element of value matrix
+ wo++; // weight of the value
+ wt++; // element of weight matrix
+ } // end i for loop
+
+ rslt_OK = G_matrix_product(w,ins); // interpolated value
+
+ G_matrix_free(w);
+ G_matrix_free(ins);
+
+ return rslt_OK->vals[0];
+}
+
+// validation
+void crossvalidation(struct int_par *xD, struct points *pnts, struct parameters *var_par)
+{
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // coordinates of points
+ int i3 = xD->i3;
+ double *vals = pnts->invals; // input values
+
+ int type = var_par->type;
+ double ratio = type == 3 ? xD->aniso_ratio : 1.; // anisotropic ratio
+ mat_struct *GM = var_par->GM; // GM = theor_var(distances)
+
+ int i, j;
+ int n_vals;
+ int n1;
+ double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
+ double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
+
+ struct ilist *list;
+ struct RTree_Rect *r_cell;
+ mat_struct *GM_sub;
+ mat_struct *GM_Inv, *g0, *w0;
+ double rslt_OK;
+
+ struct write *crossvalid = &xD->crossvalid;
+ struct write *report = &xD->report;
+ double *normal, *absval, *norm, *av;
+
+ normal = (double *) G_malloc(n * sizeof(double));
+ absval = (double *) G_malloc(n * sizeof(double));
+ norm = &normal[0];
+ av = &absval[0];
+
+ FILE *fp;
+ if (crossvalid->name) { // CV report file available:
+ fp = fopen(crossvalid->name, "w");
+ }
+
+ for (i=0; i<n; i++) { // for each input point [r0]:
+ list = G_new_ilist(); // create list of overlapping rectangles
+
+ if (i3 == TRUE) {
+ list = find_NNs_within(3, i, pnts, max_dist, max_dist_vert);
+ }
+ else {
+ list = find_NNs_within(2, i, pnts, max_dist, max_dist_vert);
+ }
+
+ n_vals = list->n_values; // # of overlapping rectangles
+
+ if (n_vals > 0 ) { // if positive:
+ correct_indices(i3, 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);
+
+ 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
+
+ G_matrix_free(g0);
+ G_matrix_free(GM_Inv);
+
+ rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
+
+ //Create output
+ *norm = *vals - rslt_OK; // differences between input and interpolated values
+ *av = fabsf(*norm); // 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);
+ }
+ else { // 2D interpolation:
+ fprintf(fp,"%d %.3f %.3f %f %f %f\n", i, *r, *(r+1), *vals, rslt_OK, *norm);
+ }
+ } // end if: n_vals > 0
+
+ else { // n_vals == 0:
+ fprintf(fp,"%d. point does not have neighbours in given radius\n", i);
+ }
+
+ r += 3; // next rectangle
+ vals++;
+ norm++;
+ av++;
+
+ G_free_ilist(list); // free list memory
+ } // end i for loop
+
+ fclose(fp);
+ G_message(_("Cross validation results have been written into <%s>"), crossvalid->name);
+
+ if (report->name) {
+ double quant05, quant10, quant25, quant50, quant75, quant90, quant95;
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Cross validation results ***\n");
+
+ test_normality(n, normal, report);
+
+ fprintf(report->fp, "Quantile of absolute values\n");
+ quant95 = quantile(0.95, n, absval, report);
+ //quant90 = quantile(0.90, n, absval, report);
+ //quant75 = quantile(0.75, n, absval, report);
+ //quant50 = quantile(0.50, n, absval, report);
+ //quant25 = quantile(0.25, n, absval, report);
+ //quant10 = quantile(0.10, n, absval, report);
+ //quant05 = quantile(0.05, n, absval, report);
+ }
+}
Added: grass-addons/grass7/vector/v.kriging/utils_raster.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_raster.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_raster.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,60 @@
+#include "local_proto.h"
+
+/* open raster layer */
+void open_layer(struct int_par *xD, struct reg_par *reg, struct output *out)
+{
+ struct write *report = &xD->report;
+ /* 2D Raster layer */
+ if (xD->i3 == FALSE) {
+ 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);
+ }
+ 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);
+ }
+ 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)
+{
+ // Local variables
+ int i3 = xD->i3;
+ struct write *report = &xD->report;
+
+ int pass = 0; /* Number of processed cells */
+ 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;
+
+ 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);
+ }
+ G_fatal_error(_("Error writing cell (%d,%d,%d) with value %f, nrows = %d"), row, col, dep, rslt_OK, reg->nrows);
+ }
+ pass++;
+ break;
+ }
+ return pass;
+}
+
Added: grass-addons/grass7/vector/v.kriging/utils_write.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_write.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_write.c 2015-06-20 11:32:18 UTC (rev 65504)
@@ -0,0 +1,320 @@
+#include "local_proto.h"
+
+void variogram_type(int code, char *type)
+{
+ switch (code) {
+ case 0:
+ strcpy(type, "horizontal");
+ break;
+ case 1:
+ strcpy(type, "vertical");
+ break;
+ case 2:
+ strcpy(type, "bivariate");
+ break;
+ case 3:
+ strcpy(type, "anisotropic");
+ break;
+ }
+}
+
+void write2file_basics(struct int_par *xD, struct opts *opt)
+{
+ struct write *report = &xD->report;
+
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Input ***\n\n");
+ fprintf(report->fp, "- vector layer: %s\n", opt->input->answer);
+ switch (xD->v3) {
+ case TRUE:
+ fprintf(report->fp, "- is 3D: yes\n");
+ break;
+ case FALSE:
+ fprintf(report->fp, "- is 3D: no\n");
+ break;
+ }
+ // todo: move to end
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "*** Output *** \n\n");
+ switch (xD->i3) {
+ case FALSE:
+ fprintf(report->fp, "- raster layer: %s\n", opt->output->answer);
+ fprintf(report->fp, "- is 3D: no\n");
+ break;
+ case TRUE:
+ fprintf(report->fp, "- volume layer: %s\n", opt->output->answer);
+ fprintf(report->fp, "- is 3D: yes\n");
+ break;
+ }
+}
+
+void write2file_vector( struct int_par *xD, struct points *pnts)
+{
+ struct select in_reg = pnts->in_reg;
+ struct write *report = &xD->report;
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Input vector layer properties ***\n\n");
+ fprintf(report->fp, "# of points (total): %d\n", in_reg.total);
+ fprintf(report->fp, "# of points (in region): %d\n", in_reg.n);
+ fprintf(report->fp, "# of points (out of region): %d\n\n", in_reg.out);
+ fprintf(report->fp, "- extent:\n");
+ switch (xD->i3) {
+ case TRUE:
+ fprintf(report->fp, "xmin=%f ymin=%f zmin=%f\n", pnts->r_min[0], pnts->r_min[1], pnts->r_min[2]);
+ fprintf(report->fp, "xmax=%f ymax=%f zmax=%f\n", pnts->r_max[0], pnts->r_max[1], pnts->r_max[2]);
+ break;
+ case FALSE:
+ fprintf(report->fp, "xmin=%f ymin=%f\n", pnts->r_min[0], pnts->r_min[1]);
+ fprintf(report->fp, "xmax=%f ymax=%f\n", pnts->r_max[0], pnts->r_max[1]);
+ break;
+ }
+}
+
+void write2file_values(struct write *report, const char *column)
+{
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Values to be interpolated ***\n\n");
+ fprintf(report->fp, "- attribute column: %s\n", column);
+}
+
+void write2file_varSetsIntro(int code, struct write *report)
+{
+ char type[12];
+ variogram_type(code, type);
+
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Variogram settings - %s ***\n\n", type);
+}
+
+void write2file_varSets(struct write *report, struct parameters *var_pars)
+{
+ char type[12];
+ variogram_type(var_pars->type, type);
+
+ fprintf(report->fp, "- number of lags in %s direction: %d\n", type, var_pars->nLag);
+ fprintf(report->fp, "- lag distance (%s): %f\n", type, var_pars->lag);
+
+ if (var_pars->type == 2) {
+ fprintf(report->fp, "- number of lags in vertical direction: %d\n", var_pars->nLag_vert);
+ fprintf(report->fp, "- lag distance (vertical): %f\n", var_pars->lag);
+ }
+
+ fprintf(report->fp, "- azimuth: %f°\n", RAD2DEG(var_pars->dir));
+ fprintf(report->fp, "- angular tolerance: %f°\n", RAD2DEG(var_pars->td));
+
+ switch (var_pars->type) {
+ case 2:
+ fprintf(report->fp, "- maximum distance in horizontal direction (1/3): %f\n", var_pars->max_dist);
+ break;
+ default:
+ fprintf(report->fp, "- maximum distance (1/3): %f\n", var_pars->max_dist);
+ break;
+ }
+}
+
+void write2file_variogram_E(struct int_par *xD, struct parameters *var_pars, mat_struct *c_M)
+{
+ int i, j;
+ int n_h;
+ double *h, *vert;
+ double *c;
+ double *gamma;
+ struct write *report = &xD->report;
+
+ char type[12];
+ variogram_type(var_pars->type, type);
+
+ c = &c_M->vals[0];
+ gamma = &var_pars->gamma->vals[0];
+
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Experimental variogram - %s ***\n\n", type);
+
+ if (var_pars->type == 2) { // bivariate variogram:
+ vert = &var_pars->vertical.h[0];
+
+ for (i=0; i<var_pars->vertical.nLag; i++) { // write header - verts
+ if (i == 0) {
+ fprintf(report->fp, " lagV ||"); // column for h
+ }
+ fprintf(report->fp, " %f ||", *vert);
+ vert++;
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i < var_pars->vertical.nLag; i++) { // write header - h c gamma
+ if (i == 0) {
+ fprintf(report->fp, " lagHZ ||"); // column for h
+ }
+ fprintf(report->fp, " c | ave ||");
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i<var_pars->vertical.nLag; i++) { // write header - h c gamma
+ if (i == 0) {
+ fprintf(report->fp, "--------||"); // column for h
+ }
+ fprintf(report->fp, "-----------");
+ }
+ fprintf(report->fp, "\n");
+ }
+ else { // univariate variogram
+ fprintf(report->fp, " lag || # of pairs | average\n");
+ fprintf(report->fp, "------------------------------------\n");
+ }
+
+ // Write values
+ h = var_pars->type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
+ n_h = var_pars->type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
+ for (i=0; i < n_h; i++) {
+ fprintf(report->fp, "%f ||", *h);
+ if (var_pars->type == 2) { // bivariate variogram
+ for (j=0; j<var_pars->vertical.nLag; j++) {
+ fprintf(report->fp, " %d | %f ||", (int) *c, *gamma);
+ c++;
+ gamma++;
+ } // end for j
+ fprintf(report->fp, "\n");
+ for (j=0; j<var_pars->vertical.nLag; j++) {
+ fprintf(report->fp, "-----------");
+ }
+ fprintf(report->fp, "\n");
+ }
+ else { // univariate variogram
+ fprintf(report->fp, " %d | %f\n", (int) *c, *gamma);
+ c++;
+ gamma++;
+ }
+ h++;
+ }
+ fprintf(report->fp, "------------------------------------\n");
+}
+
+void write2file_variogram_T(struct write *report)
+{
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Theoretical variogram ***\n\n");
+}
+
+void write_temporary2file(struct int_par *xD, struct parameters *var_pars)
+{
+ // local variables
+ int type = var_pars->type;
+ int nLag = type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
+ int nLag_vert = type == 2 ? var_pars->vertical.nLag : var_pars->nLag;
+ double *h = type == 2 ? var_pars->horizontal.h: var_pars->h;
+ double *vert = type == 2 ? var_pars->vertical.h: var_pars->h;;
+ double *gamma = var_pars->gamma->vals;
+ double sill = var_pars->sill;
+
+ int i, j; // index
+ FILE *fp;
+ int file_length;
+
+ /* Introduction */
+ switch (type) {
+ case 0: // horizontal variogram
+ fp = fopen("variogram_hz_tmp.txt", "w");
+ 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\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 (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\n", var_pars->type); // write type
+ break;
+
+ 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 (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\n", var_pars->type); // write type
+ fprintf(fp,"%d\n", var_pars->vertical.nLag); // write # of lags
+ fprintf(fp,"%f\n", var_pars->vertical.lag); // write size of lag
+ fprintf(fp,"%f\n", var_pars->vertical.max_dist); // write maximum distance
+
+ fprintf(fp,"%d\n", var_pars->horizontal.nLag); // write # of lags
+ fprintf(fp,"%f\n", var_pars->horizontal.lag); // write size of lag
+ fprintf(fp,"%f\n", var_pars->horizontal.max_dist); // write maximum distance
+ break;
+
+ 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 (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\n", var_pars->type); // write type
+ fprintf(fp,"%f\n", xD->aniso_ratio); // write ratio of anisotropy
+ break;
+ }
+
+ if (type != 2) {
+ fprintf(fp,"%d\n", nLag); // write # of lags
+ fprintf(fp,"%f\n", var_pars->lag); // write size of lag
+ fprintf(fp,"%f\n", var_pars->max_dist); // write maximum distance
+ }
+ if (type != 1) {
+ fprintf(fp,"%f\n", var_pars->td); // write maximum distance
+ }
+
+ switch (type) {
+ case 2: // bivariate variogram
+ nLag = var_pars->horizontal.nLag;
+ nLag_vert = var_pars->vertical.nLag;
+ // write h
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *h);
+ h++;
+ }
+ // write vert
+ for (i=0; i<nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *vert);
+ vert++;
+ }
+ // write gamma
+ for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *gamma);
+ gamma++;
+ }
+ fprintf(fp,"%f\n", var_pars->horizontal.sill);// write sill
+ fprintf(fp,"%f\n", var_pars->vertical.sill);// write sill
+ break;
+ default:
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f %f\n", *h, *gamma);
+ h++;
+ gamma++;
+ }
+ fprintf(fp,"%f", sill);// write sill
+ break;
+ }
+
+ fclose(fp);
+}
More information about the grass-commit
mailing list