[GRASS-SVN] r62319 - in grass-addons/grass7/vector: . v.kriging
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Oct 21 06:21:36 PDT 2014
Author: evas
Date: 2014-10-21 06:21:36 -0700 (Tue, 21 Oct 2014)
New Revision: 62319
Added:
grass-addons/grass7/vector/v.kriging/
grass-addons/grass7/vector/v.kriging/Makefile
grass-addons/grass7/vector/v.kriging/geostat.cpp
grass-addons/grass7/vector/v.kriging/getval.cpp
grass-addons/grass7/vector/v.kriging/local_proto.h
grass-addons/grass7/vector/v.kriging/main.cpp
grass-addons/grass7/vector/v.kriging/quantile.cpp
grass-addons/grass7/vector/v.kriging/stat.cpp
grass-addons/grass7/vector/v.kriging/utils.cpp
grass-addons/grass7/vector/v.kriging/utils_kriging.cpp
grass-addons/grass7/vector/v.kriging/utils_raster.cpp
grass-addons/grass7/vector/v.kriging/utils_write.cpp
grass-addons/grass7/vector/v.kriging/v.kriging.html
Log:
new module v.kriging
Added: grass-addons/grass7/vector/v.kriging/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.kriging/Makefile 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,19 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.kriging
+
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GMATHLIB) $(GISLIB) $(IOSTREAMLIB) -lpcl_apps -lpcl_common -lpcl_features -lpcl_filters -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_visualization
+DEPENDENCIES= $(RASTER3DDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GMATHLIBDEP) $(GISDEP) $(IOSTREAMDEP)
+EXTRA_INC = $(VECT_INC) -I/usr/include/pcl-1.7 -I/usr/include/eigen3 $(PCL_INCLUDE_DIRS)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+EXTRA_CFLAGS=-fopenmp
+EXTRA_LIBS=$(GISLIB) -lgomp $(MATHLIB) -lboost_system
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+default: cmd
+endif
Added: grass-addons/grass7/vector/v.kriging/geostat.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/geostat.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/geostat.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,720 @@
+#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 pcl_utils *pclp, struct var_par *pars)
+{
+ // Variogram type
+ struct parameters *var_par, *var_par_vert;
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts;
+ switch (type) {
+ case 0:
+ var_par = &pars->hz;
+ pcl_pnts = pclp->pnts_hz;
+ break;
+ case 1:
+ var_par = &pars->vert;
+ pcl_pnts = pclp->pnts_vert;
+ break;
+ default:
+ var_par = &pars->fin;
+ pcl_pnts = pclp->pnts;
+ break;
+ }
+
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *r; // xyz coordinates of input points
+ double *vals = pnts->invals; // values to be used for interpolation
+ int phase = xD->phase;
+ int function = var_par->function; // type of theoretical variogram
+
+ int nLag = var_par->nLag; // # of horizontal bins
+ int nLag_vert;
+ double *vert;
+ double dir = var_par->dir; // azimuth
+ double td = var_par->td; // angle tolerance
+ double lag = var_par->lag; // size of horizontal bin
+ double lag_vert;
+
+ /*switch (type) {
+ case 1:
+ td = var_par->tz;
+ break;
+ case 2:
+ td = var_par->horizontal->td;
+ tz = var_par->vertical->td;
+ break;
+ default:
+ td = var_par->td;
+ break;
+ }*/
+
+ // depend on variogram type:
+ int k; // index of variogram to be computed
+ int nrows = nLag;
+
+ double max_dist = var_par->max_dist; // max distance of interpolated points
+ double radius = SQUARE(max_dist); // radius of interpolation in the horizontal plane
+ struct write *report = &xD->report;
+
+ if (type == 2) {
+ nLag_vert = var_par->nLag_vert; // # of horizontal bins
+ lag_vert = var_par->lag_vert; // size of horizontal bin
+ }
+
+ int ncols = type == 2 ? nLag_vert : 1;
+
+ // 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
+
+ double *i_vals, *j_vals; // values located on the point couples
+
+ int n1 = n+1; // number of points + 1
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+
+ pcl::PointXYZ *searchPt; // search point (each input point)
+ std::vector<int> ind; // vector of indices
+ std::vector<float> sqDist;
+ int n_ind; // number of elements (including 1st irrelevant)
+ int *indices; // save indices into matrix to use them later
+ int *idx; // index to be set (output)
+ int *ii, *i0; // difference of indices between rellevant input points
+ int *ip; // index to be set into the index matrix (input)
+
+ double gamma_sum; // 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_temp; // gamma matrix (hz, vert or bivar)
+ mat_struct *c_temp; // matrix of # of dissimilarities
+ double *gamma; // pointer to gamma matrix
+ double *c; // pointer to c matrix
+
+ unsigned int percents = 50;
+ int counter;
+ int count;
+
+ double gamma_stat; // sum of gamma elements (non-nan)
+ int gamma_n; // # of gamma elements (non-nan)
+ double gamma_sill; // sill
+
+ double nugget, sum_nugget = 0.;
+ int n_nugget = 0;
+
+ /* 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
+ var_par->h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
+
+ if (type == 2) {
+ var_par->vert = (double *) G_malloc(nLag_vert * sizeof(double)); // vector of bins
+ }
+
+ // control initialization:
+ if (dr == NULL || var_par->h == NULL || (type == 2 && var_par->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..."));
+ }
+
+ r = &pnts->r[0]; // set up pointer to input coordinates
+ indices = (int *) G_malloc(n1*n * sizeof(int)); // matrix of indices of neighbouring points
+ if (indices == 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..."));
+ }
+
+ kd_tree->setInputCloud(pcl_pnts);
+ searchPt = &pcl_pnts->points[0];
+
+ // initialize:
+ c_temp = G_matrix_init(nrows, ncols, nrows); // temporal matrix of counts
+ gamma_temp = G_matrix_init(nrows, ncols, nrows); // temporal matrix (vector) of gammas
+
+ if (c_temp == NULL || gamma_temp == 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_temp->vals[0];
+ gamma = &gamma_temp->vals[0];
+
+ // set up starting values
+ gamma_stat = 0.; //
+ gamma_n = 0; //
+
+ if (percents) {
+ G_percent_reset();
+ }
+ counter = nrows * ncols;
+ count = 0;
+
+ if (type == 2) {
+ vert = &var_par->vert[0];
+ }
+
+ for (b = 0; b < ncols; b++) {
+ if (type == 2) {
+ *vert = (b+1) * lag_vert;
+ }
+
+ // set up vector of distances...
+ h = &var_par->h[0]; // ... horizontal
+
+ for (s = 0; s < nrows; s++) { // process of the horizontal piece (isotrophy!!!)
+ *h = (s+1) * lag; // distance of horizontal lag
+
+ // for every bs cycle ...
+ i_vals = &vals[0]; // ... i-vals start at the beginning
+ gamma_sum = 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 every input point...
+ if (b == 0 && s == 0) { // ... create vector of points in neighborhood
+ if (kd_tree->radiusSearch(*searchPt, max_dist, ind, sqDist) > 0) {
+ n_ind = ind.size (); // size of indices to 1st row
+ ip = &ind.data ()[1]; // values of indices to the next rows
+ idx = &indices[n1*i]; // begin new column
+ *idx = n_ind-1; // size of index vector as 1st element of row
+ for (j = 1; j < n_ind; j++) {
+ idx++; // next row
+ *idx = *ip; // set indices to the rest of the rows
+ ip++; // next index
+ } // end setting elements
+ } // end radiusSearch
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("There are no points in the surrounding of point no. %d..."), i);
+ } // end error
+ searchPt++; // next search point
+ } // end set up indices
+
+ idx = &indices[n1*i]; // ... new column of matrix is attached
+ ii = &indices[n1*i+1];
+ i0 = &indices[n1*i+1];
+ n_ind = *idx; // number of neighbouring points
+
+ j_vals = &vals[*ii];
+ for (j = 1; j < n_ind; j++) { // points within searchRadius
+ if (*ii > i) {
+ coord_diff(i, *ii, r, dr); // coordinate differences
+
+ /* Variogram processing */
+ if (type == 1) { // in case of vertical variogram...
+ tv = 0.5 * PI - atan2(dr[2], sqrt(SQUARE(dr[0]) + SQUARE(dr[1]))); // ... compute zenith angle instead of bearing
+ td = 0.5 * PI; // 0.25 usually
+ }
+ else
+ tv = atan2(*(dr+1), *dr); // bearing
+
+ ddir = tv - dir; // difference between bearing and azimuth
+ if (fabsf(ddir) <= td) { // is the difference acceptable according to tolerance?
+ // test distance
+ rv = type == 1 ? 0. : SQUARE(dr[0]) + SQUARE(dr[1]); // horizontal distance
+ if (type == 1 || type == 3) { // anisotropic only
+ rv += SQUARE(dr[2]); // vertical or anisotropic distance
+ }
+ rvh = sqrt(rv) - *h; // difference between distance and lag
+ if (rv <= radius && 0. <= rvh && rvh <= lag) { // distance test
+ if (type == 2) { // vertical test for bivariate variogram
+ rvh = *(dr+2) - *vert;
+ if (fabsf(rvh) <= lag_vert) { // elevation test
+ goto delta_V;
+ }
+ else {
+ goto end;
+ }
+ }
+ delta_V:
+ dv = *j_vals - *i_vals; // difference of values located on pair of points i, j
+ gamma_sum += dv * dv; // sum of squared differences
+ cpls += 1.; // count of elements
+ } // end distance test
+ } // end bearing test
+ //} // end dr2 test
+ } // end point selection
+ end:
+ *i0 = *ii;
+ ii++;
+ j_vals += *ii - *i0;
+ } // end j
+ i_vals++;
+ } // end i
+
+ if (isnan(gamma_sum) || cpls == 0.0) {
+ err0++;
+ goto gamma_isnan;
+ }
+
+ gamma_E = 0.5 * gamma_sum / cpls; // element of gamma matrix
+ *c = cpls; // # of values that have been used to compute gamma_e
+ *gamma = gamma_E;
+
+ gamma_stat += gamma_E; // sum of gamma elements (non-nan)
+ gamma_n++; // # of gamma elements (non-nan)
+
+ gamma_isnan:
+ h++;
+ c++;
+ gamma++;
+ count++;
+ } // end s
+
+ if (type == 2) {
+ vert++;
+ }
+ } // end b
+
+ free(indices);
+
+ 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_par->gamma = G_matrix_copy(gamma_temp);
+ var_par->c = G_matrix_copy(c_temp);
+ plot_experimental_variogram(xD, var_par->h, gamma_temp, var_par);
+
+ if (report->name) {
+ write2file_variogram_E(xD, var_par);
+ }
+
+ // Compute sill
+ if (phase < 2) {
+ switch (type) {
+ case 2:
+ int nh, nv;
+ double c_h, c_v, gamma_h, gamma_v;
+ char type[12];
+ nh = nv = 0;
+ gamma_h = gamma_v = 0.;
+
+ for (i=0; i<nrows; i++) {
+ c_h = G_matrix_get_element(c_temp, i, 0);
+ if (c_h != 0.) {
+ gamma_h += G_matrix_get_element(gamma_temp, i, 0);
+ nh++;
+ }
+ }
+
+ for (j=0; j<ncols; j++) {
+ c_v = G_matrix_get_element(c_temp, 0, j);
+ if (c_v != 0.) {
+ gamma_v += G_matrix_get_element(gamma_temp, 0, j);
+ nv++;
+ }
+ }
+
+ var_par->horizontal.sill = gamma_h / nh;
+ var_par->vertical.sill = gamma_v / nv;
+ break;
+
+ default:
+ var_par->sill = gamma_stat / gamma_n; // mean of gamma elements
+
+ variogram_type(var_par->type, type);
+ G_message(_("Default sill of %s variogram: %f"), type, var_par->sill);
+ break;
+ }
+ } // end compute sill
+
+ write_temporary2file(xD, var_par, gamma_temp);
+
+ G_matrix_free(c_temp);
+ G_matrix_free(gamma_temp);
+}
+
+/* theoretical variogram */
+void T_variogram(int type, struct opts opt, struct parameters *var_par, struct write *report)
+{
+ // report
+ if (report->name) {
+ time(&report->now);
+ if (type != 1) {
+ fprintf(report->fp, "\nComputation of theoretical variogram started on %s\n", ctime(&report->now));
+ }
+ }
+
+ // set up:
+ var_par->type = type;
+
+ char *variogram;
+ switch (type) {
+ case 0: // horizontal
+ var_par->nugget = atof(opt.nugget_hz->answer);
+ var_par->h_range = atof(opt.range_hz->answer);
+ if (opt.sill_hz->answer) {
+ var_par->sill = atof(opt.sill_hz->answer);
+ }
+ variogram = opt.function_var_hz->answer;
+ if (report->name) {
+ fprintf(report->fp, "Parameters of horizontal variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ case 1: // vertical
+ var_par->nugget = atof(opt.nugget_vert->answer);
+ var_par->h_range = atof(opt.range_vert->answer);
+ if (opt.sill_vert->answer) {
+ var_par->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_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ case 2: // bivariate variogram
+ if (opt.function_var_hz->answer == NULL && opt.function_var_vert->answer == NULL) { // linear function
+ int nZ, nL, nr, i, j, nc;
+ double *h, *vert, *gamma, *c;
+ mat_struct *gR;
+
+ var_par->function = 5;
+
+ 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
+ c = &var_par->c->vals[0];
+
+ // 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 (*c != 0.) { // todo: upgrade: !isnan(*c) L434 too
+ nr++;
+ }
+ gamma++;
+ c++;
+ } // 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
+ c = &var_par->c->vals[0];
+
+ if (var_par->function == 5) { // in case of bivariate variogram
+ vert = &var_par->vert[0]; // use separate variable for vertical direction
+ }
+
+ for (i = 0; i < nZ; i++) {
+ h = &var_par->h[0];
+ for (j = 0; j < nL; j++) {
+ if (*c != 0.) { // 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++;
+ c++;
+ } // 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
+ G_debug(0,"%f %f %f", var_par->T->vals[0], var_par->T->vals[1], var_par->T->vals[2]);
+ 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.
+ else {
+ var_par->const_val = 0;
+ }
+
+ // 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 + c\n", var_par->T->vals[0], var_par->T->vals[1], var_par->T->vals[2]);
+ }
+ }
+
+ else {
+ var_par->horizontal.nugget = atof(opt.nugget_hz->answer);
+ var_par->horizontal.h_range = atof(opt.range_hz->answer);
+ //var_par->horizontal.variogram = opt.variogram_hz->answer;
+ if (opt.sill_hz->answer) {
+ var_par->horizontal.sill = atof(opt.sill_hz->answer);
+ }
+
+ var_par->vertical.nugget = atof(opt.nugget_vert->answer);
+ var_par->vertical.h_range = atof(opt.range_vert->answer);
+ //var_par->vertical.variogram = opt.variogram_vert->answer;
+ if (opt.sill_hz->answer) {
+ var_par->vertical.sill = atof(opt.sill_vert->answer);
+ }
+
+ var_par->horizontal.function = set_function(opt.function_var_hz->answer, var_par, report);
+ var_par->vertical.function = set_function(opt.function_var_vert->answer, var_par, report);
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Nugget effect (hz): %f\n", var_par->horizontal.nugget);
+ fprintf(report->fp, "Sill (hz): %f\n", var_par->horizontal.sill);
+ fprintf(report->fp, "Range (hz): %f\n", var_par->horizontal.h_range);
+ fprintf(report->fp, "Function: %s\n\n", opt.function_var_hz->answer);
+ fprintf(report->fp, "Nugget effect (vert): %f\n", var_par->vertical.nugget);
+ fprintf(report->fp, "Sill (vert): %f\n", var_par->vertical.sill);
+ fprintf(report->fp, "Range (vert): %f\n", var_par->vertical.h_range);
+ fprintf(report->fp, "Function: %s\n", opt.function_var_vert->answer);
+ }
+ }
+
+ plot_var(TRUE, var_par->h, var_par->gamma, var_par); // Plot variogram using gnuplot
+ break;
+ case 3: // anisotropic
+ var_par->nugget = atof(opt.nugget_final->answer);
+ var_par->h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_par->sill = atof(opt.sill_final->answer);
+ }
+ variogram = opt.function_var_final->answer;
+ if (report->name) {
+ fprintf(report->fp, "Parameters of anisotropic variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_par->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_par->sill);
+ fprintf(report->fp, "Range: %f\n", var_par->h_range);
+ }
+ break;
+ }
+
+ if (type != 2) {
+ var_par->function = set_function(variogram, var_par, report);
+ plot_var(FALSE, var_par->h, var_par->gamma, var_par); // Plot variogram using gnuplot
+ }
+}
+
+void ordinary_kriging(struct int_par *xD, struct reg_par *reg, struct points *pnts, struct pcl_utils *pclp, struct var_par *pars, struct output *out)
+{
+ // Local variables
+ int n = pnts->n;
+ double *r = pnts->r; // xyz coordinates of input poinst
+ double *vals = pnts->invals; // values to be used for interpolation
+ struct write *report = &xD->report;
+ struct write *crossvalid = &xD->crossvalid;
+ struct parameters *var_par;
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts;
+ double ratio;
+
+ var_par = &pars->fin;
+ pcl_pnts = pclp->pnts;
+ 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;
+ double *r0; // xyz coordinates of cell/voxel centre
+ mat_struct *GM;
+ mat_struct *GM_sub; // submatrix of the points that are rellevant to the interpolation due the distance
+ mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
+ mat_struct *g0; // vector of diffences between known and unknown values estimated using distances and the theoretical variogram
+ 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));
+ pcl::KdTreeFLANN<pcl::PointXYZ> kd_tree;
+ kd_tree.setInputCloud(pcl_pnts); // Set up kd-tree
+ pcl::PointXYZ cellCent;
+ std::vector<int> ind;
+ std::vector<float> sqDist;
+
+ double radius;
+ FILE *fp;
+
+ if (report->name) {
+ 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) {
+ goto constant_voxel_centre;
+ }
+
+ radius = var_par->max_dist;
+
+ GM = set_up_G(pnts, var_par, &xD->report);
+ var_par->GM = G_matrix_copy(GM);
+ if (n < 1000) {
+ GM_Inv = G_matrix_inverse(GM);
+ }
+
+ if (crossvalid->name) {
+ crossvalidation(xD, pnts, pcl_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) {
+ goto constant_voxel_val;
+ }
+
+ cell_centre(col, row, dep, xD, reg, r0, var_par); // Output point
+ if (n < 1000)
+ goto small_sample;
+
+ cellCent.x = *r0;
+ cellCent.y = *(r0+1);
+ cellCent.z = *(r0+2); // if 2D: 0.
+
+ if (kd_tree.radiusSearch(cellCent, radius, ind, sqDist) > 0 ) {
+ GM_sub = submatrix(ind, GM, report);
+ GM_Inv = G_matrix_inverse(GM_sub);
+ G_matrix_free(GM_sub);
+
+ small_sample:
+ g0 = set_up_g0(xD, ind, pnts, 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
+
+ if (n >= 1000) {
+ G_matrix_free(GM_Inv);
+ G_matrix_free(g0);
+ }
+
+ rslt_OK = result(ind, pnts, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ //if (pnts->trend != NULL) {
+ /* normal gravity: rslt_OK += pnts->trend->vals[0] * r0[0] + pnts->trend->vals[1] * r0[0]*r0[1]*r0[2]/ratio + pnts->trend->vals[2] * r0[0]*r0[1] + pnts->trend->vals[3] * r0[0]*r0[2]/ratio + pnts->trend->vals[4] * r0[1]*r0[2]/ratio + pnts->trend->vals[5] * r0[2]/ratio + pnts->trend->vals[6];*/
+ //}
+
+ // Create output
+ constant_voxel_val:
+ if (xD->const_val == 1)
+ rslt_OK = (double) *vals;
+ if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
+ if (report->name) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error writing result into output layer..."));
+ }
+ } // end if searchRadius
+ 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..."));
+ }
+ } // 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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/getval.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/getval.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,551 @@
+#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);
+
+ /* 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;
+ //if (i<78)
+ // G_debug(0, "%d %f", i, *vals);
+ }
+ 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++) {
+ /* normal gravity:
+ G_matrix_set_element(A, i, 0, *x);
+ G_matrix_set_element(A, i, 1, *x * *y * *z);
+ G_matrix_set_element(A, i, 2, *x * *y);
+ G_matrix_set_element(A, i, 3, *x * *z);
+ G_matrix_set_element(A, i, 4, *z * *y);
+ G_matrix_set_element(A, i, 5, *z);
+ G_matrix_set_element(A, i, 6, 1.);
+ */
+ 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);
+ }
+
+ return values;
+}
+
+/* get coordinates of input points */
+void read_points(struct Map_info *map, struct reg_par *reg, struct points *point, struct pcl_utils *pclp, struct int_par xD, const char *zcol, int field, struct write *report)
+{
+ int type; // check elements of vector layer
+ 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)
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+
+ double *rx, *ry, *rz, *r; // pointers to coordinates
+ pcl::PointXYZ *pclr; // pointer to PCL structure
+ 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
+ r = (double *) G_malloc(n * 3 * sizeof(double)); // x0 y0 z0 ... xn yn zn
+ indices = (int *) G_malloc(n * sizeof(int));
+ point->r_min = (double *) G_malloc(3 * sizeof(double)); // minimum
+ point->r_max = (double *) G_malloc(3 * sizeof(double)); // maximum
+
+ rx = &r[0]; // set up pointer to x coordinate
+ ry = &r[1]; // set up pointer to y coordinate
+ rz = &r[2]; // set up pointer to z coordinate
+ index = &indices[0]; // set up pointer to indices
+
+ pcl_pnts->points.resize (n); // set up size of PCL structure
+ pclr = &pcl_pnts->points[0]; // set up pointer to PCL structure
+
+ // 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
+
+ pclr->x = Points->x[0]; // set up x coordinate to PCL structure
+ pclr->y = Points->y[0]; // set up y coordinate to PCL structure
+
+ if (xD.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 (n_in_reg < 79)
+ // G_debug(0,"%d %f %f %f", n_in_reg, *rx, *ry, *rz);
+ pclr->z = Points->z[0]; // set up z coordinate to PCL structure
+ 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
+ pclr->z = *z_attr; // set up x coordinate to PCL structure
+ } // end else (2,5D layer)
+ } // end if rb, rt
+ else
+ goto out_of_region;
+ } // end if 3D interpolation
+
+ else { // 2D interpolation
+ *rz = 0.; // set up z coordinate
+ pclr->z = 0.; // set up z coordinate to PCL structure
+ }
+
+ // Find extends
+ if (n_in_reg == 0)
+ point->r_min = point->r_max = triple(*rx, *ry, *rz);
+
+ else {
+ point->r_min = triple(MIN(*point->r_min, *rx), MIN(*(point->r_min+1), *ry), MIN(*(point->r_min+2), *rz));
+ point->r_max = triple(MAX(*point->r_max, *rx), MAX(*(point->r_max+1), *ry), MAX(*(point->r_max+2), *rz));
+ } // 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
+ pclr++; // go to new PCL point
+ goto finish;
+ } // end in region
+
+ out_of_region:
+ out_reg++;
+ continue;
+
+ finish:
+ if (zcol != NULL) // todo: skontroluj, co robi pre 3d
+ 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->r = (double *) G_malloc(3 * point->n * sizeof(double));
+ memcpy(point->r, r, 3 * n_in_reg * sizeof(double));
+
+ point->in_reg.indices = (int *) G_malloc(n_in_reg * sizeof(int));
+ memcpy(point->in_reg.indices, indices, n_in_reg * sizeof(int));
+
+ pcl_pnts->points.resize (n_in_reg);
+
+ pclp->pnts = pcl_pnts;
+ kd_tree->setInputCloud(pclp->pnts);
+ pclp->kd_tree = kd_tree;
+
+ 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"));
+ }
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts_hz (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl_pnts_hz->points.resize (n);
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts_vert (new pcl::PointCloud<pcl::PointXYZ>);
+ pcl_pnts_vert->points.resize (n);
+
+ for (int i=0; i<n_in_reg; i++) {
+ pcl_pnts_hz->points[i].x = pcl_pnts->points[i].x;
+ pcl_pnts_hz->points[i].y = pcl_pnts->points[i].y;
+ pcl_pnts_hz->points[i].z = 0.;
+
+ pcl_pnts_vert->points[i].x = 0.;
+ pcl_pnts_vert->points[i].y = 0.;
+ pcl_pnts_vert->points[i].z = pcl_pnts->points[i].z;
+ } // end switch var_par->function
+
+ pclp->pnts_hz = pcl_pnts_hz;
+ pclp->pnts_vert = pcl_pnts_vert;
+
+ //write2file_vector(&xD, point);
+}
+
+extern "C" {
+ 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 */
+ }
+}
+
+void read_tmp_vals(const char *file_name, struct parameters *var_par, struct int_par *xD)
+{
+ FILE *fp;
+
+ fp = fopen(file_name, "r");
+ if (fp == NULL)
+ G_fatal_error(_("File is missing..."));
+
+ else {
+ 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)
+ 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) {
+ xD->report.name = (char *) G_malloc(file_length * sizeof(char));
+ if (fscanf(fp, "%s", xD->report.name) == 0)
+ 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
+
+ no_file:
+
+ switch (type) {
+ case 2:
+ int j, nLag_vert;
+ double lag_vert, max_dist_vert;
+ double *v_elm, *g_elm, *c_elm;
+ double sill_hz, sill_vert;
+
+ 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->h = (double *) G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->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->vert = (double *) G_malloc(nLag_vert * sizeof(double));
+ v_elm = &var_par->vert[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->c = G_matrix_init(nLag, nLag_vert, nLag);
+ c_elm = &var_par->c->vals[0];
+ for (i=0; i<nLag_vert; i++) {
+ for (j=0; j<nLag; j++) {
+ if (fscanf(fp, "%lf", c_elm) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ c_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->nLag_vert = nLag_vert;
+ var_par->lag_vert = lag_vert;
+ var_par->max_dist_vert = max_dist_vert;
+ var_par->horizontal.sill = sill_hz;
+ var_par->vertical.sill = sill_vert;
+ break;
+ default:
+ if (type == 3) {
+ if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0)
+ 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->c = G_matrix_init(nLag, 1, nLag);
+ c = &var_par->c->vals[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 %lf", h_elm, c, gamma) < 3)
+ G_fatal_error(_("Nothing to scan..."));
+ h_elm++;
+ c++;
+ gamma++;
+ }
+
+ if (fscanf(fp, "%lf", &sill) == 0)
+ G_fatal_error(_("Nothing to scan..."));
+ var_par->sill = sill;
+ break;
+ }
+
+ 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/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.kriging/local_proto.h (rev 0)
+++ grass-addons/grass7/vector/v.kriging/local_proto.h 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,256 @@
+#ifndef __LOCAL_PROTO__
+#define __LOCAL_PROTO__
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+
+#include <pcl/point_types.h>
+#include <pcl/io/pcd_io.h>
+#include <pcl/point_cloud.h>
+#include <pcl/kdtree/impl/kdtree_flann.hpp>
+
+extern "C" {
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/la.h>
+#include <grass/raster3d.h>
+#include <grass/glocale.h>
+}
+
+#ifndef PI
+#define PI M_PI
+#endif
+
+#ifndef SQUARE
+#define SQUARE(a) (a*a)
+#define POW3(a) (a*a*a)
+#define POW4(a) (a*a*a*a)
+#endif
+
+struct opts
+{
+ struct Option *input, *output, *phase, *report, *crossvalid, *function_var_hz, *function_var_vert, *function_var_final, *form_file, *field, *intpl, *zcol, *var_dir_hz, *var_dir_vert, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert, *nugget_final, *sill_hz, *sill_vert, *sill_final, *range_hz, *range_vert, *range_final;
+};
+
+struct flgs
+{
+ struct Flag *d23; /* 2D/3D interpolation */
+ struct Flag *bivariate;
+ struct Flag *univariate;
+ struct Flag *detrend;
+};
+
+struct select
+{
+ int n; // # of selected
+ int out; // # of the others
+ int total; // total number
+ int *indices; // indices of selected
+};
+
+struct points // inputs
+{
+ int n; // number of points
+ double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
+ double *r_min; // min coords
+ double *r_max; // max coords
+ double center[3]; // center
+ double *invals; // values to be interpolated
+ struct select in_reg; // points in region
+ mat_struct *trend;
+};
+
+struct pcl_utils
+{
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pnts;
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree;
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pnts_hz;
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree_A;
+
+ pcl::PointCloud<pcl::PointXYZ>::Ptr pnts_vert;
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree_xy;
+};
+
+struct bivar
+{
+ int vert;
+ char *variogram;
+ int function;
+ double sill;
+ double nugget;
+ double h_range;
+};
+
+struct parameters
+{
+ int function;
+ int type;
+ int const_val;
+ double dir;
+ double td; // range of directions
+
+ double radius; // radius
+ double max_dist; // maximum distance
+ double max_dist_vert;
+
+ int nLag; // number of length pieces
+ double lag; // maximum distance between nearest neighbours (variogram lag)
+
+ int nLag_vert;
+ double lag_vert;
+
+ double *h; // horizontal intervals used to estimate experimental variogram
+ double *vert; // vertical intervals used to estimate experimental variogram
+ mat_struct *c; // number of elements in each lag (final variogram only)
+ mat_struct *gamma; // value of experimental variogram
+ double nugget;
+ double sill;
+ double part_sill;
+ double h_range;
+
+ struct bivar horizontal;
+ struct bivar vertical;
+
+ mat_struct *A;
+ mat_struct *T; // coefficients of theoretical variogram
+ mat_struct *GM; // matrix of diffences between point values based on the distances and the theoretical variogram
+
+ char *name; // name of input vector layer
+ char term[12]; // output format - gnuplot terminal
+ char ext[4]; // output format - file format
+};
+
+struct var_par // parameters of experimental variogram
+{
+ struct parameters hz;
+ struct parameters vert;
+ struct parameters fin;
+
+ //double Lmax; // maximal length in horizontal direction
+ //double dzmax; // maximal elevation difference
+ //double radius; // radius for kd-tree
+};
+
+struct write
+{
+ int write2file; // write to file or not
+ char *name; // filename
+ FILE *fp;
+ time_t now;
+};
+
+struct int_par // Interpolation settings
+{
+ int i3; // TRUE = 3D interpolation, FALSE = 2D interpolation (user sets by flag)
+ char v3; // TRUE = 3D layer, FALSE = 2D layer (according to Vect_is_3d())
+ int phase;
+ int bivar;
+ int univar;
+ double aniso_ratio;
+ int const_val;
+ struct write report;
+ struct write crossvalid;
+};
+
+extern "C" {
+ struct reg_par // Region settings -> output extent and resolution
+ {
+ struct Cell_head reg_2d; // region for 2D interpolation
+ RASTER3D_Region reg_3d; // region for 3D interpolation
+ double west; // region.west
+ double east;
+ double north; // region.north
+ double south;
+ double bot; // region.bottom
+ double top;
+ double ew_res; // east-west resolution
+ double ns_res; // north-south resolution
+ double bt_res; // bottom-top resolution
+ int nrows; // number of rows
+ int ncols; // number of cols
+ int ndeps; // number of deps
+ int nrcd;
+ };
+
+ struct output
+ {
+ DCELL *dcell;
+ int fd_2d;
+ RASTER3D_Map *fd_3d;
+ char *name; // name of output 2D/3D raster
+ };
+}
+
+double *get_col_values(struct Map_info *, struct int_par *, struct points *, int, const char *, int);
+void test_normality(int , double *, struct write *);
+void read_points(struct Map_info *, struct reg_par *, struct points *, struct pcl_utils *, struct int_par, const char *, int, struct write *);
+double min(double *, struct points *);
+double max(double *, struct points *);
+
+void coord_diff(int, int, double *, double *);
+double distance_diff(double *);
+double radius_hz_diff(double *);
+double *triple(double, double, double);
+double lag_distance(int, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *, struct write *);
+int lag_number(double, double *);
+void variogram_restricts(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *);
+void geometric_anisotropy(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr);
+double find_intersect_x(double *, double *, double *, double *, struct write *);
+double find_intersect_y(double *, double *, double *, double *, double , struct write *);
+mat_struct *LSM(mat_struct *, mat_struct *);
+mat_struct *nonlin_LMS(int , double *, double *);
+void E_variogram(int, struct int_par *, struct points *, struct pcl_utils *, struct var_par *);
+void T_variogram(int, struct opts, struct parameters *, struct write *);
+void ordinary_kriging(struct int_par *, struct reg_par *, struct points *, struct pcl_utils *, struct var_par *, struct output *);
+
+int set_function(char *, struct parameters *, struct write *);
+double RBF(double *);
+double linear(double, double, double);
+double exponential(double, double, double, double);
+double spherical(double, double, double, double);
+double gaussian(double, double, double, double);
+double variogram_fction(struct parameters *, double *);
+void set_gnuplot(char *, struct parameters *);
+void plot_experimental_variogram(struct int_par *, double *, mat_struct *, struct parameters *);
+void plot_var(int, double *, mat_struct *, struct parameters *);
+void variogram_type(int, char *);
+void write2file_basics(struct int_par *, struct opts *);
+void write2file_vector(struct int_par *, struct points *);
+void write2file_values(struct write *, const char *);
+void write2file_varSetsIntro(int, struct write *);
+void write2file_varSets(struct write *, struct parameters *);
+void write2file_variogram_E(struct int_par *, struct parameters *);
+void write2file_variogram_T(struct write *);
+void write_temporary2file(struct int_par *, struct parameters *, mat_struct *);
+void read_tmp_vals(const char *, struct parameters *, struct int_par *);
+
+mat_struct *set_up_G(struct points *, struct parameters *, struct write *);
+mat_struct *set_up_g0(struct int_par *, std::vector<int>, struct points *, double *, struct parameters *);
+mat_struct *submatrix(std::vector<int> , mat_struct *, struct write *);
+double result(std::vector<int>, struct points *, mat_struct *);
+
+void crossvalidation(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *);
+void cell_centre(unsigned int, unsigned int, unsigned int, struct int_par *, struct reg_par *, double *, struct parameters *);
+
+extern "C" {
+void get_region_pars(struct int_par *, struct reg_par *);
+void open_layer(struct int_par *, struct reg_par *, struct output *);
+int write2layer(struct int_par *, struct reg_par *, struct output *, unsigned int, unsigned int, unsigned int, double);
+}
+
+static inline double get_quantile(int);
+static void get_slot_counts(int, double *);
+static void initialize_bins(void);
+static void fill_bins(int, double *);
+static int compare(const void *, const void *);
+static void sort_bins(void);
+static void compute_quantiles(int, double, struct write *);
+double quantile(double, int, double *, struct write *);
+#endif
Added: grass-addons/grass7/vector/v.kriging/main.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/main.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/main.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,485 @@
+
+/****************************************************************
+ *
+ * MODULE: v.kriging
+ * AUTHOR: Eva Stopková
+ * in case of functions taken from another modules,
+ * they are cited above the function
+ * or at the beginning of the file (e.g. quantile.cpp
+ * that uses slightly modified functions
+ * taken from the module r.quantile (Clemens, G.))
+ * PURPOSE: Module interpolates the values to two- or three-dimensional grid using input values
+ * located on 2D/3D point vector layer. Interpolation method
+ * Ordinary kriging has been extended for 3D points (v = f(x,y) -> v = f(x,y,z)).
+ *
+ * COPYRIGHT: (C) 2012-2014 Eva Stopková and by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ **************************************************************/
+#include "local_proto.h"
+
+
+int main(int argc, char *argv[])
+{
+ // Vector layer and module
+ struct Map_info map; // Input vector map
+ struct GModule *module; // Module
+
+ struct reg_par reg; // Region parameters
+ struct points pnts; // Input points (coordinates, extent, values, etc.)
+ struct pcl_utils pclp; // PCL structure to store coordinates of input points
+
+ // Geostatistical parameters
+ struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
+ struct var_par var_par; // Variogram (experimental and theoretical)
+
+ // Outputs
+ struct output out; // Output layer properties
+ FILE *fp;
+
+ // Settings
+ int field;
+ struct opts opt;
+ struct flgs flg;
+
+ /* ------- Module creation ------- */
+ module = G_define_module();
+ G_add_keyword(_("Raster"));
+ G_add_keyword(_("3D raster"));
+ G_add_keyword(_("Ordinary kriging - for 2D and 3D data"));
+ module->description =
+ _("Interpolates 2D or 3D raster based on input values located on 2D or 3D point vector layer (method ordinary kriging extended to 3D).");
+
+ // Setting options
+ opt.input = G_define_standard_option(G_OPT_V_INPUT); // Vector input layer
+ opt.input->label = _("Name of input vector points map");
+
+ flg.d23 = G_define_flag(); // Option to process 2D or 3D interpolation
+ flg.d23->key = '2';
+ flg.d23->description = _("Force 2D interpolation even if input is 3D");
+ flg.d23->guisection = _("3D");
+
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
+
+ opt.phase = G_define_option();
+ opt.phase->key = "phase";
+ opt.phase->options = "initial, middle, final";
+ opt.phase->description = _("Phase of interpolation");
+ opt.phase->required = YES;
+
+ opt.output = G_define_option(); // Output layer
+ opt.output->key = "output";
+ opt.output->description =
+ _("Name for output 2D/3D raster map");
+ opt.output->guisection = _("Final");
+
+ opt.report = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.report->key = "report";
+ opt.report->description = _("Report file");
+ opt.report->required = NO;
+ opt.report->guisection = _("Files");
+
+ opt.crossvalid = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.crossvalid->key = "crossvalid";
+ opt.crossvalid->description = _("File with crooss validation results");
+ opt.crossvalid->required = NO;
+ opt.crossvalid->guisection = _("Files");
+
+ flg.bivariate = G_define_flag();
+ flg.bivariate->key = 'b';
+ flg.bivariate->description = _("Compute bivariate variogram");
+ flg.bivariate->guisection = _("Middle");
+
+ flg.univariate = G_define_flag();
+ flg.univariate->key = 'u';
+ flg.univariate->description = _("Compute univariate variogram");
+ flg.univariate->guisection = _("Middle");
+
+ opt.function_var_hz = G_define_option(); // Variogram type
+ opt.function_var_hz->key = "hz_function";
+ opt.function_var_hz->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_hz->description = _("Type of horizontal variogram function");
+ opt.function_var_hz->guisection = _("Middle");
+
+ opt.function_var_vert = G_define_option(); // Variogram type
+ opt.function_var_vert->key = "vert_function";
+ opt.function_var_vert->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_vert->description = _("Type of vertical variogram function");
+ opt.function_var_vert->guisection = _("Middle");
+
+ opt.function_var_final = G_define_option(); // Variogram type
+ opt.function_var_final->key = "final_function";
+ opt.function_var_final->options = "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_final->description = _("Type of final variogram function");
+ opt.function_var_final->guisection = _("Final");
+
+ flg.detrend = G_define_flag();
+ flg.detrend->key = 't';
+ flg.detrend->description = _("Eliminate trend if variogram is parabolic...");
+ flg.detrend->guisection = _("Initial");
+
+ opt.form_file = G_define_option(); // Variogram plot - various output formats
+ opt.form_file->key = "fileformat";
+ opt.form_file->options = "cdr,dxf,eps,tex,pdf,png,svg";
+ opt.form_file->description = _("Various file formats");
+ opt.form_file->guisection = _("Middle");
+
+ opt.intpl = G_define_standard_option(G_OPT_DB_COLUMN); // Input values for interpolation
+ opt.intpl->key = "icolumn";
+ opt.intpl->description =
+ _("Column containing input values for interpolation");
+ opt.intpl->required = YES;
+
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN); // Column with z coord (2D points)
+ opt.zcol->key = "zcolumn";
+ opt.zcol->description =
+ _("Column containing z coordinates - set only if you want to process 3D interpolation based on 2D point layer");
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("3D");
+
+ opt.var_dir_hz = G_define_option();
+ opt.var_dir_hz->key = "azimuth";
+ opt.var_dir_hz->type = TYPE_DOUBLE;
+ opt.var_dir_hz->required = NO;
+ opt.var_dir_hz->answer = "0.0";
+ opt.var_dir_hz->description =
+ _("Azimuth of computing variogram (isotrophic)");
+ opt.var_dir_hz->guisection = _("Initial");
+
+ opt.var_dir_vert = G_define_option();
+ opt.var_dir_vert->key = "zenith_angle";
+ opt.var_dir_vert->type = TYPE_DOUBLE;
+ opt.var_dir_vert->required = NO;
+ opt.var_dir_vert->answer = "0.0";
+ opt.var_dir_vert->description =
+ _("Zenith angle of computing variogram (isotrophic)");
+ opt.var_dir_vert->guisection = _("Initial");
+
+ opt.nL = G_define_option();
+ opt.nL->key = "lpieces";
+ opt.nL->type = TYPE_INTEGER;
+ opt.nL->required = NO;
+ opt.nL->description = _("Count of length pieces");
+ opt.nL->guisection = _("Initial");
+
+ opt.nZ = G_define_option();
+ opt.nZ->key = "vpieces";
+ opt.nZ->type = TYPE_INTEGER;
+ opt.nZ->required = NO;
+ opt.nZ->description =
+ _("Count of vertical pieces (set only for counting 3D variogram)");
+ opt.nZ->guisection = _("Initial");
+
+ opt.td_hz = G_define_option();
+ opt.td_hz->key = "td";
+ opt.td_hz->type = TYPE_DOUBLE;
+ opt.td_hz->answer = "90.0";
+ opt.td_hz->description = _("Angle of variogram processing");
+ opt.td_hz->required = NO;
+ opt.td_hz->guisection = _("Initial");
+
+ opt.nugget_hz = G_define_option();
+ opt.nugget_hz->key = "hz_nugget";
+ opt.nugget_hz->type = TYPE_DOUBLE;
+ opt.nugget_hz->description = _("Nugget effect");
+ opt.nugget_hz->required = NO;
+ opt.nugget_hz->guisection = _("Middle");
+
+ opt.nugget_vert = G_define_option();
+ opt.nugget_vert->key = "vert_nugget";
+ opt.nugget_vert->type = TYPE_DOUBLE;
+ opt.nugget_vert->description = _("Nugget effect");
+ opt.nugget_vert->required = NO;
+ opt.nugget_vert->guisection = _("Middle");
+
+ opt.nugget_final = G_define_option();
+ opt.nugget_final->key = "final_nugget";
+ opt.nugget_final->type = TYPE_DOUBLE;
+ opt.nugget_final->description = _("Nugget effect");
+ opt.nugget_final->required = NO;
+ opt.nugget_final->guisection = _("Final");
+
+ opt.sill_hz = G_define_option();
+ opt.sill_hz->key = "hz_sill";
+ opt.sill_hz->type = TYPE_DOUBLE;
+ opt.sill_hz->description = _("Variogram sill");
+ opt.sill_hz->required = NO;
+ opt.sill_hz->guisection = _("Middle");
+
+ opt.sill_vert = G_define_option();
+ opt.sill_vert->key = "vert_sill";
+ opt.sill_vert->type = TYPE_DOUBLE;
+ opt.sill_vert->description = _("Variogram sill");
+ opt.sill_vert->required = NO;
+ opt.sill_vert->guisection = _("Middle");
+
+ opt.sill_final = G_define_option();
+ opt.sill_final->key = "final_sill";
+ opt.sill_final->type = TYPE_DOUBLE;
+ opt.sill_final->description = _("Variogram sill");
+ opt.sill_final->required = NO;
+ opt.sill_final->guisection = _("Final");
+
+ opt.range_hz = G_define_option();
+ opt.range_hz->key = "hz_range";
+ opt.range_hz->type = TYPE_DOUBLE;
+ opt.range_hz->description = _("Variogram range");
+ opt.range_hz->required = NO;
+ opt.range_hz->guisection = _("Middle");
+
+ opt.range_vert = G_define_option();
+ opt.range_vert->key = "vert_range";
+ opt.range_vert->type = TYPE_DOUBLE;
+ opt.range_vert->description = _("Variogram range");
+ opt.range_vert->required = NO;
+ opt.range_vert->guisection = _("Middle");
+
+ opt.range_final = G_define_option();
+ opt.range_final->key = "final_range";
+ opt.range_final->type = TYPE_DOUBLE;
+ opt.range_final->description = _("Variogram range");
+ opt.range_final->required = NO;
+ opt.range_final->guisection = _("Final");
+ /* --------------------------------------------------------- */
+
+ G_gisinit(argv[0]);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* Get parameters from the parser */
+ if (strcmp(opt.phase->answer, "initial") == 0)
+ xD.phase = 0; // estimate range
+ else if (strcmp(opt.phase->answer, "middle") == 0)
+ xD.phase = 1; // estimate anisotropic variogram
+ else if (strcmp(opt.phase->answer, "final") == 0)
+ xD.phase = 2; // compute kriging
+
+ // Open report file if desired
+ if (opt.report->answer) {
+ xD.report.write2file = TRUE;
+ xD.report.name = opt.report->answer;
+ xD.report.fp = fopen(xD.report.name, "w");
+ time(&xD.report.now);
+ fprintf(xD.report.fp, "v.kriging started on %s\n\n", ctime(&xD.report.now));
+ G_message(_("Report is being written to %s..."), xD.report.name);
+ }
+ else
+ xD.report.write2file = FALSE;
+
+ if (opt.crossvalid->answer) {
+ xD.crossvalid.write2file = TRUE;
+ xD.crossvalid.name = opt.crossvalid->answer;
+ xD.crossvalid.fp = fopen(xD.crossvalid.name, "w");
+ }
+ else
+ xD.crossvalid.write2file = FALSE;
+
+ var_par.hz.name = var_par.vert.name = var_par.fin.name = opt.input->answer; // set name of variogram
+ var_par.hz.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
+ var_par.vert.dir = opt.var_dir_vert->answer ? DEG2RAD(atof(opt.var_dir_vert->answer)) : 0.; // Zenith angle
+ var_par.fin.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
+
+ var_par.hz.td = DEG2RAD(atof(opt.td_hz->answer)); // Angle of variogram processing
+ //var_par.vert.td = DEG2RAD(atof(opt.td_vert->answer)); // Angle of variogram processing
+
+ /*if (phase > 0) {
+ //var_par.sill = atof(opt.sill->answer);
+ var_par.nugget = atof(opt.nugget->answer);
+ var_par.h_range = atof(opt.range->answer);
+ if (/*(var_par.sill != -1. && var_par.sill <= 0.) || *//*(var_par.nugget != -1. && var_par.nugget < 0.) || (var_par.h_range != -1. && var_par.h_range <= 0.)) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("Variogram parameters should be positive..."));
+ } // error
+ } */
+
+ /*if (opt.nL->answer) { // Test if nL have been set up (optional)
+ if (var_par.nL < 1) // Invalid value
+ G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
+ else
+ var_par.nL = atof(opt.nL->answer); // todo osetrit aby sa neprepisovalo
+ }
+
+ if (opt.variogram->answer) // Function of theoretical variogram
+ set_function(opt.variogram->answer, &var_par, &xD.report);*/
+
+ if (opt.form_file->answer) { // Plotting variogram
+ set_gnuplot(opt.form_file->answer, &var_par.hz);
+ set_gnuplot(opt.form_file->answer, &var_par.vert);
+ set_gnuplot(opt.form_file->answer, &var_par.fin);
+ }
+ else {
+ strcpy(var_par.hz.term, ""); /* TO DO: ked pouzivatel zada hlupost */
+ strcpy(var_par.vert.term, "");
+ strcpy(var_par.fin.term, "");
+ }
+
+ xD.bivar = flg.bivariate->answer == TRUE ? TRUE : FALSE;
+ xD.univar = flg.univariate->answer == TRUE ? TRUE : FALSE;
+ if (xD.bivar == TRUE && xD.univar == TRUE) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("You should mark either univariate, or bivariate variogram, not both of them..."));
+ } // error
+
+ /* ---------------------------------------------------------- */
+ Vect_set_open_level(2); // Open input vector map
+
+ if (0 > Vect_open_old2(&map, opt.input->answer, "", opt.field->answer)) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("Unable to open vector map <%s>"), opt.input->answer);
+ } // error
+ Vect_set_error_handler_io(&map, NULL);
+ /* ---------------------------------------------------------- */
+
+ /* Perform 2D or 3D interpolation ? */
+ xD.i3 = flg.d23->answer ? FALSE : TRUE; // 2D/3D interpolation
+ xD.v3 = Vect_is_3d(&map); // 2D/3D layer
+
+ /* What could happen:
+ *-------------------
+ * 3D interpolation + 3D points = 3D GIS
+ * 3D interpolation + 2D points = 2,5D -> 3D GIS (needs attribute column with z and 3D region)
+ * 2D interpolation + 3D points = 3D -> 2,5D GIS
+ * 2D interpolation + 2D points = 2,5D GIS */
+
+ // 3D interpolation
+ if (xD.i3 == TRUE) {
+ if (xD.v3 == FALSE) { // 2D input
+ if (!opt.zcol->answer) {
+ 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 based on 2D input, please set attribute column containing z coordinates or switch to 2D interpolation."));
+ }
+ } // end if zcol == NULL
+ // 3D or 2,5D input
+ if (opt.nZ->answer) { // Test if nZ have been set up (optional)
+ if (var_par.vert.nLag < 1) // Invalid value
+ G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
+ else
+ var_par.vert.nLag = atof(opt.nZ->answer);
+ } // end if nZ->answer
+ /* TO DO: chyba, ak nie je cele cislo */
+ } // end if 3D interpolation
+
+ /* 2D interpolation */
+ else
+ var_par.vert.nLag = -1; // abs will be used in next steps
+
+ field = Vect_get_field_number(&map, opt.field->answer);
+ if (xD.report.write2file == TRUE)
+ write2file_basics(&xD, &opt);
+ /* ---------------------------------------------------------- */
+
+ /* Get... */
+ get_region_pars(&xD, ®); // ... region parameters
+ read_points(&map, ®, &pnts, &pclp, xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
+ pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
+
+ /* Estimate 2D/3D variogram */
+ switch (xD.phase) {
+ case 0:
+ /* Determine maximal horizontal (and vertical) distance + lags */
+ var_par.hz.type = 0;
+ var_par.vert.type = 1;
+
+ variogram_restricts(&xD, &pnts, pclp.pnts_hz, &var_par.hz);
+ variogram_restricts(&xD, &pnts, pclp.pnts_vert, &var_par.vert);
+
+ if (var_par.vert.nLag > var_par.hz.nLag) {
+ var_par.vert.nLag = var_par.hz.nLag;
+ var_par.vert.lag = var_par.vert.max_dist / var_par.vert.nLag;
+ write2file_varSets(&xD.report, &var_par.vert);
+ }
+
+ E_variogram(0, &xD, &pnts, &pclp, &var_par);
+ E_variogram(1, &xD, &pnts, &pclp, &var_par);
+ G_message(_("You may continue to computing theoretical variograms (middle phase)..."));
+ goto end;
+ case 1:
+ read_tmp_vals("variogram_hz_tmp.txt", &var_par.hz, &xD);
+ read_tmp_vals("variogram_vert_tmp.txt", &var_par.vert, &xD);
+
+ if (xD.report.name) {
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ }
+
+ T_variogram(0, opt, &var_par.hz, &xD.report); // Theoretical variogram
+ T_variogram(1, opt, &var_par.vert, &xD.report); // Theoretical variogram
+
+
+ /* difference between ranges
+ - if greater than 5% - bivariate
+ - if smaller than 5% - anisotropic
+ */
+ double sill95, diff_sill;
+ sill95 = var_par.hz.sill > var_par.vert.sill ? 0.05 * var_par.hz.sill : 0.05 * var_par.vert.sill;
+ diff_sill = fabsf(var_par.hz.sill - var_par.vert.sill);
+ if (xD.bivar == TRUE || (!flg.univariate->answer && diff_sill > sill95)) { // zonal anisotropy
+ var_par.fin.type = 2;
+ var_par.fin.td = var_par.hz.td;
+ var_par.fin.max_dist = var_par.hz.max_dist;
+ var_par.fin.max_dist_vert = var_par.vert.max_dist;
+ variogram_restricts(&xD, &pnts, pclp.pnts, &var_par.fin);
+ E_variogram(2, &xD, &pnts, &pclp, &var_par);
+ }
+
+ else if (xD.univar == TRUE || (!flg.bivariate->answer && diff_sill <= sill95)) { // geometric anisotropy
+ var_par.fin.type = 3;
+ var_par.fin.max_dist = var_par.hz.max_dist;
+ var_par.fin.td = var_par.hz.td;
+
+ xD.aniso_ratio = var_par.hz.h_range / var_par.vert.h_range;
+ geometric_anisotropy(&xD, &pnts, pclp.pnts);
+
+ variogram_restricts(&xD, &pnts, pclp.pnts, &var_par.fin);
+
+ E_variogram(3, &xD, &pnts, &pclp, &var_par);
+ }
+ G_message(_("You may continue to interpolate values (final phase)..."));
+ goto end;
+ case 2:
+ if (!opt.output->answer)
+ G_fatal_error(_("Please set up name of output layer..."));
+ out.name = opt.output->answer; // Output layer name
+
+ read_tmp_vals("variogram_final_tmp.txt", &var_par.fin, &xD);
+
+ if (xD.report.name) {
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ if (xD.report.fp == NULL)
+ G_fatal_error(_("Cannot open the file..."));
+ }
+
+ if (var_par.fin.type == 3)
+ geometric_anisotropy(&xD, &pnts, pclp.pnts);
+
+ T_variogram(var_par.fin.type, opt, &var_par.fin, &xD.report); // Theoretical variogram
+ break;
+ }
+
+ /* Ordinary kriging (including 2D/3D raster creation) */
+ ordinary_kriging(&xD, ®, &pnts, &pclp, &var_par, &out);
+ /* ---------------------------------------------------------- */
+
+ end:
+ Vect_close(&map); // Close vector map
+
+ exit(EXIT_SUCCESS);
+}
+
Added: grass-addons/grass7/vector/v.kriging/quantile.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/quantile.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/quantile.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -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 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 (int 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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/stat.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/stat.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,790 @@
+#include "local_proto.h"
+
+#define GNUPLOT "gnuplot -persist"
+
+// make coordinate triples xyz
+double *triple(double x, double y, double z)
+{
+ double *t;
+ t = (double *) G_malloc(3 * sizeof(double));
+
+ *t = x;
+ *(t+1) = y;
+ *(t+2) = z;
+
+ return t;
+}
+
+// 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 size of the lag
+double lag_distance(int direction, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par, struct write *report)
+{
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz of input points
+ int fction = var_par->function;
+ double max_dist = (var_par->type == 2 && direction == 3) ? var_par->max_dist_vert : var_par->max_dist; // maximum horizontal distance
+
+ int i, j, k;
+ pcl::PointXYZ *searchPt;
+ searchPt = &pcl_pnts->points[0];
+
+ pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree (new pcl::KdTreeFLANN<pcl::PointXYZ>);
+ kd_tree->setInputCloud(pcl_pnts);
+
+ std::vector<int> ind; // indices of NN
+ std::vector<float> sqDist; // sqDistances of NN
+
+ double dr[3]; // coordinate differences
+ double quant_dist2 = var_par->radius; // quantile
+ double lagNN; // squared quantile - lag
+ double *dist2, *d; // radius
+ dist2 = (double *) G_malloc(n * sizeof(double));
+ d = &dist2[0];
+
+ int perc5, zeros;
+ int nN;
+ perc5 = (int) round(0.05 * n);
+ double dist0;
+ int j0 = 0;
+
+ for (i=0; i < n; i++) {
+ zeros = 0; // # of zero distances
+ dist0 = 0.; // initial value of 5% distance
+ while (dist0 != 0.) { // while 5% distance is zero
+ perc5 += zeros; // 5% plus points in zero distance
+ // find 5% nearest points
+ if (kd_tree->nearestKSearch(*searchPt, perc5, ind, sqDist) > 0) {
+ j = ind.size () - 1; // store index of 5%-th point
+ dist0 = sqDist.data ()[j]; // 5% distance
+
+ zeros = 0; // count zero distances
+ for (k=j0; k <= j; k++) {
+ if (sqDist.data ()[k] == 0.)
+ zeros++;
+ }
+
+ if (zeros > 0) {
+ dist0 = 0.;
+ j0 = j+1;
+ }
+ } // end nearestKSearch
+
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp, "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error searching nearest neighbours of point %d..."), i);
+ } // end error
+ } // end while dist0 == 0
+
+ coord_diff(i, j, r, dr); // coordinate differences
+
+ // squared distance
+ switch (direction) {
+ case 12: // horizontal
+ *d = SQUARE(dr[0]) + SQUARE(dr[1]);
+ break;
+ case 3: // vertical
+ *d = SQUARE(dr[2]);
+ break;
+ case 0: // all
+ *d = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
+ break;
+ }
+ if (*d != 0.)
+ quant_dist2 = MIN(quant_dist2, *d);
+ d++;
+ searchPt++;
+ }
+
+ lagNN = sqrt(quant_dist2);
+
+ 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, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par)
+{
+ 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)
+
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ // 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_par->type) {
+ case 1:
+ var_par->max_dist = dr[2]; // zmax - zmin
+ var_par->radius = SQUARE(var_par->max_dist); // anisotropic distance (todo: try also dz)
+ break;
+ default:
+ var_par->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
+ var_par->max_dist = sqrt(var_par->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
+ break;
+ }
+
+ if (report->name)
+ write2file_varSetsIntro(var_par->type, report);
+
+ int dimension;
+ switch (var_par->type) {
+ case 0:
+ dimension = 12;
+ break;
+ case 1:
+ dimension = 3;
+ break;
+ case 3:
+ dimension = 0;
+ break;
+ }
+
+ if (var_par->type == 2) {
+ var_par->lag = lag_distance(12, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag = lag_number(var_par->lag, &var_par->max_dist);
+ var_par->max_dist = var_par->nLag * var_par->lag;
+
+ var_par->lag_vert = lag_distance(3, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag_vert = lag_number(var_par->lag_vert, &var_par->max_dist_vert);
+ var_par->max_dist_vert = var_par->nLag_vert * var_par->lag_vert;
+ }
+
+ else {
+ var_par->lag = lag_distance(dimension, pnts, pcl_pnts, var_par, report); // the shortest distance between NN in horizontal direction
+ var_par->nLag = lag_number(var_par->lag, &var_par->max_dist);
+ }
+
+ write2file_varSets(report, var_par);
+}
+
+void geometric_anisotropy(struct int_par *xD, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts)
+{
+ int i;
+ double ratio = xD->aniso_ratio;
+ double *z;
+ z = &pnts->r[2];
+
+ for (i=0; i<pnts->n; i++) {
+ *z = pcl_pnts->points[i].z = ratio * *z;
+ z += 3;
+ }
+}
+
+// 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;
+}
+
+mat_struct *nonlin_LMS(int n, double *dist, double *gamma)
+{
+ int i=0, j;
+ double ctrl=1.;
+ double a, b, c, mean_dy = 1., sum_dy, *h, *g, *add;
+ mat_struct *y, *param, *param0, *J, *yf, *JT, *JTJ, *JT_Inv, *JTdy, *dy, *delta;
+ y = G_matrix_init(n, 1, n);
+ param = G_matrix_init(3, 1, 3);
+ param0 = G_matrix_init(3, 1, 3);
+ J = G_matrix_init(n, 3, n);
+ yf = G_matrix_init(n, 1, n);
+
+ // matrix of input elements
+ g = &gamma[0];
+ for (j=0; j<n; j++) {
+ G_matrix_set_element(y, j, 0, *g);
+ g++;
+ }
+
+ while (fabsf(mean_dy) >= 0.0001) {
+ if (i==0) { // set initial parameters
+ G_matrix_set_element(param, 0, 0, 0.);
+ G_matrix_set_element(param, 1, 0, 45.);
+ G_matrix_set_element(param, 2, 0, 45.);
+ }
+ a = param->vals[0];
+ b = param->vals[1];
+ c = param->vals[2];
+ G_debug(0,"%d a=%f b=%f c=%f", i, a, b, c);
+
+ h = &dist[0];
+ for (j=0; j<n; j++) {
+ // Jacobi matrix
+ // todo: rozdelit na variogramy
+ G_matrix_set_element(J, j, 0, 1.);
+ G_matrix_set_element(J, j, 1, 1. - exp(-3. * *h / c));
+ G_matrix_set_element(J, j, 2, -3. * b * *h/SQUARE(c) * exp(-3. * *h*c));
+
+ // f(a,b,c)
+ G_matrix_set_element(yf, j, 0, b * (1. - exp(-3. * *h/c)));
+ G_debug(0, "%f %f %f", G_matrix_get_element(J, j, 0), G_matrix_get_element(J, j, 1), G_matrix_get_element(J, j, 2));
+ h++;
+ }
+
+ // differences
+ dy = G_matrix_subtract(y, yf);
+ sum_dy = 0.;
+ add = &dy->vals[0];
+ for (j=0; j<n; j++) {
+ sum_dy += SQUARE(*add);
+ add++;
+ }
+ mean_dy = sum_dy / n;
+ G_debug(0,"%d mean=%f", i, mean_dy);
+
+ JT = G_matrix_transpose(J);
+ JTJ = G_matrix_product(JT, J);
+ JT_Inv = G_matrix_inverse(JTJ);
+ JTdy = G_matrix_product(JT, dy);
+ delta = G_matrix_product(JT_Inv, JTdy);
+ G_debug(0,"a=%f b=%f c=%f", delta->vals[0], delta->vals[1], delta->vals[2]);
+
+ param0 = G_matrix_copy(param);
+ param = G_matrix_subtract(param0, delta);
+ i++;
+ if (i>3)
+ G_fatal_error(_("stoj"));
+ }
+ G_debug(0,"nugget=%f sill=%f range=%f", param->vals[0], param->vals[1], param->vals[2]);
+ return param;
+}
+
+// set up type of function according to GUI (user)
+int set_function(char *variogram, struct parameters *var_par, struct write *report)
+{
+ int function;
+ if (strcmp(variogram, "linear") == 0) {
+ function = 0;
+ // var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "parabolic") == 0) {
+ function = 1;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "exponential") == 0) {
+ function = 2;
+ // var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "spherical") == 0) {
+ function = 3;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "gaussian") == 0) {
+ function = 4;
+ //var_par->bivar = FALSE;
+ }
+ else if (strcmp(variogram, "bivariate") == 0) {
+ function = 5;
+ //var_par->bivar = TRUE;
+ }
+ 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_par)
+{
+ if (strcmp(fileformat,"cdr") == 0) {
+ strcpy(var_par->term, "corel");
+ strcpy(var_par->ext, "cdr");
+ }
+ if (strcmp(fileformat,"dxf") == 0) {
+ strcpy(var_par->term, "dxf");
+ strcpy(var_par->ext, "dxf");
+ }
+ if (strcmp(fileformat,"eps") == 0) {
+ strcpy(var_par->term, "postscript");
+ strcpy(var_par->ext, "eps");
+ }
+ if (strcmp(fileformat,"pdf") == 0) {
+ strcpy(var_par->term, "pdfcairo");
+ strcpy(var_par->ext, "pdf");
+ }
+ if (strcmp(fileformat,"png") == 0) {
+ strcpy(var_par->term, "png");
+ strcpy(var_par->ext, "png");
+ }
+ if (strcmp(fileformat,"svg") == 0) {
+ strcpy(var_par->term, "svg");
+ strcpy(var_par->ext, "svg");
+ }
+}
+
+// plot experimental variogram
+void plot_experimental_variogram(struct int_par *xD, double *h_exp, mat_struct *gamma_exp, 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;
+ double *gamma_var; // pointer to gamma matrix
+ double *c; // pointer to number of elements
+ FILE *gp; // pointer to file
+
+ nr = gamma_exp->rows; // # of rows of gamma matrix
+ nc = gamma_exp->cols; // # of cols of gamma matrix
+
+ gamma_var = &gamma_exp->vals[0]; // values of gamma matrix
+ c = &var_par->c->vals[0];
+
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (gp == NULL)
+ G_fatal_error(_("Error writing file"));
+
+ /* 3D experimental variogram */
+ if (bivar == TRUE) {
+ vert = &var_par->vert[0];
+ for (i=0; i < nc+1; i++) { // for each row (nZ)
+ h = &var_par->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 (*c == 0.)
+ fprintf(gp, "NaN"); // write other elements
+ else
+ fprintf(gp, "%f", *gamma_var); // write other elements
+
+ if (j != nr)
+ fprintf(gp, " "); // spaces between elements
+ gamma_var++;
+ c++;
+ }
+ } // 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 = &h_exp[0];
+ for (i=0; i < nr; ++i) {
+ if (*c == 0.)
+ fprintf(gp, "%f NaN\n", *h);
+ else
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ c++;
+ h++;
+ gamma_var++;
+ }
+ }
+ 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);
+
+ //remove("dataE.dat");
+}
+
+// plot experimental and theoretical variogram
+void plot_var(int bivar, double *h_exp, mat_struct *gamma_exp, struct parameters *var_par)
+{
+ 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;
+
+ switch (bivar) {
+ case FALSE:
+ function = var_par->function;
+ nugget = var_par->nugget;
+ sill = var_par->sill - nugget;
+ h_range = var_par->h_range;
+ break;
+ case TRUE:
+ if (var_par->function != 5) {
+ func_h = var_par->horizontal.function;
+ func_v = var_par->vertical.function;
+
+ nugget_h = var_par->horizontal.nugget;
+ sill_h = var_par->horizontal.sill - nugget_h;
+ h_range_h = var_par->horizontal.h_range;
+
+ nugget_v = var_par->vertical.nugget;
+ sill_v = var_par->vertical.sill - nugget_v;
+ h_range_v = var_par->vertical.h_range;
+ }
+ break;
+ }
+
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert; // pointer to vertical bins
+ //double *T_cf; // coefficients of theoretical variogram - vals
+ double h_ratio;
+ double g_teor;
+ double *c;
+ double *gamma_var; // pointer to gamma matrix
+ double *T; // coefficients of theoretical variogram - matrix
+ FILE *gp; // pointer to file
+
+ nr = gamma_exp->rows; // # of rows of gamma matrix
+ nc = gamma_exp->cols; // # of cols of gamma matrix
+
+ c = &var_par->c->vals[0];
+ gamma_var = &gamma_exp->vals[0]; // values of gamma matrix
+
+ if (var_par->type != 2) { // univariate variogram
+ gp = fopen("dataE.dat","w"); // open file to write experimental variogram
+ if (gp == NULL)
+ G_fatal_error(_("Error writing file"));
+
+ h = &h_exp[0];
+ for (i=0; i < nr; i++) {
+ if (*c == 0.)
+ fprintf(gp, "%f NaN\n", *h); // write other elements
+ else
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ h++;
+ c++;
+ gamma_var++;
+ }
+ }
+ else {
+ gp = fopen("dataE.dat","r"); // open file to write experimental variogram
+ if (gp == NULL)
+ 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 (gp == NULL)
+ G_fatal_error(_("Error opening file..."));
+
+ double dd;
+ double hh;
+ double hv[2];
+ h = &h_exp[0];
+
+
+
+ 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_par->T->vals[0]; // values
+ vert = &var_par->vert[0];
+ for (i=0; i < nc+1; i++) { // for each row
+ h = &var_par->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_par, 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_par->term,"") != 0) {
+ //fprintf(gp, "set terminal %s size 750,450\n", var_par->term);
+ fprintf(gp, "set terminal wxt size 750,450\n", var_par->type);
+ switch (var_par->type) {
+ case 0:
+ if (bivar == TRUE)
+ fprintf(gp, "set output \"variogram_bivar.%s\"\n", var_par->ext);
+ else
+ fprintf(gp, "set output \"variogram_horizontal.%s\"\n", var_par->ext);
+ break;
+ case 1:
+ fprintf(gp, "set output \"variogram_vertical.%s\"\n", var_par->ext);
+ break;
+ case 2:
+ if (bivar == TRUE)
+ fprintf(gp, "set output \"variogram_bivariate.%s\"\n", var_par->ext);
+ else
+ fprintf(gp, "set output \"variogram_anisotropic.%s\"\n", var_par->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_par->type);
+ }
+
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp, "set title \"Experimental and theoretical 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\", 'dataT.dat' every ::1:1 matrix with lines title \"theoretical variogram\" palette\n");
+ }
+
+ else { // univariate variogram*/
+ //if (xD->i3 == TRUE) { // 3D
+ if (var_par->type == 0) { // horizontal variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D hz) of <%s>\"\n", var_par->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_par->type == 1) { // vertical variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D vert) of <%s>\"\n", var_par->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_par->type == 3) // anisotropic variogram
+ fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_par->name);
+ //}
+
+/*else // 2D
+ fprintf(gp, "set title \"Experimental and theoretical variogram (2D) of <%s>\"\n", var_par->name);*/
+
+ //if (var_par->type == 2/* || xD->i3 == FALSE*/) {
+ switch (var_par->function) {
+ case 0: // linear
+ 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]);
+ break;
+ case 1: // parabolic
+ fprintf(gp, "set label \"parabolic: gamma(h) = %f*h^2 + %f\" at screen 0.25,0.90\n", var_par->T->vals[0], var_par->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_par->nugget, var_par->sill - var_par->nugget, var_par->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_par->nugget, var_par->sill - var_par->nugget, var_par->h_range, var_par->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_par->nugget, var_par->sill - var_par->nugget, var_par->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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_kriging.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_kriging.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,473 @@
+#include "local_proto.h"
+
+ __inline double square(double x)
+ {
+ return x*x;
+ }
+
+// formulation of variogram functions
+double variogram_fction(struct parameters *var_par, double *dr)
+{
+ // Local variables
+ int variogram = var_par->function; // theoretical variogram
+ int type = var_par->type;
+ double nugget;
+ double part_sill;
+ double h_range;
+ double *T;
+
+ if (type == 2) {
+ 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 (bi- or univariate)
+
+ for (int 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 {
+ 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;
+ }
+
+ 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)
+ switch (var_par->function) {
+ case 5:
+ r0[2] = reg->bot + (dep + 0.5) * reg->bt_res; // z0
+ break;
+ default:
+ r0[2] = xD->aniso_ratio * (reg->bot + (dep + 0.5) * reg->bt_res); // z0
+ break;
+ }
+ else
+ 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; // number of input points
+ double *r = pnts->r; // xyz coordinates of input points
+ int type = var_par->type;
+
+ int i, j; // indices of matrix rows/cols
+ int n1 = n+1; // number of matrix rows/cols
+ 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 *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++) { // for elements in each col 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) {
+ dr[0] = sqrt(SQUARE(dr[0]) + SQUARE(dr[1]));
+ dr[1] = dr[2];
+ }
+ teor_var = variogram_fction(var_par, dr);
+
+ if (isnan(teor_var)) {
+ 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) teor_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(std::vector<int> index, mat_struct *GM_all, struct write *report)
+{
+ // Local variables
+ int n = index.size ();
+ int *lind;
+ lind = &index.data ()[0];
+ mat_struct *GM = GM_all;
+
+ int i, j, k, N1 = GM->rows, n1 = n+1, *dinR0, *dinR, *dini, *dini0, *dinj0, *dinj;
+ doublereal *dbo, *dbx, *dbu, *dbl, *md, *mu, *ml, *m1r, *m1c;
+
+ mat_struct *sub;
+ 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 = &lind[0]; // i - set processing column
+ dini0 = &lind[0]; // i-1 - previous column
+ // initialize new col
+ dinR0 = &lind[1]; // return to first processing row - lower GM
+ dinR = &lind[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
+ dinj0 = dinR; // orig: inicialize element (row - 1) in (col)
+ // 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
+ //G_debug(0,"GM=%f", *dbx);
+ mu += n1; // go to element in next column
+ ml++; // go to element in next row
+ // Original matrix
+ dinj0 = dinj; // save current ind
+ dinj++; // set next ind
+ if (j < n-1)
+ dbx += *dinj - *dinj0;
+ } // end j
+ // Original matrix
+ dini0 = dini; // save current ind
+ dini++; // set next ind
+ if (i < n-1) {
+ dbu += n1+1; // new col in GM
+ dbl += n1+1; // new row in GM
+ }
+ dinR0 = dinR;
+ dinR++;
+ dbo += (*dini - *dini0) * N1 + (*dinR - *dinR0);
+ // 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, std::vector<int> ind, struct points *pnts, double *r0, struct parameters *var_par)
+{
+ // Local variables
+ int i3 = xD->i3; // 2D or 3D interpolation
+ int bivar = xD->bivar;
+ int type = var_par->type;
+ int n;
+ double *r; // xyz coordinates of input points
+
+ int *lind, *lind0;
+
+ if (ind.empty()) {
+ n = pnts->n;
+ r = &pnts->r[0];
+ }
+ else {
+ n = ind.size (); // number of input points
+ lind = &ind[0];
+ lind0 = &ind[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);
+
+ double *ro;
+ doublereal *g;
+ g = &g0->vals[0];
+
+ for (i=0; i<n; i++) { // count of input points
+ // Coord diffs (input points and cell/voxel center)
+ dr[0] = *r - r0[0]; // 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) {
+ dr[0] = sqrt(SQUARE(dr[0]) + SQUARE(dr[1]));
+ dr[1] = dr[2];
+ }
+
+ *g = variogram_fction(var_par, dr);
+ g++;
+
+ if (ind.empty())
+ r += 3;
+ else {
+ *lind0 = *lind;
+ lind++;
+ r += 3 * (*lind - *lind0);
+ }
+ } // 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(std::vector<int> ind, struct points *pnts, mat_struct *w0)
+{
+ int n;
+ double *vo;
+ mat_struct *loc_w = w0;
+
+ int *indo, *indo0, *lind;
+
+ if (ind.empty())
+ n = pnts->n;
+ else {
+ n = ind.size ();
+ lind = &ind[0];
+ }
+
+ int i;
+ mat_struct *ins, *w, *rslt_OK;
+ doublereal *vt, *wo, *wt;
+ ins = G_matrix_init(n, 1, n);
+ w = G_matrix_init(1, n, 1);
+
+ if (ind.empty())
+ vo = &pnts->invals[0];// original input values
+ else {
+ indo = &lind[0]; // original indexes
+ indo0 = &lind[0];
+ vo = &pnts->invals[*indo];// original input values
+ }
+
+ vt = &ins->vals[0]; // selected inputs values
+ wo = &loc_w->vals[0]; // selected inputs values
+ wt = &w->vals[0]; // weights without last one
+
+ for (i=0; i<n; i++) {
+ *vt = *vo;
+ *wt = *wo;
+ if (i<n-1) {
+ indo0 = indo;
+ indo++;
+ if (ind.empty())
+ vo++;
+ else
+ vo += *indo - *indo0;
+ vt++;
+ wo++;
+ wt++;
+ }
+ }
+ rslt_OK = G_matrix_product(w,ins);
+
+ G_matrix_free(w);
+ G_matrix_free(ins);
+
+ return rslt_OK->vals[0];
+}
+
+// validation
+void crossvalidation(struct int_par *xD, struct points *pnts, pcl::PointCloud<pcl::PointXYZ>::Ptr pcl_pnts, struct parameters *var_par)
+{
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // coordinates of points
+ double *vals = pnts->invals;
+ double ratio;
+ mat_struct *GM = var_par->GM;
+
+ ratio = var_par->type == 3 ? xD->aniso_ratio : 1.;
+
+ int i, j;
+ int n1;
+ double radius;
+ double r0[3];
+ mat_struct *GM_sub;
+ mat_struct *GM_Inv, *g0, *w0;
+ double rslt_OK;
+ pcl::KdTreeFLANN<pcl::PointXYZ> kd_tree;
+ kd_tree.setInputCloud(pcl_pnts); // Set up kd-tree
+ pcl::PointXYZ cellCent;
+ std::vector<int> ind;
+ std::vector<float> sqDist;
+ 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)
+ fp = fopen(crossvalid->name, "w");
+
+ radius = var_par->max_dist;
+
+ for (i=0; i<n; i++) {
+ // set up processing point
+ *r0 = cellCent.x = *r;
+ *(r0+1) = cellCent.y = *(r+1);
+ *(r0+2) = cellCent.z = *(r+2); // if 2D: 0.
+
+ if (kd_tree.radiusSearch(cellCent, radius, ind, sqDist) > 0 ) {
+ n1 = ind.size () + 1; // number of rellevant points + weights
+ GM_sub = submatrix(ind, GM, &xD->report);
+ GM_Inv = G_matrix_inverse(GM_sub);
+ G_matrix_free(GM_sub);
+
+ g0 = set_up_g0(xD, ind, pnts, 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(g0);
+ G_matrix_free(GM_Inv);
+ rslt_OK = result(ind, pnts, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
+
+ //Create output
+ *norm = *vals - rslt_OK;
+ *av = fabsf(*norm);
+ if (xD->i3 == TRUE)
+ fprintf(fp,"%d %.3f %.3f %.2f %f %f %f\n", i, *r, *(r+1), *(r+2) / ratio, pnts->invals[i], rslt_OK, *norm);
+ else
+ fprintf(fp,"%d %.3f %.3f %f %f %f\n", i, *r, *(r+1), *vals, rslt_OK, *norm);
+ }
+ else
+ fprintf(fp,"%d. point does not have neighbours in given radius\n", i);
+ r += 3;
+ vals++;
+ norm++;
+ av++;
+ }
+ 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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_raster.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_raster.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,61 @@
+#include "local_proto.h"
+
+extern "C" {
+ /* 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.cpp
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_write.cpp (rev 0)
+++ grass-addons/grass7/vector/v.kriging/utils_write.cpp 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1,305 @@
+#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_par)
+{
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ fprintf(report->fp, "- number of lags in %s direction: %d\n", type, var_par->nLag);
+ fprintf(report->fp, "- lag distance (%s): %f\n", type, var_par->lag);
+
+ if (var_par->function == 5) {
+ fprintf(report->fp, "- number of lags in vertical direction: %d\n", var_par->nLag_vert);
+ fprintf(report->fp, "- lag distance (vertical): %f\n", var_par->lag);
+ }
+
+ fprintf(report->fp, "- azimuth: %f°\n", RAD2DEG(var_par->dir));
+ fprintf(report->fp, "- angular tolerance: %f°\n", RAD2DEG(var_par->td));
+ //fprintf(report->fp, "- variogram plotted: %s\n", );
+
+ switch (var_par->function) {
+ case 5:
+ fprintf(report->fp, "- maximum distance in horizontal direction (1/3): %f\n", var_par->max_dist);
+ break;
+ default:
+ fprintf(report->fp, "- maximum distance (1/3): %f\n", var_par->max_dist);
+ break;
+ }
+}
+
+void write2file_variogram_E(struct int_par *xD, struct parameters *var_par)
+{
+ int i, j;
+ int fction = var_par->function;
+ double *h, *vert;
+ double *c = var_par->c->vals;
+ double *gamma = var_par->gamma->vals;
+ struct write *report = &xD->report;
+
+ char type[12];
+ variogram_type(var_par->type, type);
+
+ fprintf(report->fp, "\n************************************************\n");
+ fprintf(report->fp, "*** Experimental variogram - %s ***\n\n", type);
+
+ if (fction == 5) { // if variogram is bivariate
+ vert = &var_par->vert[0]; // use vertical variable
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - verts
+ if (i == 0)
+ fprintf(report->fp, " lagV ||"); // column for h
+ fprintf(report->fp, " %f ||", *vert);
+ //fprintf(report->fp, " # of pairs | average ||");
+ //fprintf(report->fp, "------------------------");
+ vert++;
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - h c gamma
+ if (i == 0)
+ fprintf(report->fp, " lagHZ ||"); // column for h
+ fprintf(report->fp, " c | ave ||", *vert);
+ }
+ fprintf(report->fp, "\n");
+
+ for (i=0; i<var_par->nLag_vert; i++) { // write header - h c gamma
+ if (i == 0)
+ fprintf(report->fp, "--------||"); // column for h
+ fprintf(report->fp, "-----------", *vert);
+ }
+ fprintf(report->fp, "\n");
+ }
+ else {
+ fprintf(report->fp, " lag || # of pairs | average\n");
+ fprintf(report->fp, "------------------------------------\n");
+ }
+
+ // Write values
+ h = &var_par->h[0];
+ for (i=0; i<var_par->nLag; i++) {
+ fprintf(report->fp, "%f ||", *h);
+ if (fction == 5) {
+ //for (j=0; j<var_par->nZ; j++) {
+ fprintf(report->fp, " %d | %f ||", (int) *c, *gamma);
+ c++;
+ gamma++;
+ //} // end for j
+ fprintf(report->fp, "\n");
+ //for (j=0; j<var_par->nZ; j++)
+ fprintf(report->fp, "-----------", (int) *c, *gamma);
+ fprintf(report->fp, "\n");
+ }
+ else {
+ 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_par, mat_struct *gammaMT)
+{
+ // local variables
+ int type = var_par->type;
+ int nLag = var_par->nLag;
+ int nLag_vert;
+ double *h = var_par->h;
+ double *vert;
+ double *c = var_par->c->vals;
+ double *gamma = gammaMT->vals;
+ double sill = var_par->sill;
+
+ int i, j; // index
+ FILE *fp;
+ int file_length;
+
+ 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_par->type); // write # of lags
+ break;
+ case 1: // vertical variogram
+ fp = fopen("variogram_vert_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 3)
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp,"%d\n", var_par->type); // write type
+ break;
+ case 2:
+ 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 2 characters...")); // todo: error
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+
+ fprintf(fp,"%d\n", var_par->type); // write # of lags
+ fprintf(fp,"%d\n", var_par->nLag_vert); // write # of lags
+ fprintf(fp,"%f\n", var_par->lag_vert); // write size of lag
+ fprintf(fp,"%f\n", var_par->max_dist_vert); // 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_par->type); // write type
+ fprintf(fp,"%f\n", xD->aniso_ratio); // write ratio of anisotropy
+ break;
+ }
+
+ fprintf(fp,"%d\n", nLag); // write # of lags
+ fprintf(fp,"%f\n", var_par->lag); // write size of lag
+ fprintf(fp,"%f\n", var_par->max_dist); // write maximum distance
+ if (type != 1)
+ fprintf(fp,"%f\n", var_par->td); // write maximum distance
+
+ switch (type) {
+ case 2:
+ nLag_vert = var_par->nLag_vert;
+ // write h
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *h);
+ h++;
+ }
+ // write vert
+ vert = &var_par->vert[0];
+ for (i=0; i<nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *vert);
+ vert++;
+ }
+ // write c
+ for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *c);
+ c++;
+ }
+ // write gamma
+ for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp,"%f\n", *gamma);
+ gamma++;
+ }
+ fprintf(fp,"%f\n", var_par->horizontal.sill);// write sill
+ fprintf(fp,"%f\n", var_par->vertical.sill);// write sill
+ break;
+ default:
+ for (i=0; i<nLag; i++) { // write experimental variogram
+ fprintf(fp,"%f %f %f\n", *h, *c, *gamma);
+ h++;
+ c++;
+ gamma++;
+ }
+ fprintf(fp,"%f", sill);// write sill
+ break;
+ }
+
+ fclose(fp);
+}
Added: grass-addons/grass7/vector/v.kriging/v.kriging.html
===================================================================
--- grass-addons/grass7/vector/v.kriging/v.kriging.html (rev 0)
+++ grass-addons/grass7/vector/v.kriging/v.kriging.html 2014-10-21 13:21:36 UTC (rev 62319)
@@ -0,0 +1 @@
+TO DO
More information about the grass-commit
mailing list