[GRASS-SVN] r65551 - grass-addons/grass7/vector/v.kriging
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jul 9 02:17:29 PDT 2015
Author: evas
Date: 2015-07-09 02:17:29 -0700 (Thu, 09 Jul 2015)
New Revision: 65551
Modified:
grass-addons/grass7/vector/v.kriging/Makefile
grass-addons/grass7/vector/v.kriging/geostat.c
grass-addons/grass7/vector/v.kriging/getval.c
grass-addons/grass7/vector/v.kriging/main.c
grass-addons/grass7/vector/v.kriging/quantile.c
grass-addons/grass7/vector/v.kriging/spatial_index.c
grass-addons/grass7/vector/v.kriging/stat.c
grass-addons/grass7/vector/v.kriging/utils.c
grass-addons/grass7/vector/v.kriging/utils_kriging.c
grass-addons/grass7/vector/v.kriging/utils_raster.c
grass-addons/grass7/vector/v.kriging/utils_write.c
Log:
v.kriging: indented
Modified: grass-addons/grass7/vector/v.kriging/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/Makefile 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/Makefile 2015-07-09 09:17:29 UTC (rev 65551)
@@ -1,4 +1,4 @@
-MODULE_TOPDIR = ../..
+MODULE_TOPDIR = /home/evergreen/grass7
PGM = v.kriging
Modified: grass-addons/grass7/vector/v.kriging/geostat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/geostat.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/geostat.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -2,618 +2,650 @@
/* experimental variogram
* based on 2D variogram (alghalandis.com/?page_id=463)) */
- void E_variogram(int type, struct int_par *xD, struct points *pnts, struct var_par *pars)
+void E_variogram(int type, struct int_par *xD, struct points *pnts,
+ struct var_par *pars)
{
- // Variogram properties
- struct parameters *var_pars;
- double max_dist, max_dist_vert; // max distance of interpolated points
- double radius, radius_vert; // radius of interpolation in the horizontal plane
+ // Variogram properties
+ struct parameters *var_pars;
+ double max_dist, max_dist_vert; // max distance of interpolated points
+ double radius, radius_vert; // radius of interpolation in the horizontal plane
- switch (type) {
- case 0: // horizontal variogram
- var_pars = &pars->hz; // 0: horizontal variogram
- break;
- case 1: // vertica variogram
- var_pars = &pars->vert; // 1: vertical variogram
- break;
- case 2: // bivariate variogram
- var_pars = &pars->fin;
+ switch (type) {
+ case 0: // horizontal variogram
+ var_pars = &pars->hz; // 0: horizontal variogram
+ break;
+ case 1: // vertica variogram
+ var_pars = &pars->vert; // 1: vertical variogram
+ break;
+ case 2: // bivariate variogram
+ var_pars = &pars->fin;
- max_dist = var_pars->horizontal.max_dist;
- max_dist_vert = var_pars->vertical.max_dist;
- radius_vert = SQUARE(max_dist_vert);
- break;
- case 3: //anisotropic variogram
- var_pars = &pars->fin; // default: final variogram
- max_dist = max_dist_vert = var_pars->max_dist;
- radius = radius_vert = SQUARE(max_dist);
- break;
- }
+ max_dist = var_pars->horizontal.max_dist;
+ max_dist_vert = var_pars->vertical.max_dist;
+ radius_vert = SQUARE(max_dist_vert);
+ break;
+ case 3: //anisotropic variogram
+ var_pars = &pars->fin; // default: final variogram
+ max_dist = max_dist_vert = var_pars->max_dist;
+ radius = radius_vert = SQUARE(max_dist);
+ break;
+ }
- if (type < 2) {
- max_dist = var_pars->max_dist;
- }
- radius = SQUARE(max_dist);
+ if (type < 2) {
+ max_dist = var_pars->max_dist;
+ }
+ radius = SQUARE(max_dist);
- // Local variables
- int n = pnts->n; // # of input points
- double *pnts_r = pnts->r; // xyz coordinates of input points
- double *r; // pointer to xyz coordinates
- double *search; // pointer to search point coordinates
- double *vals = pnts->invals; // values to be used for interpolation
- int phase = xD->phase; // phase: initial / middle / final
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *pnts_r = pnts->r; // xyz coordinates of input points
+ double *r; // pointer to xyz coordinates
+ double *search; // pointer to search point coordinates
+ double *vals = pnts->invals; // values to be used for interpolation
+ int phase = xD->phase; // phase: initial / middle / final
- int nLag; // # of horizontal bins
- int nLag_vert; // # of vertical bins
- double *vert; // pointer to vertical bins
- double dir = var_pars->dir; // azimuth
- double td = var_pars->td; // angle tolerance
- double lag = var_pars->lag; // size of horizontal bin
- double lag_vert;
+ int nLag; // # of horizontal bins
+ int nLag_vert; // # of vertical bins
+ double *vert; // pointer to vertical bins
+ double dir = var_pars->dir; // azimuth
+ double td = var_pars->td; // angle tolerance
+ double lag = var_pars->lag; // size of horizontal bin
+ double lag_vert;
- if (type == 2) { // bivariate variogram
- nLag = var_pars->horizontal.nLag;
- nLag_vert = var_pars->vertical.nLag;
+ if (type == 2) { // bivariate variogram
+ nLag = var_pars->horizontal.nLag;
+ nLag_vert = var_pars->vertical.nLag;
- dir = var_pars->dir;
- td = var_pars->td;
- lag = var_pars->horizontal.lag;
- lag_vert = var_pars->vertical.lag;
- }
+ dir = var_pars->dir;
+ td = var_pars->td;
+ lag = var_pars->horizontal.lag;
+ lag_vert = var_pars->vertical.lag;
+ }
- else { // univariate variogram
- nLag = var_pars->nLag;
- dir = var_pars->dir;
- td = var_pars->td;
- lag = var_pars->lag;
- }
+ else { // univariate variogram
+ nLag = var_pars->nLag;
+ dir = var_pars->dir;
+ td = var_pars->td;
+ lag = var_pars->lag;
+ }
- // depend on variogram type:
- int nrows = nLag;
- int ncols = type == 2 ? nLag_vert : 1;
+ // depend on variogram type:
+ int nrows = nLag;
+ int ncols = type == 2 ? nLag_vert : 1;
- struct write *report = &xD->report;
-
- // Variogram processing variables
- int s; // index of horizontal segment
- int b; // index of verical segment (bivariate only)
- int i, j; // indices of the input points
- int err0=0; // number of invalid values
+ struct write *report = &xD->report;
- 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 ddir1, ddir2; // the azimuth of computing variogram
- double rv; // radius of the couple of 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
+ // 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
- struct ilist *list; // list of selected nearest neighbours
- int n_vals; // # of selected NN
+ 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 ddir1, ddir2; // the azimuth of computing variogram
+ double rv; // radius of the couple of 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
+ struct ilist *list; // list of selected nearest neighbours
+ int n_vals; // # of selected NN
- int *ii; // difference of indices between rellevant input points
-
- double gamma_lag; // sum of dissimilarities in one bin
- double cpls; // # of dissimilarities in one bin
- mat_struct *gamma_M; // gamma matrix (hz, vert or bivar)
- mat_struct *c_M; // matrix of # of dissimilarities
- double *gamma; // pointer to gamma matrix
- double *c; // pointer to c matrix
+ double *i_vals, *j_vals; // values located on the point couples
- unsigned int percents = 50;
+ int *ii; // difference of indices between rellevant input points
- double gamma_sum; // sum of gamma elements (non-nan)
- int gamma_n; // # of gamma elements (non-nan)
-
- /* 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
- */
+ double gamma_lag; // sum of dissimilarities in one bin
+ double cpls; // # of dissimilarities in one bin
+ mat_struct *gamma_M; // gamma matrix (hz, vert or bivar)
+ mat_struct *c_M; // matrix of # of dissimilarities
+ double *gamma; // pointer to gamma matrix
+ double *c; // pointer to c matrix
- // allocate:
- dr = (double *) G_malloc(3 * sizeof(double)); // vector of coordinate differences
- search = (double *) G_malloc(3 * sizeof(double)); // vector of search point coordinates
+ unsigned int percents = 50;
- if (type != 2) {
- var_pars->h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
- }
+ double gamma_sum; // sum of gamma elements (non-nan)
+ int gamma_n; // # of gamma elements (non-nan)
- if (type == 2) {
- var_pars->horizontal.h = (double *) G_malloc(nLag * sizeof(double)); // vector of bins
- var_pars->vertical.h = (double *) G_malloc(nLag_vert * sizeof(double)); // vector of bins
- }
+ /* 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
+ */
- // control initialization:
- if (dr == NULL || search == NULL || var_pars->h == NULL || (type == 2 && var_pars->vert == NULL)) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp, "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ // allocate:
+ dr = (double *)G_malloc(3 * sizeof(double)); // vector of coordinate differences
+ search = (double *)G_malloc(3 * sizeof(double)); // vector of search point coordinates
+
+ if (type != 2) {
+ var_pars->h = (double *)G_malloc(nLag * sizeof(double)); // vector of bins
}
- G_fatal_error(_("Memory allocation failed..."));
- }
-
- // initialize:
- c_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix of counts
- gamma_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix (vector) of gammas
- if (c_M == NULL || gamma_M == NULL) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp, "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ if (type == 2) {
+ var_pars->horizontal.h = (double *)G_malloc(nLag * sizeof(double)); // vector of bins
+ var_pars->vertical.h = (double *)G_malloc(nLag_vert * sizeof(double)); // vector of bins
}
- G_fatal_error(_("Memory allocation failed..."));
- }
-
- // set up pointers
- c = &c_M->vals[0];
- gamma = &gamma_M->vals[0];
- // set up starting values
- gamma_sum = 0.;
- gamma_n = 0;
+ // control initialization:
+ if (dr == NULL || search == NULL || var_pars->h == NULL ||
+ (type == 2 && var_pars->vert == NULL)) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
+ }
- if (percents) {
- G_percent_reset();
- }
+ // initialize:
+ c_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix of counts
+ gamma_M = G_matrix_init(nrows, ncols, nrows); // temporal matrix (vector) of gammas
- if (type == 2) {
- vert = &var_pars->vertical.h[0];
- }
-
- /* *** Experimental variogram computation *** */
- for (b = 0; b < ncols; b++) { // for each vertical lag...
- if (type == 2) { // just in case that the variogram is vertical
- *vert = (b+1) * lag_vert; // lag distance
+ if (c_M == NULL || gamma_M == NULL) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Memory allocation failed..."));
}
-
- h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0]; // ... horizontal lags
- for (s = 0; s < nrows; s++) { // for each horizontal lag (isotrophy!!!)...
- *h = (s+1) * lag; // lag distance
- r = &pnts->r[0]; // pointer to the input coordinates
-
- // for every bs cycle ...
- i_vals = &vals[0]; // pointer to input values to be interpolated
- gamma_lag = 0.; // gamma in dir direction and h distance
- cpls = 0.; // count of couples in dir direction and h distance
+ // set up pointers
+ c = &c_M->vals[0];
+ gamma = &gamma_M->vals[0];
- /* Compute variogram for points in relevant neighbourhood */
- for (i = 0; i < n-1; i++) { // for each input point...
- search = &pnts->r[3 * i];
- switch (type) { // find NNs according to variogram type
- case 0: // horizontal variogram
- list = find_NNs_within(2, search, pnts, max_dist, -1);
- break;
- case 1: // vertical variogram
- list = find_NNs_within(1, search, pnts, -1, max_dist);
- break;
- default: // anisotropic or bivariate variogram
- list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
- break;
- }
-
- n_vals = list->n_values; // # of input values located on NN
- if (n_vals > 0) {
- correct_indices(list, r, pnts, var_pars);
- ii = &list->value[0]; // indices of these input values (note: increased by 1)
- j_vals = &vals[*ii]; // pointer to input values
-
- for (j = 1; j < n_vals; j++) { // for each point within overlapping rectangle
- if (*ii > i) { // use the points just once
- coord_diff(i, *ii, pnts_r, dr); // compute coordinate differences
-
- // Variogram processing
- if (type == 1) { // vertical variogram:
- tv = 0.5 * PI - zenith_angle(dr); // compute zenith angle instead of bearing
- td = 0.5 * PI; // todo: repair, 0.25 usually
- }
- else { // hz / aniso / bivar variogram:
- tv = atan2(*(dr+1), *dr); // bearing
- if (tv < 0.) {
- tv += 2. * PI;
- }
- }
-
- ddir1 = dir - tv; // difference between bearing and azimuth
- ddir2 = (dir + PI) - tv;
-
- if (fabsf(ddir1) <= td || fabsf(ddir2) <= td) { // angle test: compare the diff with critical value
- // test squared distance: vertical variogram => 0., ...
- rv = type == 1 ? 0. : radius_hz_diff(dr); // ... otherwise horizontal distance
+ // set up starting values
+ gamma_sum = 0.;
+ gamma_n = 0;
- if (type == 1 || type == 3) { // vertical or anisotropic variogram:
- rv += SQUARE(*(dr+2)); // consider also vertical direction
- }
+ if (percents) {
+ G_percent_reset();
+ }
- rvh = sqrt(rv) - *h; // the difference between distance and lag
- if (rv <= radius && fabsf(rvh) <= lag) { // distance test: compare the distance with critical value and find out if the j-point is located within i-lag
- if (type == 2) { // vertical test for bivariate variogram:
- rvh = *(dr+2) - *vert; // compare vertical
-
- if (fabsf(rvh) <= lag_vert) { // elevation test: vertical lag
- goto delta_V;
- }
- else {
- goto end;
- }
- }
- delta_V:
- dv = *j_vals - *i_vals; // difference of values located on pair of points i, j
- gamma_lag += SQUARE(dv); // sum of squared differences
- cpls += 1.; // count of elements
- } // end distance test: rv <= radius ^ |rvh| <= lag
- } // end angle test: |ddir| <= td
- } // end i test: *ii > i
- end:
- ii++; // go to the next index
- j_vals += *ii - *(ii-1); // go to the next value
- } // end j for loop: points within overlapping rectangles
- } // end test: n_vals > 0
- else {
- if (report->name) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("This point does not have neighbours in given radius..."));
- }
+ if (type == 2) {
+ vert = &var_pars->vertical.h[0];
+ }
- r += 3; // go to the next search point
- i_vals++; // and to its value
- } // end for loop i: variogram computation for each i-th search point
-
- if (isnan(gamma_lag) || cpls == 0.0) { // empty lags:
- err0++; // error indicator
- goto gamma_isnan;
- }
+ /* *** Experimental variogram computation *** */
+ for (b = 0; b < ncols; b++) { // for each vertical lag...
+ if (type == 2) { // just in case that the variogram is vertical
+ *vert = (b + 1) * lag_vert; // lag distance
+ }
- *gamma = 0.5 * gamma_lag / cpls; // element of gamma matrix
- *c = cpls; // # of values that have been used to compute gamma_e
-
- gamma_sum += *gamma; // sum of gamma elements (non-nan)
- gamma_n++; // # of gamma elements (non-nan)
+ h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0]; // ... horizontal lags
- gamma_isnan: // there are no available point couples to compute the dissimilarities in the lag:
- h++;
- c++;
- gamma++;
- } // end for loop s
+ for (s = 0; s < nrows; s++) { // for each horizontal lag (isotrophy!!!)...
+ *h = (s + 1) * lag; // lag distance
+ r = &pnts->r[0]; // pointer to the input coordinates
- if (type == 2) { // vertical variogram:
- vert++;
- }
- } // end for loop b
+ // for every bs cycle ...
+ i_vals = &vals[0]; // pointer to input values to be interpolated
+ gamma_lag = 0.; // gamma in dir direction and h distance
+ cpls = 0.; // count of couples in dir direction and h distance
- 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
+ /* Compute variogram for points in relevant neighbourhood */
+ for (i = 0; i < n - 1; i++) { // for each input point...
+ search = &pnts->r[3 * i];
+ switch (type) { // find NNs according to variogram type
+ case 0: // horizontal variogram
+ list = find_NNs_within(2, search, pnts, max_dist, -1);
+ break;
+ case 1: // vertical variogram
+ list = find_NNs_within(1, search, pnts, -1, max_dist);
+ break;
+ default: // anisotropic or bivariate variogram
+ list =
+ find_NNs_within(3, search, pnts, max_dist,
+ max_dist_vert);
+ break;
+ }
- var_pars->gamma = G_matrix_copy(gamma_M);
- var_pars->gamma_sum = gamma_sum;
- var_pars->gamma_n = gamma_n;
+ n_vals = list->n_values; // # of input values located on NN
+ if (n_vals > 0) {
+ correct_indices(list, r, pnts, var_pars);
+ ii = &list->value[0]; // indices of these input values (note: increased by 1)
+ j_vals = &vals[*ii]; // pointer to input values
- plot_experimental_variogram(xD, var_pars);
+ for (j = 1; j < n_vals; j++) { // for each point within overlapping rectangle
+ if (*ii > i) { // use the points just once
+ coord_diff(i, *ii, pnts_r, dr); // compute coordinate differences
- if (phase < 2) { // initial and middle phase:
- sill(var_pars); // compute sill
- }
+ // Variogram processing
+ if (type == 1) { // vertical variogram:
+ tv = 0.5 * PI - zenith_angle(dr); // compute zenith angle instead of bearing
+ td = 0.5 * PI; // todo: repair, 0.25 usually
+ }
+ else { // hz / aniso / bivar variogram:
+ tv = atan2(*(dr + 1), *dr); // bearing
+ if (tv < 0.) {
+ tv += 2. * PI;
+ }
+ }
- if (report->write2file == TRUE) { // report file available:
- write2file_variogram_E(xD, var_pars, c_M); // write to file
- }
+ ddir1 = dir - tv; // difference between bearing and azimuth
+ ddir2 = (dir + PI) - tv;
- write_temporary2file(xD, var_pars);
-
- G_free_ilist(list); // free list memory
- G_matrix_free(c_M);
- G_matrix_free(gamma_M);
-}
+ if (fabsf(ddir1) <= td || fabsf(ddir2) <= td) { // angle test: compare the diff with critical value
+ // test squared distance: vertical variogram => 0., ...
+ rv = type == 1 ? 0. : radius_hz_diff(dr); // ... otherwise horizontal distance
-/* theoretical variogram */
-void T_variogram(int type, int i3, struct opts opt, struct parameters *var_pars, struct write *report)
-{
- char *variogram;
+ if (type == 1 || type == 3) { // vertical or anisotropic variogram:
+ rv += SQUARE(*(dr + 2)); // consider also vertical direction
+ }
- // report
- if (report->write2file == TRUE) { // report file available:
- time(&report->now); // write down time of start
- if (type != 1) {
- fprintf(report->fp, "\nComputation of theoretical variogram started on %s\n", ctime(&report->now));
- }
- }
+ rvh = sqrt(rv) - *h; // the difference between distance and lag
+ if (rv <= radius && fabsf(rvh) <= lag) { // distance test: compare the distance with critical value and find out if the j-point is located within i-lag
+ if (type == 2) { // vertical test for bivariate variogram:
+ rvh = *(dr + 2) - *vert; // compare vertical
- // set up:
- var_pars->type = type; // hz / vert / bivar / aniso
- var_pars->const_val = 0; // input values are not constants
+ if (fabsf(rvh) <= lag_vert) { // elevation test: vertical lag
+ goto delta_V;
+ }
+ else {
+ goto end;
+ }
+ }
+ delta_V:
+ dv = *j_vals - *i_vals; // difference of values located on pair of points i, j
+ gamma_lag += SQUARE(dv); // sum of squared differences
+ cpls += 1.; // count of elements
+ } // end distance test: rv <= radius ^ |rvh| <= lag
+ } // end angle test: |ddir| <= td
+ } // end i test: *ii > i
+ end:
+ ii++; // go to the next index
+ j_vals += *ii - *(ii - 1); // go to the next value
+ } // end j for loop: points within overlapping rectangles
+ } // end test: n_vals > 0
+ else {
+ if (report->name) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("This point does not have neighbours in given radius..."));
+ }
- switch (type) {
- case 0: // horizontal variogram
- if (i3 == TRUE) { // 3D interpolation (middle phase)
- variogram = opt.function_var_hz->answer; // function type available:
- var_pars->function = set_function(variogram, report);
+ r += 3; // go to the next search point
+ i_vals++; // and to its value
+ } // end for loop i: variogram computation for each i-th search point
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
- var_pars->nugget = atof(opt.nugget_hz->answer);
- var_pars->h_range = atof(opt.range_hz->answer);
- if (opt.sill_hz->answer) {
- var_pars->sill = atof(opt.sill_hz->answer);
- }
- }
- else { // function type not available:
- LMS_variogram(var_pars, report);
- }
- }
- else { // 2D interpolation (final phase)
- variogram = opt.function_var_final->answer; // function type available:
- var_pars->function = set_function(variogram, report);
+ if (isnan(gamma_lag) || cpls == 0.0) { // empty lags:
+ err0++; // error indicator
+ goto gamma_isnan;
+ }
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
- var_pars->nugget = atof(opt.nugget_final->answer);
- var_pars->h_range = atof(opt.range_final->answer);
- if (opt.sill_final->answer) {
- var_pars->sill = atof(opt.sill_final->answer);
- }
- }
- else { // function type not available:
- LMS_variogram(var_pars, report);
- }
- }
-
- if (report->name) {
- fprintf(report->fp, "Parameters of horizontal variogram:\n");
- fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
- fprintf(report->fp, "Sill: %f\n", var_pars->sill);
- fprintf(report->fp, "Range: %f\n", var_pars->h_range);
- }
- break;
+ *gamma = 0.5 * gamma_lag / cpls; // element of gamma matrix
+ *c = cpls; // # of values that have been used to compute gamma_e
- case 1: // vertical variogram
- var_pars->nugget = atof(opt.nugget_vert->answer);
- var_pars->h_range = atof(opt.range_vert->answer);
- if (opt.sill_vert->answer) {
- var_pars->sill = atof(opt.sill_vert->answer);
+ gamma_sum += *gamma; // sum of gamma elements (non-nan)
+ gamma_n++; // # of gamma elements (non-nan)
+
+ gamma_isnan: // there are no available point couples to compute the dissimilarities in the lag:
+ h++;
+ c++;
+ gamma++;
+ } // end for loop s
+
+ if (type == 2) { // vertical variogram:
+ vert++;
+ }
+ } // end for loop b
+
+ if (err0 == nLag) { // todo> kedy nie je riesitelny teoreticky variogram?
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to compute experimental variogram..."));
+ } // end error
+
+ var_pars->gamma = G_matrix_copy(gamma_M);
+ var_pars->gamma_sum = gamma_sum;
+ var_pars->gamma_n = gamma_n;
+
+ plot_experimental_variogram(xD, var_pars);
+
+ if (phase < 2) { // initial and middle phase:
+ sill(var_pars); // compute sill
}
- variogram = opt.function_var_vert->answer;
- var_pars->function = set_function(variogram, report);
- if (report->name) {
- fprintf(report->fp, "Parameters of vertical variogram:\n");
- fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
- fprintf(report->fp, "Sill: %f\n", var_pars->sill);
- fprintf(report->fp, "Range: %f\n", var_pars->h_range);
+ if (report->write2file == TRUE) { // report file available:
+ write2file_variogram_E(xD, var_pars, c_M); // write to file
}
- break;
- case 2: // bivariate variogram (just final phase)
- if (!(opt.function_var_final->answer && opt.function_var_final_vert->answer) || strcmp(opt.function_var_final->answer, "linear") == 0) { // planar function:
- var_pars->function = 5; // planar variogram (3D)
- LMS_variogram(var_pars, report);
- }
- else { // function type was set up by the user:
- var_pars->horizontal.nugget = atof(opt.nugget_final->answer);
- var_pars->horizontal.h_range = atof(opt.range_final->answer);
- if (opt.sill_final->answer) {
- var_pars->horizontal.sill = atof(opt.sill_final->answer);
- }
+ write_temporary2file(xD, var_pars);
- var_pars->vertical.nugget = atof(opt.nugget_final_vert->answer);
- var_pars->vertical.h_range = atof(opt.range_final_vert->answer);
- if (opt.sill_final_vert->answer) {
- var_pars->vertical.sill = atof(opt.sill_final_vert->answer);
- }
-
- var_pars->horizontal.function = set_function(opt.function_var_final->answer, report);
- var_pars->vertical.function = set_function(opt.function_var_final_vert->answer, report);
+ G_free_ilist(list); // free list memory
+ G_matrix_free(c_M);
+ G_matrix_free(gamma_M);
+}
- if (report->name) {
- fprintf(report->fp, "Parameters of bivariate variogram:\n");
- fprintf(report->fp, "Nugget effect (hz): %f\n", var_pars->horizontal.nugget);
- fprintf(report->fp, "Sill (hz): %f\n", var_pars->horizontal.sill);
- fprintf(report->fp, "Range (hz): %f\n", var_pars->horizontal.h_range);
- fprintf(report->fp, "Function: %s\n\n", opt.function_var_hz->answer);
- fprintf(report->fp, "Nugget effect (vert): %f\n", var_pars->vertical.nugget);
- fprintf(report->fp, "Sill (vert): %f\n", var_pars->vertical.sill);
- fprintf(report->fp, "Range (vert): %f\n", var_pars->vertical.h_range);
- fprintf(report->fp, "Function: %s\n", opt.function_var_vert->answer);
- }
+/* theoretical variogram */
+void T_variogram(int type, int i3, struct opts opt,
+ struct parameters *var_pars, struct write *report)
+{
+ char *variogram;
+
+ // report
+ if (report->write2file == TRUE) { // report file available:
+ time(&report->now); // write down time of start
+ if (type != 1) {
+ fprintf(report->fp,
+ "\nComputation of theoretical variogram started on %s\n",
+ ctime(&report->now));
+ }
}
- plot_var(i3, TRUE, var_pars); // Plot variogram using gnuplot
- break;
+ // set up:
+ var_pars->type = type; // hz / vert / bivar / aniso
+ var_pars->const_val = 0; // input values are not constants
- case 3: // univariate (just final phase)
- variogram = opt.function_var_final->answer;
- var_pars->function = set_function(variogram, report);
+ switch (type) {
+ case 0: // horizontal variogram
+ if (i3 == TRUE) { // 3D interpolation (middle phase)
+ variogram = opt.function_var_hz->answer; // function type available:
+ var_pars->function = set_function(variogram, report);
- if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear and not parabolic variogram:
- var_pars->nugget = atof(opt.nugget_final->answer);
- var_pars->h_range = atof(opt.range_final->answer);
- if (opt.sill_final->answer) {
- var_pars->sill = atof(opt.sill_final->answer);
- }
- variogram = opt.function_var_final->answer;
- if (report->write2file == TRUE) {
- if (i3 == TRUE) { // 3D interpolation:
- fprintf(report->fp, "Parameters of anisotropic variogram:\n");
- }
- else {
- fprintf(report->fp, "Parameters of 2D variogram:\n");
- }
- fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
- fprintf(report->fp, "Sill: %f\n", var_pars->sill);
- fprintf(report->fp, "Range: %f\n", var_pars->h_range);
- }
+ if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
+ var_pars->nugget = atof(opt.nugget_hz->answer);
+ var_pars->h_range = atof(opt.range_hz->answer);
+ if (opt.sill_hz->answer) {
+ var_pars->sill = atof(opt.sill_hz->answer);
+ }
+ }
+ else { // function type not available:
+ LMS_variogram(var_pars, report);
+ }
+ }
+ else { // 2D interpolation (final phase)
+ variogram = opt.function_var_final->answer; // function type available:
+ var_pars->function = set_function(variogram, report);
+
+ if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear or not parabolic
+ var_pars->nugget = atof(opt.nugget_final->answer);
+ var_pars->h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_pars->sill = atof(opt.sill_final->answer);
+ }
+ }
+ else { // function type not available:
+ LMS_variogram(var_pars, report);
+ }
+ }
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of horizontal variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_pars->sill);
+ fprintf(report->fp, "Range: %f\n", var_pars->h_range);
+ }
+ break;
+
+ case 1: // vertical variogram
+ var_pars->nugget = atof(opt.nugget_vert->answer);
+ var_pars->h_range = atof(opt.range_vert->answer);
+ if (opt.sill_vert->answer) {
+ var_pars->sill = atof(opt.sill_vert->answer);
+ }
+ variogram = opt.function_var_vert->answer;
+ var_pars->function = set_function(variogram, report);
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of vertical variogram:\n");
+ fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_pars->sill);
+ fprintf(report->fp, "Range: %f\n", var_pars->h_range);
+ }
+ break;
+ case 2: // bivariate variogram (just final phase)
+ if (!(opt.function_var_final->answer && opt.function_var_final_vert->answer) || strcmp(opt.function_var_final->answer, "linear") == 0) { // planar function:
+ var_pars->function = 5; // planar variogram (3D)
+ LMS_variogram(var_pars, report);
+ }
+
+ else { // function type was set up by the user:
+ var_pars->horizontal.nugget = atof(opt.nugget_final->answer);
+ var_pars->horizontal.h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_pars->horizontal.sill = atof(opt.sill_final->answer);
+ }
+
+ var_pars->vertical.nugget = atof(opt.nugget_final_vert->answer);
+ var_pars->vertical.h_range = atof(opt.range_final_vert->answer);
+ if (opt.sill_final_vert->answer) {
+ var_pars->vertical.sill = atof(opt.sill_final_vert->answer);
+ }
+
+ var_pars->horizontal.function =
+ set_function(opt.function_var_final->answer, report);
+ var_pars->vertical.function =
+ set_function(opt.function_var_final_vert->answer, report);
+
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Nugget effect (hz): %f\n",
+ var_pars->horizontal.nugget);
+ fprintf(report->fp, "Sill (hz): %f\n",
+ var_pars->horizontal.sill);
+ fprintf(report->fp, "Range (hz): %f\n",
+ var_pars->horizontal.h_range);
+ fprintf(report->fp, "Function: %s\n\n",
+ opt.function_var_hz->answer);
+ fprintf(report->fp, "Nugget effect (vert): %f\n",
+ var_pars->vertical.nugget);
+ fprintf(report->fp, "Sill (vert): %f\n",
+ var_pars->vertical.sill);
+ fprintf(report->fp, "Range (vert): %f\n",
+ var_pars->vertical.h_range);
+ fprintf(report->fp, "Function: %s\n",
+ opt.function_var_vert->answer);
+ }
+ }
+
+ plot_var(i3, TRUE, var_pars); // Plot variogram using gnuplot
+ break;
+
+ case 3: // univariate (just final phase)
+ variogram = opt.function_var_final->answer;
+ var_pars->function = set_function(variogram, report);
+
+ if (strcmp(variogram, "linear") != 0 && strcmp(variogram, "parabolic") != 0) { // nonlinear and not parabolic variogram:
+ var_pars->nugget = atof(opt.nugget_final->answer);
+ var_pars->h_range = atof(opt.range_final->answer);
+ if (opt.sill_final->answer) {
+ var_pars->sill = atof(opt.sill_final->answer);
+ }
+ variogram = opt.function_var_final->answer;
+ if (report->write2file == TRUE) {
+ if (i3 == TRUE) { // 3D interpolation:
+ fprintf(report->fp,
+ "Parameters of anisotropic variogram:\n");
+ }
+ else {
+ fprintf(report->fp, "Parameters of 2D variogram:\n");
+ }
+ fprintf(report->fp, "Nugget effect: %f\n", var_pars->nugget);
+ fprintf(report->fp, "Sill: %f\n", var_pars->sill);
+ fprintf(report->fp, "Range: %f\n", var_pars->h_range);
+ }
+ }
+ else { // linear or parabolic variogram
+ LMS_variogram(var_pars, report);
+ }
+ break;
}
- else { // linear or parabolic variogram
- LMS_variogram(var_pars, report);
+
+ if (type != 2) {
+ plot_var(i3, FALSE, var_pars); // Plot variogram using gnuplot
}
- break;
- }
-
- if (type != 2) {
- plot_var(i3, FALSE, var_pars); // Plot variogram using gnuplot
- }
}
-void ordinary_kriging(struct int_par *xD, struct reg_par *reg, struct points *pnts, struct var_par *pars, struct output *out)
+void ordinary_kriging(struct int_par *xD, struct reg_par *reg,
+ struct points *pnts, struct var_par *pars,
+ struct output *out)
{
- // Local variables
- int i3 = xD->i3;
- double *vals = pnts->invals; // values to be used for interpolation
- struct write *report = &xD->report;
- struct write *crossvalid = &xD->crossvalid;
- struct parameters *var_par = &pars->fin;
+ // Local variables
+ int i3 = xD->i3;
+ double *vals = pnts->invals; // values to be used for interpolation
+ struct write *report = &xD->report;
+ struct write *crossvalid = &xD->crossvalid;
+ struct parameters *var_par = &pars->fin;
- int type = var_par->type;
- double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
- //max_dist = sqrt(0.5 * SQUARE(max_dist));
+ int type = var_par->type;
+ double max_dist =
+ type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
+ //max_dist = sqrt(0.5 * SQUARE(max_dist));
- double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
-
- unsigned int percents=50; // counter
- unsigned int row, col, dep; // indices of cells/voxels
- double rslt_OK; // interpolated value located on r0
+ double max_dist_vert =
+ type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
- pnts->max_dist = var_par->lag;
- struct ilist *list;
-
- double *r0; // xyz coordinates of cell/voxel centre
- mat_struct *GM;
- mat_struct *GM_sub; // submatrix of selected points
- mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
- mat_struct *g0; // diffences between known and unknown values = theor_var(dist)
- mat_struct *w0; // weights of values located on the input points
+ unsigned int percents = 50; // counter
+ unsigned int row, col, dep; // indices of cells/voxels
+ double rslt_OK; // interpolated value located on r0
- // Cell/voxel center coords (location of interpolated value)
- r0 = (double *) G_malloc(3 * sizeof(double));
-
- if (report->name) { // report file available:
- report->write2file = TRUE;
- report->fp = fopen(report->name, "a");
- time(&report->now);
- fprintf(report->fp, "Interpolating values started on %s\n\n", ctime(&report->now));
- }
+ pnts->max_dist = var_par->lag;
+ struct ilist *list;
- G_message(_("Interpolating unknown values..."));
- if (percents) {
- G_percent_reset();
- }
+ double *r0; // xyz coordinates of cell/voxel centre
+ mat_struct *GM;
+ mat_struct *GM_sub; // submatrix of selected points
+ mat_struct *GM_Inv; // inverted GM (GM_sub) matrix
+ mat_struct *g0; // diffences between known and unknown values = theor_var(dist)
+ mat_struct *w0; // weights of values located on the input points
- open_layer(xD, reg, out); // open 2D/3D raster
+ // Cell/voxel center coords (location of interpolated value)
+ r0 = (double *)G_malloc(3 * sizeof(double));
- if (var_par->const_val == 1) { // input values are constant:
- goto constant_voxel_centre;
- }
+ if (report->name) { // report file available:
+ report->write2file = TRUE;
+ report->fp = fopen(report->name, "a");
+ time(&report->now);
+ fprintf(report->fp, "Interpolating values started on %s\n\n",
+ ctime(&report->now));
+ }
- GM = set_up_G(pnts, var_par, &xD->report); // set up matrix of dissimilarities of input points
- var_par->GM = G_matrix_copy(GM); // copy matrix because of cross validation
+ G_message(_("Interpolating unknown values..."));
+ if (percents) {
+ G_percent_reset();
+ }
- // perform cross validation...
- if (crossvalid->name) { // ... if desired
- crossvalidation(xD, pnts, var_par);
- }
-
- constant_voxel_centre:
- for (dep=0; dep < reg->ndeps; dep++) {
- if (xD->i3 == TRUE) {
- if (percents) {
- G_percent(dep, reg->ndeps, 1);
- }
+ open_layer(xD, reg, out); // open 2D/3D raster
+
+ if (var_par->const_val == 1) { // input values are constant:
+ goto constant_voxel_centre;
}
- 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, GM, GM_Inv, g0, w0, rslt_OK)
- for (col=0; col < reg->ncols; col++) {
-
- if (var_par->const_val == 1) { // constant input values
- goto constant_voxel_val;
- }
- cell_centre(col, row, dep, xD, reg, r0, var_par); // coords of output point
+ GM = set_up_G(pnts, var_par, &xD->report); // set up matrix of dissimilarities of input points
+ var_par->GM = G_matrix_copy(GM); // copy matrix because of cross validation
- // add cell centre to the R-tree
- list = G_new_ilist(); // create list of overlapping rectangles
-
- if (i3 == TRUE) { // 3D kriging:
- list = find_NNs_within(3, r0, pnts, max_dist, max_dist_vert);
- }
- else { // 2D kriging:
- list = find_NNs_within(2, r0, pnts, max_dist, max_dist_vert);
- }
+ // perform cross validation...
+ if (crossvalid->name) { // ... if desired
+ crossvalidation(xD, pnts, var_par);
+ }
- if (list->n_values > 1) { // positive # of selected points:
- correct_indices(list, r0, 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, GM, GM_Inv, g0, w0, rslt_OK)
+ for (col = 0; col < reg->ncols; col++) {
- GM_sub = submatrix(list, GM, report); // make submatrix for selected points
- GM_Inv = G_matrix_inverse(GM_sub); // invert submatrix
- G_matrix_free(GM_sub);
+ if (var_par->const_val == 1) { // constant input values
+ goto constant_voxel_val;
+ }
- g0 = set_up_g0(xD, pnts, list, r0, var_par); // Diffs inputs - unknowns (incl. cond. 1))
- w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
-
- G_matrix_free(GM_Inv);
- G_matrix_free(g0);
+ cell_centre(col, row, dep, xD, reg, r0, var_par); // coords of output point
- rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
- G_matrix_free(w0);
- }
- else if (list->n_values == 1) {
- rslt_OK = vals[list->value[0] - 1]; // Estimated cell/voxel value rslt_OK = w x inputs
- }
- else if (list->n_values == 0) {
- if (report->name) { // report file available:
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("This point does not have neighbours in given radius..."));
- } // end else: error
+ // add cell centre to the R-tree
+ list = G_new_ilist(); // create list of overlapping rectangles
- G_free_ilist(list); // free list memory
+ if (i3 == TRUE) { // 3D kriging:
+ list =
+ find_NNs_within(3, r0, pnts, max_dist, max_dist_vert);
+ }
+ else { // 2D kriging:
+ list =
+ find_NNs_within(2, r0, pnts, max_dist, max_dist_vert);
+ }
- // Create output
- constant_voxel_val:
- if (var_par->const_val == 1) { // constant input values:
- rslt_OK = (double) *vals; // setup input as output
- }
+ if (list->n_values > 1) { // positive # of selected points:
+ correct_indices(list, r0, pnts, var_par);
- // write output to the (3D) raster layer
- if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
- if (report->name) { // report file available
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp); // close report file
- }
- G_fatal_error(_("Error writing result into output layer..."));
- }
- } // end col
- } // end row
- } // end dep
+ GM_sub = submatrix(list, GM, report); // make submatrix for selected points
+ GM_Inv = G_matrix_inverse(GM_sub); // invert submatrix
+ G_matrix_free(GM_sub);
- 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);
- }
+ g0 = set_up_g0(xD, pnts, list, r0, var_par); // Diffs inputs - unknowns (incl. cond. 1))
+ w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
- 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;
- }
+ G_matrix_free(GM_Inv);
+ G_matrix_free(g0);
+
+ rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
+ }
+ else if (list->n_values == 1) {
+ rslt_OK = vals[list->value[0] - 1]; // Estimated cell/voxel value rslt_OK = w x inputs
+ }
+ else if (list->n_values == 0) {
+ if (report->name) { // report file available:
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("This point does not have neighbours in given radius..."));
+ } // end else: error
+
+ G_free_ilist(list); // free list memory
+
+ // Create output
+ constant_voxel_val:
+ if (var_par->const_val == 1) { // constant input values:
+ rslt_OK = (double)*vals; // setup input as output
+ }
+
+ // write output to the (3D) raster layer
+ if (write2layer(xD, reg, out, col, row, dep, rslt_OK) == 0) {
+ if (report->name) { // report file available
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp); // close report file
+ }
+ G_fatal_error(_("Error writing result into output layer..."));
+ }
+ } // 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;
+ }
}
Modified: grass-addons/grass7/vector/v.kriging/getval.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/getval.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/getval.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -6,549 +6,583 @@
#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)
+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;
+ struct select *in_reg = &pnts->in_reg;
- int i, n, nrec, ctype;
- struct field_info *Fi;
+ int i, n, nrec, ctype;
+ struct field_info *Fi;
- dbCatValArray cvarr;
- dbDriver *Driver;
+ dbCatValArray cvarr;
+ dbDriver *Driver;
- int *index;
- double *values, *vals;
+ int *index;
+ double *values, *vals;
- int save;
- if (xD == NULL) {
- save = 0;
- }
- else {
- save = 1;
- }
-
- G_message(_("Reading values from attribute column <%s>..."), column);
+ int save;
- 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);
+ if (xD == NULL) {
+ save = 0;
}
- 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);
+ else {
+ save = 1;
}
- 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 */
+ G_message(_("Reading values from attribute column <%s>..."), column);
- /* 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);
+ 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);
}
- 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);
+ 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);
}
- G_fatal_error(_("Column must be numeric"));
- }
- db_close_database_shutdown_driver(Driver);
- G_free(Fi);
+ /* Note do not check if the column exists in the table because it may be expression */
- /* 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;
- }
+ /* TODO: only select values we need instead of all in column */
- values = (double *) G_malloc(n * sizeof(double));
- vals = &values[0];
+ /* 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);
+ }
- for (i = 0; i < n; i++) {
- /* TODO: Only for point data:
- * - store indexes of skipped points in read_points and compare them with database
- */
- if (index == NULL || (index != NULL && *index == i)) {
- if (ctype == DB_C_TYPE_INT)
- *vals = (double) cvarr.value[i].val.i;
- else if (ctype == DB_C_TYPE_DOUBLE) {
- *vals = cvarr.value[i].val.d;
- }
- vals++;
- if (index != NULL)
- index++;
+ 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"));
}
- }
- 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
+ db_close_database_shutdown_driver(Driver);
+ G_free(Fi);
- for (i = 0; i < n; i++) {
- G_matrix_set_element(A, i, 0, *x);
- G_matrix_set_element(A, i, 1, *y);
- G_matrix_set_element(gR, i, 0, *vals);
- x+=3;
- y+=3;
- z+=3;
- vals++;
+ /* 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;
}
- 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");
-
+ values = (double *)G_malloc(n * sizeof(double));
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++;
+ /* TODO: Only for point data:
+ * - store indexes of skipped points in read_points and compare them with database
+ */
+ if (index == NULL || (index != NULL && *index == i)) {
+ if (ctype == DB_C_TYPE_INT)
+ *vals = (double)cvarr.value[i].val.i;
+ else if (ctype == DB_C_TYPE_DOUBLE) {
+ *vals = cvarr.value[i].val.d;
+ }
+ vals++;
+ if (index != NULL)
+ index++;
+ }
}
- fclose(fp);
- //G_debug(0,"a=%f b=%f c=%f d=%f", T->vals[0], T->vals[1], T->vals[2], T->vals[3]);
- pnts->trend = T;
- }
- if (xD->phase == 0 && xD->report.name) {
- write2file_values(&xD->report, column);
- test_normality(n, values, &xD->report);
- }
+ if (detrend == 1) {
+ double *x, *y, *z;
- db_CatValArray_free(&cvarr); // free memory of the array of categories and values
-
- return values;
+ x = (double *)G_malloc(n * sizeof(double));
+ y = (double *)G_malloc(n * sizeof(double));
+ z = (double *)G_malloc(n * sizeof(double));
+
+ mat_struct *A, *gR, *T, *res;
+
+ x = &pnts->r[0];
+ y = &pnts->r[1];
+ z = &pnts->r[2];
+ vals = &values[0];
+ // Number of columns of design matrix A
+ A = G_matrix_init(n, 2, n); // initialise design matrix, normal grav: 7
+ gR = G_matrix_init(n, 1, n); // initialise vector of observations
+
+ for (i = 0; i < n; i++) {
+ G_matrix_set_element(A, i, 0, *x);
+ G_matrix_set_element(A, i, 1, *y);
+ G_matrix_set_element(gR, i, 0, *vals);
+ x += 3;
+ y += 3;
+ z += 3;
+ vals++;
+ }
+ T = LSM(A, gR); // Least Square Method
+ res = G_matrix_product(A, T);
+
+ FILE *fp;
+
+ x = &pnts->r[0];
+ y = &pnts->r[1];
+ z = &pnts->r[2];
+ fp = fopen("trend.txt", "w");
+
+ vals = &values[0];
+ double *resid;
+
+ resid = &res->vals[0];
+ for (i = 0; i < n; i++) {
+ *vals = *vals - *resid;
+ fprintf(fp, "%f %f %f %f\n", *x, *y, *z, *vals);
+ x += 3;
+ y += 3;
+ z += 3;
+ vals++;
+ resid++;
+ }
+ fclose(fp);
+ //G_debug(0,"a=%f b=%f c=%f d=%f", T->vals[0], T->vals[1], T->vals[2], T->vals[3]);
+ pnts->trend = T;
+ }
+
+ if (xD->phase == 0 && xD->report.name) {
+ write2file_values(&xD->report, column);
+ test_normality(n, values, &xD->report);
+ }
+
+ db_CatValArray_free(&cvarr); // free memory of the array of categories and values
+
+ return values;
}
/* get coordinates of input points */
-void read_points(struct Map_info *map, struct reg_par *reg, struct points *point, struct int_par *xD, const char *zcol, int field, struct write *report)
+void read_points(struct Map_info *map, struct reg_par *reg,
+ struct points *point, struct int_par *xD, const char *zcol,
+ int field, struct write *report)
{
- int type; // check elements of vector layer
- int i3 = xD->i3; // 2D / 3D interpolation
- int n; // # of all elements in vector layer
- int nskipped; // # of skipped elements (all except points)
- int n_in_reg = 0; // # of points within a region
- int out_reg = 0; // # of points out of region
- double *z_attr; // vertical coordinate from attribute value
-
- struct line_pnts *Points; // structures to hold line *Points (map)
+ int type; // check elements of vector layer
+ int i3 = xD->i3; // 2D / 3D interpolation
+ int n; // # of all elements in vector layer
+ int nskipped; // # of skipped elements (all except points)
+ int n_in_reg = 0; // # of points within a region
+ int out_reg = 0; // # of points out of region
+ double *z_attr; // vertical coordinate from attribute value
- // create spatial index R-tree
- point->Rtree_hz = RTreeCreateTree(-1, 0, 2); // create 2D spatial index
- if (i3 == TRUE) { // 3D:
- point->R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
- point->Rtree_vert = RTreeCreateTree(-1, 0, 1); // create 1D spatial index
- }
+ struct line_pnts *Points; // structures to hold line *Points (map)
- int sid=0; // id of the rectangle (will be increased by 1)
-
- double *rx, *ry, *rz; // pointers to coordinates
- int ind = 0; // index of the point
- int *indices; // indices of the points within the region
- int *index; // pointer to the vector of indices
+ // create spatial index R-tree
+ point->Rtree_hz = RTreeCreateTree(-1, 0, 2); // create 2D spatial index
+ if (i3 == TRUE) { // 3D:
+ point->R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
+ point->Rtree_vert = RTreeCreateTree(-1, 0, 1); // create 1D spatial index
+ }
- G_message(_("Reading coordinates..."));
+ int sid = 0; // id of the rectangle (will be increased by 1)
- Points = Vect_new_line_struct(); // Make point structure
- n = Vect_get_num_primitives(map, GV_POINT); // Number of points - topology required
+ double *rx, *ry, *rz; // pointers to coordinates
+ int ind = 0; // index of the point
+ int *indices; // indices of the points within the region
+ int *index; // pointer to the vector of indices
- point->r = (double *) G_malloc(n * 3 * sizeof(double)); // coords: x0 y0 z0 ... xn yn zn
- indices = (int *) G_malloc(n * sizeof(int)); // indices of pts within region
+ G_message(_("Reading coordinates..."));
- point->r_min = (double *) G_malloc(3 * sizeof(double)); // minimum
- point->r_max = (double *) G_malloc(3 * sizeof(double)); // maximum
+ Points = Vect_new_line_struct(); // Make point structure
+ n = Vect_get_num_primitives(map, GV_POINT); // Number of points - topology required
- rx = &point->r[0]; // pointer to x coordinate
- ry = &point->r[1]; // pointer to y coordinate
- rz = &point->r[2]; // pointer to z coordinate
- index = &indices[0]; // pointer to indices
-
- // Get 3rd coordinate of 2D points from attribute column -> 3D interpolation
- if (xD->v3 == FALSE && zcol != NULL) {
- z_attr = get_col_values(map, NULL, point, field, zcol, FALSE);
- }
+ point->r = (double *)G_malloc(n * 3 * sizeof(double)); // coords: x0 y0 z0 ... xn yn zn
+ indices = (int *)G_malloc(n * sizeof(int)); // indices of pts within region
- 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;
- }
+ point->r_min = (double *)G_malloc(3 * sizeof(double)); // minimum
+ point->r_max = (double *)G_malloc(3 * sizeof(double)); // maximum
- // skip everything except points
- if (type != GV_POINT) {
- nskipped++;
- continue;
- }
+ rx = &point->r[0]; // pointer to x coordinate
+ ry = &point->r[1]; // pointer to y coordinate
+ rz = &point->r[2]; // pointer to z coordinate
+ index = &indices[0]; // pointer to indices
- 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);
+ // 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);
}
- // 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
+ nskipped = 0; // # of skipped elements (lines, areas)
- insert_rectangle(2, sid, point); // R-tree node 2D (2D and 3D kriging)
-
- if (i3 == TRUE) { // 3D interpolation:
- if (reg->bot <= Points->z[0] && Points->z[0] <= reg->top) {
- if (zcol == NULL) { // 3D layer:
- *rz = Points->z[0]; // set up z coordinate
+ // 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;
+ }
- if (*rz != Points->z[0]) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("Error reading input coordinates z..."));
- }
- }
- else { // 2D layer with z-column:
- *rz = *z_attr; // set up x coordinate
- } // end else (2,5D layer)
+ // skip everything except points
+ if (type != GV_POINT) {
+ nskipped++;
+ continue;
+ }
- insert_rectangle(3, sid, point); // R-tree node 3D (just 3D kriging)
- insert_rectangle(1, sid, point); // R-tree node 1D (just 3D kriging)
- } // end if rb, rt
+ 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);
+ }
- else { // point is out of region:
- goto out_of_region;
- }
- } // end if 3D interpolation
-
- else { // 2D interpolation:
- *rz = 0.; // set up z coordinate
- }
+ // 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
- // Find extends
- if (n_in_reg == 0) {
- triple(*rx, *ry, *rz, point->r_min);
- triple(*rx, *ry, *rz, point->r_max);
- }
+ insert_rectangle(2, sid, point); // R-tree node 2D (2D and 3D kriging)
- else {
- triple(MIN(*point->r_min, *rx), MIN(*(point->r_min+1), *ry), MIN(*(point->r_min+2), *rz), point->r_min);
- triple(MAX(*point->r_max, *rx), MAX(*(point->r_max+1), *ry), MAX(*(point->r_max+2), *rz), point->r_max);
- } // end else (extent test)
-
- n_in_reg++; // # of points in region
- *index = ind; // store index of point within the region
- index++; // go to next index
+ if (i3 == TRUE) { // 3D interpolation:
+ if (reg->bot <= Points->z[0] && Points->z[0] <= reg->top) {
+ if (zcol == NULL) { // 3D layer:
+ *rz = Points->z[0]; // set up z coordinate
- rx += 3; // go to new x coordinate
- ry += 3; // go to new y coordinate
- rz += 3; // go to new z coordinate
+ if (*rz != Points->z[0]) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error reading input coordinates z..."));
+ }
+ }
+ else { // 2D layer with z-column:
+ *rz = *z_attr; // set up x coordinate
+ } // end else (2,5D layer)
- sid++; // go to new ID for spatial indexing
-
- goto finish;
- } // end in region
+ insert_rectangle(3, sid, point); // R-tree node 3D (just 3D kriging)
+ insert_rectangle(1, sid, point); // R-tree node 1D (just 3D kriging)
+ } // end if rb, rt
- out_of_region:
- out_reg++;
- continue;
+ else { // point is out of region:
+ goto out_of_region;
+ }
+ } // end if 3D interpolation
- finish:
- if (zcol != NULL) {
- z_attr++;
- }
- ind++;
- } // end while (TRUE)
-
- point->in_reg.total = n;
- point->in_reg.n = n_in_reg;
- point->in_reg.out = out_reg;
- point->n = n_in_reg;
+ else { // 2D interpolation:
+ *rz = 0.; // set up z coordinate
+ }
- point->in_reg.indices = (int *) G_malloc(n_in_reg * sizeof(int));
- memcpy(point->in_reg.indices, indices, n_in_reg * sizeof(int));
- G_free(indices);
+ // Find extends
+ if (n_in_reg == 0) {
+ triple(*rx, *ry, *rz, point->r_min);
+ triple(*rx, *ry, *rz, point->r_max);
+ }
- if (nskipped > 0) {
- G_warning(_("%d features skipped, only points are accepted"), nskipped);
- }
+ else {
+ triple(MIN(*point->r_min, *rx), MIN(*(point->r_min + 1), *ry),
+ MIN(*(point->r_min + 2), *rz), point->r_min);
+ triple(MAX(*point->r_max, *rx), MAX(*(point->r_max + 1), *ry),
+ MAX(*(point->r_max + 2), *rz), point->r_max);
+ } // end else (extent test)
- Vect_destroy_line_struct(Points);
+ n_in_reg++; // # of points in region
+ *index = ind; // store index of point within the region
+ index++; // go to next index
- if (n_in_reg == 0) {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
+ rx += 3; // go to new x coordinate
+ ry += 3; // go to new y coordinate
+ rz += 3; // go to new z coordinate
+
+ sid++; // go to new ID for spatial indexing
+
+ goto finish;
+ } // end in region
+
+ out_of_region:
+ out_reg++;
+ continue;
+
+ finish:
+ if (zcol != NULL) {
+ z_attr++;
+ }
+ ind++;
+ } // end while (TRUE)
+
+ point->in_reg.total = n;
+ point->in_reg.n = n_in_reg;
+ point->in_reg.out = out_reg;
+ point->n = n_in_reg;
+
+ point->in_reg.indices = (int *)G_malloc(n_in_reg * sizeof(int));
+ memcpy(point->in_reg.indices, indices, n_in_reg * sizeof(int));
+ G_free(indices);
+
+ if (nskipped > 0) {
+ G_warning(_("%d features skipped, only points are accepted"),
+ nskipped);
}
- G_fatal_error(_("Unable to read coordinates of point layer (all points are out of region)"));
- }
- if (out_reg > 0) {
- G_message(_("Unused points: %d (out of region)"), out_reg);
- }
+ Vect_destroy_line_struct(Points);
- if (xD->phase == 0) { // initial phase:
- write2file_vector(xD, point); // describe properties
- }
+ if (n_in_reg == 0) {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Unable to read coordinates of point layer (all points are out of region)"));
+ }
+
+ if (out_reg > 0) {
+ G_message(_("Unused points: %d (out of region)"), out_reg);
+ }
+
+ if (xD->phase == 0) { // initial phase:
+ write2file_vector(xD, point); // describe properties
+ }
}
void get_region_pars(struct int_par *xD, struct reg_par *reg)
{
- /* Get extent and resolution of output raster 2D/3D */
- if (xD->i3 == FALSE) { /* 2D region */
- /* Get database region parameters */
- Rast_get_window(®->reg_2d); /* stores the current default window in region */
- reg->west = reg->reg_2d.west;
- reg->east = reg->reg_2d.east;
- reg->north = reg->reg_2d.north;
- reg->south = reg->reg_2d.south;
- reg->ew_res = reg->reg_2d.ew_res;
- reg->ns_res = reg->reg_2d.ns_res;
- reg->nrows = reg->reg_2d.rows;
- reg->ncols = reg->reg_2d.cols;
- reg->ndeps = 1;
- }
- if (xD->i3 == TRUE) { /* 3D region */
- Rast3d_init_defaults();
- /* Get database region parameters */
- Rast3d_get_window(®->reg_3d); /* stores the current default window in region */
- if (reg->reg_3d.bottom == 0 && (reg->reg_3d.tb_res == 0 || reg->reg_3d.depths == 0)) {
- if (xD->report.write2file == TRUE) { // close report file
- fprintf(xD->report.fp, "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
- }
- G_fatal_error("To process 3D interpolation, please set 3D region settings.");
+ /* 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;
}
- 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 */
+ if (xD->i3 == TRUE) { /* 3D region */
+ Rast3d_init_defaults();
+ /* Get database region parameters */
+ Rast3d_get_window(®->reg_3d); /* stores the current default window in region */
+ if (reg->reg_3d.bottom == 0 &&
+ (reg->reg_3d.tb_res == 0 || reg->reg_3d.depths == 0)) {
+ if (xD->report.write2file == TRUE) { // close report file
+ fprintf(xD->report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error
+ ("To process 3D interpolation, please set 3D region settings.");
+ }
+ reg->west = reg->reg_3d.west;
+ reg->east = reg->reg_3d.east;
+ reg->north = reg->reg_3d.north;
+ reg->south = reg->reg_3d.south;
+ reg->bot = reg->reg_3d.bottom; /* output 3D raster */
+ reg->top = reg->reg_3d.top; /* output 3D raster */
+ reg->ew_res = reg->reg_3d.ew_res;
+ reg->ns_res = reg->reg_3d.ns_res;
+ reg->bt_res = reg->reg_3d.tb_res;
+ reg->nrows = reg->reg_3d.rows;
+ reg->ncols = reg->reg_3d.cols;
+ reg->ndeps = xD->i3 == TRUE ? reg->reg_3d.depths : 1; /* 2D/3D interpolation */
+ }
+ reg->nrcd = reg->nrows * reg->ncols * reg->ndeps; /* number of cells */
}
// read values from temporary files
-void read_tmp_vals(const char *file_name, struct parameters *var_par, struct int_par *xD)
+void read_tmp_vals(const char *file_name, struct parameters *var_par,
+ struct int_par *xD)
{
- FILE *fp;
+ FILE *fp;
- int j, nLag_vert;
- double lag_vert, max_dist_vert;
- double *v_elm, *g_elm;
- double sill_hz, sill_vert;
+ int j, nLag_vert;
+ double lag_vert, max_dist_vert;
+ double *v_elm, *g_elm;
+ double sill_hz, sill_vert;
-
- fp = fopen(file_name, "r");
- if (fp == NULL) {
- G_fatal_error(_("Temporary file is missing, please repeat an initial phase..."));
- }
-
- else { // file exists:
- int i, type;
- int nLag;
- double lag, max_dist, td_hz, sill;
- double *h_elm, *gamma;
- int file, file_length;
- for (i=0; i<2; i++) {
- if (fscanf(fp, "%d", &file_length) == 0) { // filename length
- G_fatal_error(_("Nothing to scan..."));
- }
- if (file_length > 3) {
- if (fscanf(fp, "%d", &file) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- if (file == 9) { // filetype code (9 - report, 8 - crossval)
- xD->report.name = (char *) G_malloc(file_length * sizeof(char));
- if (fscanf(fp, "%s", xD->report.name) == 0) { // filename
- G_fatal_error(_("Nothing to scan..."));
- }
- continue;
- }
- else if (file == 8) {
- xD->crossvalid.name = (char *) G_malloc(file_length * sizeof(char));
- if (fscanf(fp, "%s", xD->crossvalid.name) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- continue;
- }
- }
- else {
- type = file_length;
- goto no_file;
- }
- } // todo: test without report and crossval
+ fp = fopen(file_name, "r");
+ if (fp == NULL) {
+ G_fatal_error(_("Temporary file is missing, please repeat an initial phase..."));
+ }
- if (fscanf(fp, "%d", &type) == 0) { // read type
- G_fatal_error(_("Nothing to scan..."));
- }
-
- no_file:
-
- switch (type) {
- case 2: // bivariate variogram
- var_par->type = type;
- if (fscanf(fp, "%d %lf %lf %d %lf %lf %lf", &nLag_vert, &lag_vert, &max_dist_vert, &nLag, &lag, &max_dist, &td_hz) < 7) {
- G_fatal_error(_("Nothing to scan..."));
- }
+ else { // file exists:
+ int i, type;
+ int nLag;
+ double lag, max_dist, td_hz, sill;
+ double *h_elm, *gamma;
+ int file, file_length;
- var_par->horizontal.h = (double *) G_malloc(nLag * sizeof(double));
- h_elm = &var_par->horizontal.h[0];
- for (i=0; i<nLag; i++) {
- if (fscanf(fp, "%lf", h_elm) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- h_elm++;
- }
+ for (i = 0; i < 2; i++) {
+ if (fscanf(fp, "%d", &file_length) == 0) { // filename length
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ if (file_length > 3) {
+ if (fscanf(fp, "%d", &file) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ if (file == 9) { // filetype code (9 - report, 8 - crossval)
+ xD->report.name =
+ (char *)G_malloc(file_length * sizeof(char));
+ if (fscanf(fp, "%s", xD->report.name) == 0) { // filename
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ continue;
+ }
+ else if (file == 8) {
+ xD->crossvalid.name =
+ (char *)G_malloc(file_length * sizeof(char));
+ if (fscanf(fp, "%s", xD->crossvalid.name) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ continue;
+ }
+ }
+ else {
+ type = file_length;
+ goto no_file;
+ }
+ } // todo: test without report and crossval
- var_par->vertical.h = (double *) G_malloc(nLag_vert * sizeof(double));
- v_elm = &var_par->vertical.h[0];
- for (i=0; i<nLag_vert; i++) {
- if (fscanf(fp, "%lf", v_elm) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- v_elm++;
- }
+ if (fscanf(fp, "%d", &type) == 0) { // read type
+ G_fatal_error(_("Nothing to scan..."));
+ }
- 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++;
- }
- }
+ no_file:
- if (fscanf(fp, "%lf %lf", &sill_hz, &sill_vert) < 2) {
- G_fatal_error(_("Nothing to scan..."));
- }
+ switch (type) {
+ case 2: // bivariate variogram
+ var_par->type = type;
+ if (fscanf
+ (fp, "%d %lf %lf %d %lf %lf %lf", &nLag_vert, &lag_vert,
+ &max_dist_vert, &nLag, &lag, &max_dist, &td_hz) < 7) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
- var_par->vertical.nLag = nLag_vert;
- var_par->vertical.lag = lag_vert;
- var_par->vertical.max_dist = max_dist_vert;
- var_par->vertical.sill = sill_vert;
+ var_par->horizontal.h = (double *)G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->horizontal.h[0];
+ for (i = 0; i < nLag; i++) {
+ if (fscanf(fp, "%lf", h_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ h_elm++;
+ }
- var_par->horizontal.nLag = nLag;
- var_par->horizontal.lag = lag;
- var_par->horizontal.max_dist = max_dist;
- var_par->horizontal.sill = sill_hz;
- break;
+ var_par->vertical.h =
+ (double *)G_malloc(nLag_vert * sizeof(double));
+ v_elm = &var_par->vertical.h[0];
+ for (i = 0; i < nLag_vert; i++) {
+ if (fscanf(fp, "%lf", v_elm) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ v_elm++;
+ }
- default:
- if (type == 3) { // anisotropic variogram:
- if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0) { // anisotropic ratio
- G_fatal_error(_("Nothing to scan..."));
- }
- }
+ 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, "%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..."));
- }
- }
+ if (fscanf(fp, "%lf %lf", &sill_hz, &sill_vert) < 2) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
- var_par->h = (double *) G_malloc(nLag * sizeof(double));
- h_elm = &var_par->h[0];
+ var_par->vertical.nLag = nLag_vert;
+ var_par->vertical.lag = lag_vert;
+ var_par->vertical.max_dist = max_dist_vert;
+ var_par->vertical.sill = sill_vert;
- var_par->gamma = G_matrix_init(nLag, 1, nLag);
- gamma = &var_par->gamma->vals[0];
-
- for (i=0; i<nLag; i++) {
- if (fscanf(fp, "%lf %lf", h_elm, gamma) < 2) {
- G_fatal_error(_("Nothing to scan..."));
- }
- h_elm++;
- gamma++;
- }
+ var_par->horizontal.nLag = nLag;
+ var_par->horizontal.lag = lag;
+ var_par->horizontal.max_dist = max_dist;
+ var_par->horizontal.sill = sill_hz;
+ break;
- if (fscanf(fp, "%lf", &sill) == 0) {
- G_fatal_error(_("Nothing to scan..."));
- }
- var_par->sill = sill;
- break;
+ default:
+ if (type == 3) { // anisotropic variogram:
+ if (fscanf(fp, "%lf", &xD->aniso_ratio) == 0) { // anisotropic ratio
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ }
+
+ if (fscanf(fp, "%d %lf %lf", &nLag, &lag, &max_dist) < 3) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+
+ if (type != 1) {
+ if (fscanf(fp, "%lf", &td_hz) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ }
+
+ var_par->h = (double *)G_malloc(nLag * sizeof(double));
+ h_elm = &var_par->h[0];
+
+ var_par->gamma = G_matrix_init(nLag, 1, nLag);
+ gamma = &var_par->gamma->vals[0];
+
+ for (i = 0; i < nLag; i++) {
+ if (fscanf(fp, "%lf %lf", h_elm, gamma) < 2) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ h_elm++;
+ gamma++;
+ }
+
+ if (fscanf(fp, "%lf", &sill) == 0) {
+ G_fatal_error(_("Nothing to scan..."));
+ }
+ var_par->sill = sill;
+ break;
+ }
+
+ if (type != 2) {
+ var_par->type = type;
+ var_par->nLag = nLag;
+ var_par->lag = lag;
+ var_par->max_dist = max_dist;
+ var_par->td = td_hz;
+ }
}
+ fclose(fp);
- if (type != 2) {
- var_par->type = type;
- var_par->nLag = nLag;
- var_par->lag = lag;
- var_par->max_dist = max_dist;
- var_par->td = td_hz;
+ if (xD->phase == 2) {
+ remove(file_name);
}
- }
- fclose(fp);
-
- if (xD->phase == 2) {
- remove(file_name);
- }
}
Modified: grass-addons/grass7/vector/v.kriging/main.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/main.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/main.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -24,579 +24,609 @@
int main(int argc, char *argv[])
{
- // Vector layer and module
- struct Map_info map; // Input vector map
- struct GModule *module; // Module
+ // 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.)
-
- // Geostatistical parameters
- struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
- struct var_par var_pars; // Variogram (experimental and theoretical)
+ struct reg_par reg; // Region parameters
+ struct points pnts; // Input points (coordinates, extent, values, etc.)
- // Outputs
- struct output out; // Output layer properties
-
- // Settings
- int field;
- struct opts opt;
- struct flgs flg;
+ // Geostatistical parameters
+ struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
+ struct var_par var_pars; // Variogram (experimental and theoretical)
- /* ------- 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).");
+ // Outputs
+ struct output out; // Output layer properties
- // Setting options
- opt.input = G_define_standard_option(G_OPT_V_INPUT); // Vector input layer
- opt.input->label = _("Name of input vector points map");
+ // Settings
+ int field;
+ struct opts opt;
+ struct flgs flg;
- 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");
+ /* ------- 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).");
- opt.field = G_define_standard_option(G_OPT_V_FIELD);
+ // Setting options
+ opt.input = G_define_standard_option(G_OPT_V_INPUT); // Vector input layer
+ opt.input->label = _("Name of input vector points map");
- opt.phase = G_define_option();
- opt.phase->key = "phase";
- opt.phase->options = "initial, middle, final";
- opt.phase->description = _("Phase of interpolation. In the initial phase, there is empirical variogram computed. In the middle phase, function of theoretical variogram is chosen by the user and its coefficients are estimated empirically. In the final phase, unknown values are interpolated using theoretical variogram from previous phase.");
- opt.phase->required = YES;
+ 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.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.field = G_define_standard_option(G_OPT_V_FIELD);
- opt.report = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
- opt.report->key = "report";
- opt.report->description = _("File to write the report");
- opt.report->required = NO;
- opt.report->guisection = _("Initial");
+ opt.phase = G_define_option();
+ opt.phase->key = "phase";
+ opt.phase->options = "initial, middle, final";
+ opt.phase->description =
+ _("Phase of interpolation. In the initial phase, there is empirical variogram computed. In the middle phase, function of theoretical variogram is chosen by the user and its coefficients are estimated empirically. In the final phase, unknown values are interpolated using theoretical variogram from previous phase.");
+ opt.phase->required = YES;
- opt.crossvalid = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
- opt.crossvalid->key = "crossvalid";
- opt.crossvalid->description = _("File to write the results of cross validation");
- opt.crossvalid->required = NO;
- opt.crossvalid->guisection = _("Final");
+ 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");
- flg.bivariate = G_define_flag();
- flg.bivariate->key = 'b';
- flg.bivariate->description = _("Compute bivariate variogram (3D interpolation only)");
- flg.bivariate->guisection = _("Middle");
+ opt.report = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.report->key = "report";
+ opt.report->description = _("File to write the report");
+ opt.report->required = NO;
+ opt.report->guisection = _("Initial");
- flg.univariate = G_define_flag();
- flg.univariate->key = 'u';
- flg.univariate->description = _("Compute univariate variogram (3D interpolation only)");
- flg.univariate->guisection = _("Middle");
+ opt.crossvalid = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.crossvalid->key = "crossvalid";
+ opt.crossvalid->description =
+ _("File to write the results of cross validation");
+ opt.crossvalid->required = NO;
+ opt.crossvalid->guisection = _("Final");
- 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 = _("Horizontal variogram function");
- opt.function_var_hz->guisection = _("Middle");
+ flg.bivariate = G_define_flag();
+ flg.bivariate->key = 'b';
+ flg.bivariate->description =
+ _("Compute bivariate variogram (3D interpolation only)");
+ flg.bivariate->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 = _("Vertical variogram function");
- opt.function_var_vert->guisection = _("Middle");
+ flg.univariate = G_define_flag();
+ flg.univariate->key = 'u';
+ flg.univariate->description =
+ _("Compute univariate variogram (3D interpolation only)");
+ flg.univariate->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 = _("Final variogram function (anisotropic or horizontal component of bivariate variogram)");
- opt.function_var_final->guisection = _("Final");
+ 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 = _("Horizontal variogram function");
+ opt.function_var_hz->guisection = _("Middle");
- opt.function_var_final_vert = G_define_option(); // Variogram type
- opt.function_var_final_vert->key = "final_vert_function";
- opt.function_var_final_vert->options = "linear, exponential, spherical, gaussian, bivariate";
- opt.function_var_final_vert->description = _("Final variogram function (vertical component of bivariate variogram)");
- opt.function_var_final_vert->guisection = _("Final");
+ 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 = _("Vertical variogram function");
+ opt.function_var_vert->guisection = _("Middle");
- flg.detrend = G_define_flag();
- flg.detrend->key = 't';
- flg.detrend->description = _("Eliminate trend if variogram is parabolic");
- flg.detrend->guisection = _("Initial");
+ 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 =
+ _("Final variogram function (anisotropic or horizontal component of bivariate variogram)");
+ opt.function_var_final->guisection = _("Final");
- opt.form_file = G_define_option(); // Variogram plot - various output formats
- opt.form_file->key = "fileformat";
- opt.form_file->options = "cdr,dxf,eps,tex,pdf,png,svg";
- opt.form_file->description = _("File format to save variogram plot (empty: preview in Gnuplot terminal)");
- opt.form_file->guisection = _("Middle");
+ opt.function_var_final_vert = G_define_option(); // Variogram type
+ opt.function_var_final_vert->key = "final_vert_function";
+ opt.function_var_final_vert->options =
+ "linear, exponential, spherical, gaussian, bivariate";
+ opt.function_var_final_vert->description =
+ _("Final variogram function (vertical component of bivariate variogram)");
+ opt.function_var_final_vert->guisection = _("Final");
- opt.intpl = G_define_standard_option(G_OPT_DB_COLUMN); // Input values for interpolation
- opt.intpl->key = "icolumn";
- opt.intpl->description =
- _("Attribute column containing input values for interpolation");
- opt.intpl->required = YES;
+ flg.detrend = G_define_flag();
+ flg.detrend->key = 't';
+ flg.detrend->description = _("Eliminate trend if variogram is parabolic");
+ flg.detrend->guisection = _("Initial");
- opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN); // Column with z coord (2D points)
- opt.zcol->key = "zcolumn";
- opt.zcol->description =
- _("Attribute column containing z coordinates (only for 3D interpolation based on 2D point layer)");
- opt.zcol->required = NO;
- opt.zcol->guisection = _("3D");
+ 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 =
+ _("File format to save variogram plot (empty: preview in Gnuplot terminal)");
+ opt.form_file->guisection = _("Middle");
- 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 = "45.0";
- opt.var_dir_hz->description =
- _("Azimuth of variogram computing (isotrophic)");
- opt.var_dir_hz->guisection = _("Initial");
+ opt.intpl = G_define_standard_option(G_OPT_DB_COLUMN); // Input values for interpolation
+ opt.intpl->key = "icolumn";
+ opt.intpl->description =
+ _("Attribute column containing input values for interpolation");
+ opt.intpl->required = YES;
- 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 variogram computing (isotrophic)");
- opt.var_dir_vert->guisection = _("Initial");
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN); // Column with z coord (2D points)
+ opt.zcol->key = "zcolumn";
+ opt.zcol->description =
+ _("Attribute column containing z coordinates (only for 3D interpolation based on 2D point layer)");
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("3D");
- opt.nL = G_define_option();
- opt.nL->key = "lpieces";
- opt.nL->type = TYPE_INTEGER;
- opt.nL->required = NO;
- opt.nL->description = _("Number of horizontal lags");
- opt.nL->guisection = _("Initial");
+ 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 = "45.0";
+ opt.var_dir_hz->description =
+ _("Azimuth of variogram computing (isotrophic)");
+ opt.var_dir_hz->guisection = _("Initial");
- opt.nZ = G_define_option();
- opt.nZ->key = "vpieces";
- opt.nZ->type = TYPE_INTEGER;
- opt.nZ->required = NO;
- opt.nZ->description =
- _("Number of vertical lags (only for 3D variogram)");
- opt.nZ->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 variogram computing (isotrophic)");
+ opt.var_dir_vert->guisection = _("Initial");
- opt.td_hz = G_define_option();
- opt.td_hz->key = "td";
- opt.td_hz->type = TYPE_DOUBLE;
- opt.td_hz->answer = "45.0";
- opt.td_hz->description = _("Angle of variogram processing");
- opt.td_hz->required = NO;
- opt.td_hz->guisection = _("Initial");
+ opt.nL = G_define_option();
+ opt.nL->key = "lpieces";
+ opt.nL->type = TYPE_INTEGER;
+ opt.nL->required = NO;
+ opt.nL->description = _("Number of horizontal lags");
+ opt.nL->guisection = _("Initial");
- opt.nugget_hz = G_define_option();
- opt.nugget_hz->key = "hz_nugget";
- opt.nugget_hz->type = TYPE_DOUBLE;
- opt.nugget_hz->answer = "0.0";
- opt.nugget_hz->description = _("Nugget effect of horizontal variogram");
- opt.nugget_hz->required = NO;
- opt.nugget_hz->guisection = _("Middle");
+ opt.nZ = G_define_option();
+ opt.nZ->key = "vpieces";
+ opt.nZ->type = TYPE_INTEGER;
+ opt.nZ->required = NO;
+ opt.nZ->description =
+ _("Number of vertical lags (only for 3D variogram)");
+ opt.nZ->guisection = _("Initial");
- opt.nugget_vert = G_define_option();
- opt.nugget_vert->key = "vert_nugget";
- opt.nugget_vert->type = TYPE_DOUBLE;
- opt.nugget_vert->answer = "0.0";
- opt.nugget_vert->description = _("Nugget effect of vertical variogram");
- opt.nugget_vert->required = NO;
- opt.nugget_vert->guisection = _("Middle");
+ opt.td_hz = G_define_option();
+ opt.td_hz->key = "td";
+ opt.td_hz->type = TYPE_DOUBLE;
+ opt.td_hz->answer = "45.0";
+ opt.td_hz->description = _("Angle of variogram processing");
+ opt.td_hz->required = NO;
+ opt.td_hz->guisection = _("Initial");
- opt.nugget_final = G_define_option();
- opt.nugget_final->key = "final_nugget";
- opt.nugget_final->type = TYPE_DOUBLE;
- opt.nugget_final->answer = "0.0";
- opt.nugget_final->description = _("Nugget effect of anisotropic variogram (or horizontal component of bivariate variogram)");
- opt.nugget_final->required = NO;
- opt.nugget_final->guisection = _("Final");
+ opt.nugget_hz = G_define_option();
+ opt.nugget_hz->key = "hz_nugget";
+ opt.nugget_hz->type = TYPE_DOUBLE;
+ opt.nugget_hz->answer = "0.0";
+ opt.nugget_hz->description = _("Nugget effect of horizontal variogram");
+ opt.nugget_hz->required = NO;
+ opt.nugget_hz->guisection = _("Middle");
- opt.nugget_final_vert = G_define_option();
- opt.nugget_final_vert->key = "final_vert_nugget";
- opt.nugget_final_vert->type = TYPE_DOUBLE;
- opt.nugget_final_vert->answer = "0.0";
- opt.nugget_final_vert->description = _("For bivariate variogram only: nuget effect of vertical component");
- opt.nugget_final_vert->required = NO;
- opt.nugget_final_vert->guisection = _("Final");
+ opt.nugget_vert = G_define_option();
+ opt.nugget_vert->key = "vert_nugget";
+ opt.nugget_vert->type = TYPE_DOUBLE;
+ opt.nugget_vert->answer = "0.0";
+ opt.nugget_vert->description = _("Nugget effect of vertical variogram");
+ opt.nugget_vert->required = NO;
+ opt.nugget_vert->guisection = _("Middle");
- opt.sill_hz = G_define_option();
- opt.sill_hz->key = "hz_sill";
- opt.sill_hz->type = TYPE_DOUBLE;
- opt.sill_hz->description = _("Sill of horizontal variogram");
- opt.sill_hz->required = NO;
- opt.sill_hz->guisection = _("Middle");
+ opt.nugget_final = G_define_option();
+ opt.nugget_final->key = "final_nugget";
+ opt.nugget_final->type = TYPE_DOUBLE;
+ opt.nugget_final->answer = "0.0";
+ opt.nugget_final->description =
+ _("Nugget effect of anisotropic variogram (or horizontal component of bivariate variogram)");
+ opt.nugget_final->required = NO;
+ opt.nugget_final->guisection = _("Final");
- opt.sill_vert = G_define_option();
- opt.sill_vert->key = "vert_sill";
- opt.sill_vert->type = TYPE_DOUBLE;
- opt.sill_vert->description = _("Sill of vertical variogram");
- opt.sill_vert->required = NO;
- opt.sill_vert->guisection = _("Middle");
+ opt.nugget_final_vert = G_define_option();
+ opt.nugget_final_vert->key = "final_vert_nugget";
+ opt.nugget_final_vert->type = TYPE_DOUBLE;
+ opt.nugget_final_vert->answer = "0.0";
+ opt.nugget_final_vert->description =
+ _("For bivariate variogram only: nuget effect of vertical component");
+ opt.nugget_final_vert->required = NO;
+ opt.nugget_final_vert->guisection = _("Final");
- opt.sill_final = G_define_option();
- opt.sill_final->key = "final_sill";
- opt.sill_final->type = TYPE_DOUBLE;
- opt.sill_final->description = _("Sill of anisotropic variogram (or horizontal component of bivariate variogram)");
- opt.sill_final->required = NO;
- opt.sill_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 = _("Sill of horizontal variogram");
+ opt.sill_hz->required = NO;
+ opt.sill_hz->guisection = _("Middle");
- opt.sill_final_vert = G_define_option();
- opt.sill_final_vert->key = "final_vert_sill";
- opt.sill_final_vert->type = TYPE_DOUBLE;
- opt.sill_final_vert->description = _("For bivariate variogram only: sill of vertical component");
- opt.sill_final_vert->required = NO;
- opt.sill_final_vert->guisection = _("Final");
+ opt.sill_vert = G_define_option();
+ opt.sill_vert->key = "vert_sill";
+ opt.sill_vert->type = TYPE_DOUBLE;
+ opt.sill_vert->description = _("Sill of vertical variogram");
+ opt.sill_vert->required = NO;
+ opt.sill_vert->guisection = _("Middle");
- opt.range_hz = G_define_option();
- opt.range_hz->key = "hz_range";
- opt.range_hz->type = TYPE_DOUBLE;
- opt.range_hz->description = _("Range of horizontal variogram");
- opt.range_hz->required = NO;
- opt.range_hz->guisection = _("Middle");
+ opt.sill_final = G_define_option();
+ opt.sill_final->key = "final_sill";
+ opt.sill_final->type = TYPE_DOUBLE;
+ opt.sill_final->description =
+ _("Sill of anisotropic variogram (or horizontal component of bivariate variogram)");
+ opt.sill_final->required = NO;
+ opt.sill_final->guisection = _("Final");
- opt.range_vert = G_define_option();
- opt.range_vert->key = "vert_range";
- opt.range_vert->type = TYPE_DOUBLE;
- opt.range_vert->description = _("Range of vertical variogram");
- opt.range_vert->required = NO;
- opt.range_vert->guisection = _("Middle");
+ opt.sill_final_vert = G_define_option();
+ opt.sill_final_vert->key = "final_vert_sill";
+ opt.sill_final_vert->type = TYPE_DOUBLE;
+ opt.sill_final_vert->description =
+ _("For bivariate variogram only: sill of vertical component");
+ opt.sill_final_vert->required = NO;
+ opt.sill_final_vert->guisection = _("Final");
- opt.range_final = G_define_option();
- opt.range_final->key = "final_range";
- opt.range_final->type = TYPE_DOUBLE;
- opt.range_final->description = _("Range of anisotropic variogram (or horizontal component of bivariate variogram)");
- opt.range_final->required = NO;
- opt.range_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 = _("Range of horizontal variogram");
+ opt.range_hz->required = NO;
+ opt.range_hz->guisection = _("Middle");
- opt.range_final_vert = G_define_option();
- opt.range_final_vert->key = "final_vert_range";
- opt.range_final_vert->type = TYPE_DOUBLE;
- opt.range_final_vert->description = _("Range of final variogram: one value for anisotropic, two values for bivariate (hz and vert component)");
- opt.range_final_vert->required = NO;
- opt.range_final_vert->guisection = _("Final");
- /* --------------------------------------------------------- */
+ opt.range_vert = G_define_option();
+ opt.range_vert->key = "vert_range";
+ opt.range_vert->type = TYPE_DOUBLE;
+ opt.range_vert->description = _("Range of vertical variogram");
+ opt.range_vert->required = NO;
+ opt.range_vert->guisection = _("Middle");
- G_gisinit(argv[0]);
+ opt.range_final = G_define_option();
+ opt.range_final->key = "final_range";
+ opt.range_final->type = TYPE_DOUBLE;
+ opt.range_final->description =
+ _("Range of anisotropic variogram (or horizontal component of bivariate variogram)");
+ opt.range_final->required = NO;
+ opt.range_final->guisection = _("Final");
- if (G_parser(argc, argv)) {
- exit(EXIT_FAILURE);
- }
+ opt.range_final_vert = G_define_option();
+ opt.range_final_vert->key = "final_vert_range";
+ opt.range_final_vert->type = TYPE_DOUBLE;
+ opt.range_final_vert->description =
+ _("Range of final variogram: one value for anisotropic, two values for bivariate (hz and vert component)");
+ opt.range_final_vert->required = NO;
+ opt.range_final_vert->guisection = _("Final");
+ /* --------------------------------------------------------- */
- /* 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
- }
+ G_gisinit(argv[0]);
- // 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 (G_parser(argc, argv)) {
+ exit(EXIT_FAILURE);
+ }
- 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_pars.hz.name = var_pars.vert.name = var_pars.fin.name = opt.input->answer; // set name of variogram
- var_pars.hz.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
- var_pars.vert.dir = opt.var_dir_vert->answer ? DEG2RAD(atof(opt.var_dir_vert->answer)) : 0.; // Zenith angle
- var_pars.fin.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
-
- var_pars.hz.td = DEG2RAD(atof(opt.td_hz->answer)); // Angle of variogram processing
-
- if (opt.nL->answer) { // Test if nL have been set up (optional)
- if (var_pars.hz.nLag < 1) // Invalid value
- G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
- else {
- var_pars.hz.nLag = atof(opt.nL->answer);
- }
- }
-
- if (opt.form_file->answer) { // Plotting variogram
- set_gnuplot(opt.form_file->answer, &var_pars.hz);
- set_gnuplot(opt.form_file->answer, &var_pars.vert);
- set_gnuplot(opt.form_file->answer, &var_pars.fin);
- }
- else {
- strcpy(var_pars.hz.term, "");
- strcpy(var_pars.vert.term, "");
- strcpy(var_pars.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);
+ /* Get parameters from the parser */
+ if (strcmp(opt.phase->answer, "initial") == 0) {
+ xD.phase = 0; // estimate range
}
- G_fatal_error(_("You should mark either univariate, or bivariate variogram, not both of them..."));
- } // error
+ 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
+ }
- /* ---------------------------------------------------------- */
- 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);
+ // 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);
}
- 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 */
+ else
+ xD.report.write2file = FALSE;
- // 3D interpolation
- if (xD.i3 == TRUE) { // 3D interpolation:
- if (xD.v3 == FALSE) { // 2D input:
- if (!opt.zcol->answer) { // zcolumn not available:
- if (xD.report.write2file == TRUE) { // close report file
- fprintf(xD.report.fp, "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_pars.vert.nLag < 1) { // Invalid value
- G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
- }
- else {
- var_pars.vert.nLag = atof(opt.nZ->answer);
- }
- } // end if nZ->answer
- } // end if 3D interpolation
-
- else { // 2D interpolation:
- G_warning(_("Not recommended to process for sparse or spatially heterogeneous data. The result can be inaccurate - trying to solve asap..."));
- var_pars.vert.nLag = -1; // abs will be used in next steps
- if (xD.v3 == TRUE) {
- if (xD.report.write2file == TRUE) { // close report file
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("2D interpolation based on 3D input has been temporarily disabled. Please select another option..."));
+ 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;
+ }
- field = Vect_get_field_number(&map, opt.field->answer);
- if (xD.report.write2file == TRUE)
- write2file_basics(&xD, &opt);
- /* ---------------------------------------------------------- */
+ var_pars.hz.name = var_pars.vert.name = var_pars.fin.name = opt.input->answer; // set name of variogram
+ var_pars.hz.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
+ var_pars.vert.dir = opt.var_dir_vert->answer ? DEG2RAD(atof(opt.var_dir_vert->answer)) : 0.; // Zenith angle
+ var_pars.fin.dir = opt.var_dir_hz->answer ? DEG2RAD(atof(opt.var_dir_hz->answer)) : 0.; // Azimuth
- /* Get... */
- get_region_pars(&xD, ®); // ... region parameters
- read_points(&map, ®, &pnts, &xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
- pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
+ var_pars.hz.td = DEG2RAD(atof(opt.td_hz->answer)); // Angle of variogram processing
- /* Estimate 2D/3D variogram */
- switch (xD.phase) {
- case 0: // initial phase
- var_pars.hz.type = 0; // horizontal variogram
- if (xD.i3 == TRUE) { // 3D interpolation:
- var_pars.vert.type = 1; // vertical variogram
+ if (opt.nL->answer) { // Test if nL have been set up (optional)
+ if (var_pars.hz.nLag < 1) // Invalid value
+ G_message(_("Number of horizontal pieces must be at least 1. Default value will be used..."));
+ else {
+ var_pars.hz.nLag = atof(opt.nL->answer);
+ }
}
-
- variogram_restricts(&xD, &pnts, &var_pars.hz); // estimate lag size and number of lags (hz)
- if (xD.i3 == TRUE) { // 3D interpolation:
- variogram_restricts(&xD, &pnts, &var_pars.vert); // estimate lag size and number of lags (vert)
- if (var_pars.vert.nLag > var_pars.hz.nLag) { // more lags in vertical than in horizontal direction:
- var_pars.vert.nLag = var_pars.hz.nLag; // set the numbers to be equal
- var_pars.vert.lag = var_pars.vert.max_dist / var_pars.vert.nLag; // modify lag size
- }
+ if (opt.form_file->answer) { // Plotting variogram
+ set_gnuplot(opt.form_file->answer, &var_pars.hz);
+ set_gnuplot(opt.form_file->answer, &var_pars.vert);
+ set_gnuplot(opt.form_file->answer, &var_pars.fin);
}
- E_variogram(0, &xD, &pnts, &var_pars); // horizontal variogram (for both 2D and 3D interpolation)
- if (xD.i3 == TRUE) { // 3D interpolation:
- E_variogram(1, &xD, &pnts, &var_pars); // vertical variogram
- G_message(_("You may continue to computing theoretical variograms (middle phase)..."));
- }
else {
- G_message(_("You may continue to computing theoretical variograms (final phase)..."));
+ strcpy(var_pars.hz.term, "");
+ strcpy(var_pars.vert.term, "");
+ strcpy(var_pars.fin.term, "");
}
- goto end;
- case 1:
- read_tmp_vals("variogram_hz_tmp.txt", &var_pars.hz, &xD); // read properties of horizontal variogram from temp file
- read_tmp_vals("variogram_vert_tmp.txt", &var_pars.vert, &xD); // read properties of vertical variogram from temp file
+ 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
- if (xD.report.name) { // report file available:
- xD.report.write2file = TRUE;
- xD.report.fp = fopen(xD.report.name, "a");
- }
+ /* ---------------------------------------------------------- */
+ Vect_set_open_level(2); // Open input vector map
- T_variogram(0, TRUE, opt, &var_pars.hz, &xD.report); // compute theoretical variogram - hz
- T_variogram(1, TRUE, opt, &var_pars.vert, &xD.report); // compute theoretical variogram - vert
+ 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);
+ /* ---------------------------------------------------------- */
- /* compare range of hz and vert variogram:
- - if the difference is greater than 5% -> bivariate variogram
- - if the difference is smaller than 5% - anisotropic variogram
+ /* Perform 2D or 3D interpolation ? */
+ xD.i3 = flg.d23->answer ? FALSE : TRUE; // 2D/3D interpolation
+ xD.v3 = Vect_is_3d(&map); // 2D/3D layer
- => choose variogram type:
- */
+ /* 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 */
- sill_compare(&xD, &flg, &var_pars, &pnts);
- variogram_restricts(&xD, &pnts, &var_pars.fin);
+ // 3D interpolation
+ if (xD.i3 == TRUE) { // 3D interpolation:
+ if (xD.v3 == FALSE) { // 2D input:
+ if (!opt.zcol->answer) { // zcolumn not available:
+ if (xD.report.write2file == TRUE) { // close report file
+ fprintf(xD.report.fp,
+ "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_pars.vert.nLag < 1) { // Invalid value
+ G_message(_("Number of vertical pieces must be at least 1. Default value will be used..."));
+ }
+ else {
+ var_pars.vert.nLag = atof(opt.nZ->answer);
+ }
+ } // end if nZ->answer
+ } // end if 3D interpolation
- E_variogram(var_pars.fin.type, &xD, &pnts, &var_pars);
- G_message(_("You may continue to interpolate values (final phase)..."));
- goto end;
-
- case 2: // final phase:
- // Module should crash if:
- if (!opt.output->answer) { // output name not available:
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up name of output layer..."));
+ else { // 2D interpolation:
+ G_warning(_("Not recommended to process for sparse or spatially heterogeneous data. The result can be inaccurate - trying to solve asap..."));
+ var_pars.vert.nLag = -1; // abs will be used in next steps
+ if (xD.v3 == TRUE) {
+ if (xD.report.write2file == TRUE) { // close report file
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("2D interpolation based on 3D input has been temporarily disabled. Please select another option..."));
+ }
}
-
- // read properties of final variogram from temp file
- if (xD.i3 == TRUE) { // 3D kriging:
- read_tmp_vals("variogram_final_tmp.txt", &var_pars.fin, &xD);
- }
- else { // 2D kriging
- read_tmp_vals("variogram_hz_tmp.txt", &var_pars.fin, &xD);
- }
- if (xD.report.name) { // if report name available:
- xD.report.write2file = TRUE;
- xD.report.fp = fopen(xD.report.name, "a");
- if (xD.report.fp == NULL) // ... the file does not exist:
- G_fatal_error(_("Cannot open the file..."));
- }
-
- // check variogram settings
- if (var_pars.fin.type == 2 && strcmp(opt.function_var_final->answer, "linear") != 0) { // bivariate nonlinear variogram:
- // just one function type is set up (none or both should be)
- if (!(opt.function_var_final->answer && opt.function_var_final_vert->answer) &&
- (opt.function_var_final->answer || opt.function_var_final_vert->answer)) {
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("If you wish to specify components of bivariate variogram please set up function type for both of them..."));
- }
-
- // if both of the function types are set up:
- if (opt.function_var_final->answer && opt.function_var_final_vert->answer) {
- if (!opt.nugget_final->answer) { // horizontal nugget effect should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up nugget effect of horizontal component of bivariate variogram..."));
- }
+ field = Vect_get_field_number(&map, opt.field->answer);
+ if (xD.report.write2file == TRUE)
+ write2file_basics(&xD, &opt);
+ /* ---------------------------------------------------------- */
- if (!opt.nugget_final_vert->answer) { // vertical nugget effect should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up nugget effect of vertical component of bivariate variogram..."));
- }
+ /* Get... */
+ get_region_pars(&xD, ®); // ... region parameters
+ read_points(&map, ®, &pnts, &xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
+ pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
- if (!opt.range_final->answer) { // horizontal range should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up range of horizontal component of bivariate variogram..."));
- }
+ /* Estimate 2D/3D variogram */
+ switch (xD.phase) {
+ case 0: // initial phase
+ var_pars.hz.type = 0; // horizontal variogram
+ if (xD.i3 == TRUE) { // 3D interpolation:
+ var_pars.vert.type = 1; // vertical variogram
+ }
- if (!opt.range_final_vert->answer) { // vertical range should not be missing
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up range of vertical component of bivariate variogram..."));
- }
- }
- }
+ variogram_restricts(&xD, &pnts, &var_pars.hz); // estimate lag size and number of lags (hz)
- else { // univariate variogram
- if (xD.i3 == TRUE && (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0)) { // anisotropic 3D:
- if (opt.function_var_final_vert->answer || atof(opt.nugget_final_vert->answer) != 0. || opt.range_final_vert->answer) { // vertical settings available:
- G_fatal_error(_("Not necessary to set up vertical components properties. Anisotropic variogram will be used..."));
- } // end if vert settings available
- } // end if 3D
+ if (xD.i3 == TRUE) { // 3D interpolation:
+ variogram_restricts(&xD, &pnts, &var_pars.vert); // estimate lag size and number of lags (vert)
+ if (var_pars.vert.nLag > var_pars.hz.nLag) { // more lags in vertical than in horizontal direction:
+ var_pars.vert.nLag = var_pars.hz.nLag; // set the numbers to be equal
+ var_pars.vert.lag = var_pars.vert.max_dist / var_pars.vert.nLag; // modify lag size
+ }
+ }
+ E_variogram(0, &xD, &pnts, &var_pars); // horizontal variogram (for both 2D and 3D interpolation)
+ if (xD.i3 == TRUE) { // 3D interpolation:
+ E_variogram(1, &xD, &pnts, &var_pars); // vertical variogram
+ G_message(_("You may continue to computing theoretical variograms (middle phase)..."));
+ }
+ else {
+ G_message(_("You may continue to computing theoretical variograms (final phase)..."));
+ }
+ goto end;
- if (!opt.function_var_final->answer) { // missing function
- opt.function_var_final->answer = "linear";
- G_message(_("Linear variogram will be computed..."));
- }
+ case 1:
+ read_tmp_vals("variogram_hz_tmp.txt", &var_pars.hz, &xD); // read properties of horizontal variogram from temp file
+ read_tmp_vals("variogram_vert_tmp.txt", &var_pars.vert, &xD); // read properties of vertical variogram from temp file
- if (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0) { // nonlinear and not parabolic:
- if (!opt.nugget_final->answer) { // missing horizontal nugget effect
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up nugget effect..."));
- } // end if nugget missing
+ if (xD.report.name) { // report file available:
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ }
- if (!opt.range_final->answer) { // missing horizontal range:
- if (xD.report.write2file == TRUE) { // report file available:
- fprintf(xD.report.fp, "Error (see standard output). Process killed...");
- fclose(xD.report.fp);
- }
- G_fatal_error(_("Please set up range..."));
- } // end if range missing
+ T_variogram(0, TRUE, opt, &var_pars.hz, &xD.report); // compute theoretical variogram - hz
+ T_variogram(1, TRUE, opt, &var_pars.vert, &xD.report); // compute theoretical variogram - vert
- if (opt.sill_final->answer) { // missing horizontal range:
- var_pars.fin.sill = atof(opt.sill_final->answer);
- } // end if sill has been changed by the user
- } // end nonlinear and not parabolic variogram
- } // end if univariate variogram (2D or 3D)
-
- out.name = opt.output->answer; // Output layer name
-
- if (var_pars.fin.type == 3) { // if variogram is anisotropic:
- geometric_anisotropy(&xD, &pnts); // exaggerate z coord and rebuild the spatial index
+ /* compare range of hz and vert variogram:
+ - if the difference is greater than 5% -> bivariate variogram
+ - if the difference is smaller than 5% - anisotropic variogram
+
+ => choose variogram type:
+ */
+
+ sill_compare(&xD, &flg, &var_pars, &pnts);
+ variogram_restricts(&xD, &pnts, &var_pars.fin);
+
+ E_variogram(var_pars.fin.type, &xD, &pnts, &var_pars);
+ G_message(_("You may continue to interpolate values (final phase)..."));
+ goto end;
+
+ case 2: // final phase:
+ // Module should crash if:
+ if (!opt.output->answer) { // output name not available:
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up name of output layer..."));
+ }
+
+ // read properties of final variogram from temp file
+ if (xD.i3 == TRUE) { // 3D kriging:
+ read_tmp_vals("variogram_final_tmp.txt", &var_pars.fin, &xD);
+ }
+ else { // 2D kriging
+ read_tmp_vals("variogram_hz_tmp.txt", &var_pars.fin, &xD);
+ }
+
+ if (xD.report.name) { // if report name available:
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ if (xD.report.fp == NULL) // ... the file does not exist:
+ G_fatal_error(_("Cannot open the file..."));
+ }
+
+ // check variogram settings
+ if (var_pars.fin.type == 2 && strcmp(opt.function_var_final->answer, "linear") != 0) { // bivariate nonlinear variogram:
+ // just one function type is set up (none or both should be)
+ if (!
+ (opt.function_var_final->answer &&
+ opt.function_var_final_vert->answer) &&
+ (opt.function_var_final->answer ||
+ opt.function_var_final_vert->answer)) {
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("If you wish to specify components of bivariate variogram please set up function type for both of them..."));
+ }
+
+ // if both of the function types are set up:
+ if (opt.function_var_final->answer &&
+ opt.function_var_final_vert->answer) {
+ if (!opt.nugget_final->answer) { // horizontal nugget effect should not be missing
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up nugget effect of horizontal component of bivariate variogram..."));
+ }
+
+ if (!opt.nugget_final_vert->answer) { // vertical nugget effect should not be missing
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up nugget effect of vertical component of bivariate variogram..."));
+ }
+
+ if (!opt.range_final->answer) { // horizontal range should not be missing
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up range of horizontal component of bivariate variogram..."));
+ }
+
+ if (!opt.range_final_vert->answer) { // vertical range should not be missing
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up range of vertical component of bivariate variogram..."));
+ }
+ }
+ }
+
+ else { // univariate variogram
+ if (xD.i3 == TRUE && (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0)) { // anisotropic 3D:
+ if (opt.function_var_final_vert->answer || atof(opt.nugget_final_vert->answer) != 0. || opt.range_final_vert->answer) { // vertical settings available:
+ G_fatal_error(_("Not necessary to set up vertical components properties. Anisotropic variogram will be used..."));
+ } // end if vert settings available
+ } // end if 3D
+
+ if (!opt.function_var_final->answer) { // missing function
+ opt.function_var_final->answer = "linear";
+ G_message(_("Linear variogram will be computed..."));
+ }
+
+ if (strcmp(opt.function_var_final->answer, "linear") != 0 && strcmp(opt.function_var_final->answer, "parabolic") != 0) { // nonlinear and not parabolic:
+ if (!opt.nugget_final->answer) { // missing horizontal nugget effect
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up nugget effect..."));
+ } // end if nugget missing
+
+ if (!opt.range_final->answer) { // missing horizontal range:
+ if (xD.report.write2file == TRUE) { // report file available:
+ fprintf(xD.report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD.report.fp);
+ }
+ G_fatal_error(_("Please set up range..."));
+ } // end if range missing
+
+ if (opt.sill_final->answer) { // missing horizontal range:
+ var_pars.fin.sill = atof(opt.sill_final->answer);
+ } // end if sill has been changed by the user
+ } // end nonlinear and not parabolic variogram
+ } // end if univariate variogram (2D or 3D)
+
+ out.name = opt.output->answer; // Output layer name
+
+ if (var_pars.fin.type == 3) { // if variogram is anisotropic:
+ geometric_anisotropy(&xD, &pnts); // exaggerate z coord and rebuild the spatial index
+ }
+
+ T_variogram(var_pars.fin.type, xD.i3, opt, &var_pars.fin, &xD.report); // compute theoretical variogram
+
+ G_debug(0, "sill: %f", var_pars.fin.sill);
+ break;
}
-
- T_variogram(var_pars.fin.type, xD.i3, opt, &var_pars.fin, &xD.report); // compute theoretical variogram
- G_debug(0,"sill: %f", var_pars.fin.sill);
- break;
- }
+ /* Ordinary kriging (including 2D/3D raster creation) */
+ ordinary_kriging(&xD, ®, &pnts, &var_pars, &out);
+ /* ---------------------------------------------------------- */
- /* Ordinary kriging (including 2D/3D raster creation) */
- ordinary_kriging(&xD, ®, &pnts, &var_pars, &out);
- /* ---------------------------------------------------------- */
+ end:
- end:
-
- G_free(pnts.r);
- Vect_close(&map); // Close vector map
- exit(EXIT_SUCCESS);
+ G_free(pnts.r);
+ Vect_close(&map); // Close vector map
+ exit(EXIT_SUCCESS);
}
-
Modified: grass-addons/grass7/vector/v.kriging/quantile.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/quantile.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/quantile.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -4,9 +4,9 @@
struct bin
{
- unsigned long origin;
- double minimum, maximum;
- int base, count;
+ unsigned long origin;
+ double minimum, maximum;
+ int base, count;
};
static double minimum, maximum;
@@ -24,96 +24,97 @@
static inline int get_slot(double c)
{
- int i = (int) floor((c - minimum) / slot_size);
-
+ int i = (int)floor((c - minimum) / slot_size);
+
if (i < 0)
- i = 0;
+ i = 0;
if (i > num_slots - 1)
- i = num_slots - 1;
+ i = num_slots - 1;
return i;
}
double get_quantile(int n)
{
- return (double)total * quants[n]; // # of lower values
+ return (double)total *quants[n]; // # of lower values
}
void get_slot_counts(int n, double *data)
{
- int i;
-
- total = 0;
- for (i = 0; i < n; i++) {
- int j;
+ int i;
- // todo nan
- j = get_slot(data[i]);
+ total = 0;
+ for (i = 0; i < n; i++) {
+ int j;
- slots[j]++; // rise value of j-th slot
- total++;
- }
- //G_percent(i, n, 2);
+ // todo nan
+ j = get_slot(data[i]);
+
+ slots[j]++; // rise value of j-th slot
+ total++;
+ }
+ //G_percent(i, n, 2);
}
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
+ 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
+ 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
+ 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;
+ if (accum2 > next) { // # of values is greater than percentile
+ struct bin *b = &bins[bin]; // make 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);
+ slot_bins[slot] = ++bin;
- while (accum2 > next)
- next = get_quantile(++quant);
+ 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);
- num_values += count;
+ while (accum2 > next)
+ next = get_quantile(++quant);
+
+ num_values += count;
+ }
+
+ accum = accum2;
}
- accum = accum2;
- }
-
- num_bins = bin;
+ num_bins = bin;
}
void fill_bins(int n, double *data)
{
- int i;
-
- for (i = 0; i < n; i++) {
- int j, bin;
- struct bin *b;
+ int i;
- // todo nan
+ for (i = 0; i < n; i++) {
+ int j, bin;
+ struct bin *b;
- j = get_slot(data[i]);
+ // todo nan
- if (!slot_bins[j])
- continue;
+ j = get_slot(data[i]);
- bin = slot_bins[j] - 1;
- b = &bins[bin];
+ if (!slot_bins[j])
+ continue;
- values[b->base + b->count++] = data[i];
- }
+ bin = slot_bins[j] - 1;
+ b = &bins[bin];
- //G_percent(i, n, 2);
+ values[b->base + b->count++] = data[i];
+ }
+
+ //G_percent(i, n, 2);
}
int compare(const void *aa, const void *bb)
@@ -122,9 +123,9 @@
double b = *(const double *)bb;
if (a < b)
- return -1;
+ return -1;
if (a > b)
- return 1;
+ return 1;
return 0;
}
@@ -133,10 +134,10 @@
int bin;
for (bin = 0; bin < num_bins; bin++) {
- struct bin *b = &bins[bin];
+ struct bin *b = &bins[bin];
- qsort(&values[b->base], b->count, sizeof(double), compare);
- //G_percent(bin, num_bins, 2);
+ qsort(&values[b->base], b->count, sizeof(double), compare);
+ //G_percent(bin, num_bins, 2);
}
//G_percent(bin, num_bins, 2);
}
@@ -148,76 +149,76 @@
int quant;
for (quant = 0; quant < num_quants; quant++) {
- struct bin *b;
- double next = get_quantile(quant);
- double k, v;
- int i0, i1;
+ 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;
- }
+ 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 (bin < num_bins) {
+ k = next - b->origin;
+ i0 = (int)floor(k);
+ i1 = (int)ceil(k);
- if (i1 > b->count - 1)
- i1 = b->count - 1;
+ 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;
+ 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]);
+ *quantile = v;
+ prev_v = v;
+ if (report != NULL && report->write2file == TRUE)
+ fprintf(report->fp, "%f:%f:%f\n", prev_v, v, quants[quant]);
}
}
double quantile(double q, int n, double *data, struct write *report)
{
- int i, recode=FALSE;
- double quantile;
+ int i, recode = FALSE;
+ double quantile;
- num_slots = 1000000; // # of slots
-
- num_quants = 1;
- quants = (double *) G_malloc(num_quants * sizeof(double));
- quants[0] = q; // compute quantile for 0.05
+ num_slots = 1000000; // # of slots
- // 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));
+ num_quants = 1;
+ quants = (double *)G_malloc(num_quants * sizeof(double));
+ quants[0] = q; // compute quantile for 0.05
- // get min and max of the data
- minimum = data[0];
- maximum = data[0];
- for (i=0; i<n; i++) {
- minimum = MIN(data[i], minimum);
- maximum = MAX(data[i], maximum);
- }
+ // 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));
- slot_size = (maximum - minimum) / num_slots;
- get_slot_counts(n, data); // # of data in each slot
+ // get min and max of the data
+ minimum = data[0];
+ maximum = data[0];
+ for (i = 0; i < n; i++) {
+ minimum = MIN(data[i], minimum);
+ maximum = MAX(data[i], maximum);
+ }
- bins = (struct bin *) G_calloc(num_quants, sizeof(struct bin));
- initialize_bins();
- G_free(slots);
+ slot_size = (maximum - minimum) / num_slots;
+ get_slot_counts(n, data); // # of data in each slot
- values = (double *) G_calloc(num_values, sizeof(double));
- fill_bins(n, data);
- G_free(slot_bins);
+ bins = (struct bin *)G_calloc(num_quants, sizeof(struct bin));
+ initialize_bins();
+ G_free(slots);
- sort_bins();
- compute_quantiles(recode, &quantile, report);
-
- return quantile;
+ 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;
}
Modified: grass-addons/grass7/vector/v.kriging/spatial_index.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/spatial_index.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/spatial_index.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -2,124 +2,117 @@
void insert_rectangle(int dim, int i, struct points *pnts)
{
- struct RTree_Rect *rect; // rectangle
- double *r; // pointer to the input coordinates
- r = &pnts->r[3 * i];
+ struct RTree_Rect *rect; // rectangle
+ double *r; // pointer to the input coordinates
- switch (dim) {
- case 3: // 3D:
- rect = RTreeAllocRect(pnts->R_tree); // allocate the rectangle
- RTreeSetRect3D(rect, pnts->R_tree,
- *r, *r,
- *(r+1), *(r+1),
- *(r+2), *(r+2)); // set 3D coordinates
- RTreeInsertRect(rect, i + 1, pnts->R_tree); // insert rectangle to the spatial index
- break;
- case 2: // 2D:
- rect = RTreeAllocRect(pnts->Rtree_hz); // allocate the rectangle
- RTreeSetRect2D(rect, pnts->Rtree_hz,
- *r, *r,
- *(r+1), *(r+1)); // set 2D coordinates
- RTreeInsertRect(rect, i + 1, pnts->Rtree_hz); // insert rectangle to the spatial index
- break;
- case 1: // 1D
- rect = RTreeAllocRect(pnts->Rtree_vert); // allocate the rectangle
- RTreeSetRect1D(rect, pnts->Rtree_vert,
- *(r+2), *(r+2)); // set 1D coordinates
- RTreeInsertRect(rect, i + 1, pnts->Rtree_vert); // insert rectangle to the spatial index
- break;
- }
-
- RTreeFreeRect(rect); // free rectangle memory
+ r = &pnts->r[3 * i];
+
+ switch (dim) {
+ case 3: // 3D:
+ rect = RTreeAllocRect(pnts->R_tree); // allocate the rectangle
+ RTreeSetRect3D(rect, pnts->R_tree, *r, *r, *(r + 1), *(r + 1), *(r + 2), *(r + 2)); // set 3D coordinates
+ RTreeInsertRect(rect, i + 1, pnts->R_tree); // insert rectangle to the spatial index
+ break;
+ case 2: // 2D:
+ rect = RTreeAllocRect(pnts->Rtree_hz); // allocate the rectangle
+ RTreeSetRect2D(rect, pnts->Rtree_hz, *r, *r, *(r + 1), *(r + 1)); // set 2D coordinates
+ RTreeInsertRect(rect, i + 1, pnts->Rtree_hz); // insert rectangle to the spatial index
+ break;
+ case 1: // 1D
+ rect = RTreeAllocRect(pnts->Rtree_vert); // allocate the rectangle
+ RTreeSetRect1D(rect, pnts->Rtree_vert, *(r + 2), *(r + 2)); // set 1D coordinates
+ RTreeInsertRect(rect, i + 1, pnts->Rtree_vert); // insert rectangle to the spatial index
+ break;
+ }
+
+ RTreeFreeRect(rect); // free rectangle memory
}
/* find neighbors in cubic (square) surroundings */
-struct ilist *find_NNs_within(int dim, double *search_pt, struct points *pnts, double max_dist, double max_dist_vert)
+struct ilist *find_NNs_within(int dim, double *search_pt, struct points *pnts,
+ double max_dist, double max_dist_vert)
{
- // local variables
- double *r = search_pt; // pointer to vector of the coordinates
+ // local variables
+ double *r = search_pt; // pointer to vector of the coordinates
- int NN_sum = 0; // # of NN
- double dist_step = max_dist; // distance iteration in case of empty closest surrounding of search point
- double dist_step_vert = max_dist_vert; // vertical distance iteration
- struct RTree_Rect *search; // search rectangle
+ int NN_sum = 0; // # of NN
+ double dist_step = max_dist; // distance iteration in case of empty closest surrounding of search point
+ double dist_step_vert = max_dist_vert; // vertical distance iteration
+ struct RTree_Rect *search; // search rectangle
- struct ilist *list;
- list = G_new_ilist();
+ struct ilist *list;
- switch (dim) {
- case 3:
- search = RTreeAllocRect(pnts->R_tree); // allocate new rectangle
- break;
- case 2:
- search = RTreeAllocRect(pnts->Rtree_hz); // allocate new rectangle
- break;
- case 1:
- search = RTreeAllocRect(pnts->Rtree_vert); // allocate new rectangle
- break;
- }
+ list = G_new_ilist();
- while (NN_sum < 2) { // TODO: this might be tricky - what # (density) is optimal?
switch (dim) {
- case 3: // 3D:
- RTreeSetRect3D(search, pnts->R_tree,
- *r - max_dist, *r + max_dist,
- *(r+1) - max_dist, *(r+1) + max_dist,
- *(r+2) - max_dist_vert, *(r+2) + max_dist_vert); // set up searching rectangle
- RTreeSearch2(pnts->R_tree, search, list); // search the nearest rectangle
- break;
- case 2: // 2D:
- RTreeSetRect2D(search, pnts->Rtree_hz,
- *r - max_dist, *r + max_dist,
- *(r+1) - max_dist, *(r+1) + max_dist); // set up searching rectangle
- RTreeSearch2(pnts->Rtree_hz, search, list); // search the nearest rectangle
- break;
- case 1: // 1D:
- RTreeSetRect1D(search, pnts->Rtree_vert,
- *(r+2) - max_dist_vert, *(r+2) + max_dist_vert); // set up searching rectangle
-
- RTreeSearch2(pnts->Rtree_vert, search, list); // search the nearest rectangle
- break;
+ case 3:
+ search = RTreeAllocRect(pnts->R_tree); // allocate new rectangle
+ break;
+ case 2:
+ search = RTreeAllocRect(pnts->Rtree_hz); // allocate new rectangle
+ break;
+ case 1:
+ search = RTreeAllocRect(pnts->Rtree_vert); // allocate new rectangle
+ break;
}
- NN_sum = list->n_values;
- if (NN_sum < 2) {
- G_warning(_("Point \"x=%f y=%f z=%f\" has less than 2 neighbours in its closest surrounding. The perimeter of the surrounding will be increased to include more neighbouring points"), *r, *(r+1), *(r+2));
- max_dist += dist_step;
- max_dist_vert += dist_step_vert;
- }
- } // end while: nonzero # of NNs
+ while (NN_sum < 2) { // TODO: this might be tricky - what # (density) is optimal?
+ switch (dim) {
+ case 3: // 3D:
+ RTreeSetRect3D(search, pnts->R_tree, *r - max_dist, *r + max_dist, *(r + 1) - max_dist, *(r + 1) + max_dist, *(r + 2) - max_dist_vert, *(r + 2) + max_dist_vert); // set up searching rectangle
+ RTreeSearch2(pnts->R_tree, search, list); // search the nearest rectangle
+ break;
+ case 2: // 2D:
+ RTreeSetRect2D(search, pnts->Rtree_hz, *r - max_dist, *r + max_dist, *(r + 1) - max_dist, *(r + 1) + max_dist); // set up searching rectangle
+ RTreeSearch2(pnts->Rtree_hz, search, list); // search the nearest rectangle
+ break;
+ case 1: // 1D:
+ RTreeSetRect1D(search, pnts->Rtree_vert, *(r + 2) - max_dist_vert, *(r + 2) + max_dist_vert); // set up searching rectangle
- RTreeFreeRect(search);
+ RTreeSearch2(pnts->Rtree_vert, search, list); // search the nearest rectangle
+ break;
+ }
- return list;
+ NN_sum = list->n_values;
+ if (NN_sum < 2) {
+ G_warning(_("Point \"x=%f y=%f z=%f\" has less than 2 neighbours in its closest surrounding. The perimeter of the surrounding will be increased to include more neighbouring points"),
+ *r, *(r + 1), *(r + 2));
+ max_dist += dist_step;
+ max_dist_vert += dist_step_vert;
+ }
+ } // end while: nonzero # of NNs
+
+ RTreeFreeRect(search);
+
+ return list;
}
/* find neighbors in spherical (circular) surroundings */
struct ilist *find_n_NNs(int dim, int i, struct points *pnts, int n)
{
- // local variables
- double max_dist0 = pnts->max_dist; // initial maximum distance
- struct ilist *list; // list of selected rectangles (points)
+ // local variables
+ double max_dist0 = pnts->max_dist; // initial maximum distance
+ struct ilist *list; // list of selected rectangles (points)
- int n_vals = 0; // # of the nearest neighbours (NN)
- int iter = 1; // iteration number (searching NN)
- double max_dist = max_dist0; // iterated search radius
- double *search;
- search = &pnts->r[3 * i];
+ int n_vals = 0; // # of the nearest neighbours (NN)
+ int iter = 1; // iteration number (searching NN)
+ double max_dist = max_dist0; // iterated search radius
+ double *search;
- while (n_vals < n) { // until some NN is found (2 points: identical and NN)
- list = find_NNs_within(dim, search, pnts, max_dist, -1);
- n_vals = list->n_values; // set up control variable
-
- if (n_vals < n) {
- iter += 1; // go to the next iteration
- max_dist = iter * max_dist0; // enlarge search radius
+ search = &pnts->r[3 * i];
+
+ while (n_vals < n) { // until some NN is found (2 points: identical and NN)
+ list = find_NNs_within(dim, search, pnts, max_dist, -1);
+ n_vals = list->n_values; // set up control variable
+
+ if (n_vals < n) {
+ iter += 1; // go to the next iteration
+ max_dist = iter * max_dist0; // enlarge search radius
+ }
}
- }
- // check spherical (circular) surrounding
- list = find_NNs_within(dim, search, pnts, sqrt(2.) * max_dist, -1);
+ // check spherical (circular) surrounding
+ list = find_NNs_within(dim, search, pnts, sqrt(2.) * max_dist, -1);
- return list;
+ return list;
}
Modified: grass-addons/grass7/vector/v.kriging/stat.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/stat.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/stat.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -2,51 +2,51 @@
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;
+ 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 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++;
- }
+ // 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.;
+ 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);
+ // 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");
+ 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");
}
Modified: grass-addons/grass7/vector/v.kriging/utils.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/utils.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -3,853 +3,920 @@
#define GNUPLOT "gnuplot -persist"
int cmpVals(const void *v1, const void *v2)
- {
+{
int *p1, *p2;
p1 = (int *)v1;
p2 = (int *)v2;
if (p1[0] > p2[0])
- return 1;
+ return 1;
else if (p1[0] < p2[0])
- return -1;
+ return -1;
else
- return 0;
- }
+ return 0;
+}
-void correct_indices(struct ilist *list, double *r0, struct points *pnts, struct parameters *var_pars)
+void correct_indices(struct ilist *list, double *r0, struct points *pnts,
+ struct parameters *var_pars)
{
- int i;
- int n = list->n_values;
- int *vals = list->value;
- double dx, dy, dz;
- int type = var_pars->type;
-
- int n_new = 0; // # of effective indices (not identical points, nor too far)
- int *newvals, *save; // effective indices
- newvals = (int *) G_malloc(n * sizeof(int));
- save = &newvals[0];
+ int i;
+ int n = list->n_values;
+ int *vals = list->value;
+ double dx, dy, dz;
+ int type = var_pars->type;
- double *r, sqDist_i;
+ int n_new = 0; // # of effective indices (not identical points, nor too far)
+ int *newvals, *save; // effective indices
- for (i=0; i<n; i++) {
- *vals -= 1;
- r = &pnts->r[3 * *vals];
+ newvals = (int *)G_malloc(n * sizeof(int));
+ save = &newvals[0];
- dx = *r - *r0;
- dy = *(r+1) - *(r0+1);
- dz = *(r+2) - *(r0+2);
- sqDist_i = SQUARE(dx) + SQUARE(dy) + SQUARE(dz);
-
- if (type != 2 && (n == 1 || (sqDist_i != 0.))) {
- *save = *vals;
- n_new++;
- save++;
+ double *r, sqDist_i;
+
+ for (i = 0; i < n; i++) {
+ *vals -= 1;
+ r = &pnts->r[3 * *vals];
+
+ dx = *r - *r0;
+ dy = *(r + 1) - *(r0 + 1);
+ dz = *(r + 2) - *(r0 + 2);
+ sqDist_i = SQUARE(dx) + SQUARE(dy) + SQUARE(dz);
+
+ if (type != 2 && (n == 1 || (sqDist_i != 0.))) {
+ *save = *vals;
+ n_new++;
+ save++;
+ }
+ vals++;
}
- vals++;
- }
- if (type != 2 && n_new < n) {
- list->n_values = n_new;
- list->value = (int *) G_realloc(list->value, n_new * sizeof(int));
- memcpy(list->value, newvals, n_new * sizeof(int));
- }
+ if (type != 2 && n_new < n) {
+ list->n_values = n_new;
+ list->value = (int *)G_realloc(list->value, n_new * sizeof(int));
+ memcpy(list->value, newvals, n_new * sizeof(int));
+ }
- G_free(newvals);
+ G_free(newvals);
- return;
+ return;
}
// make coordinate triples xyz
void triple(double x, double y, double z, double *triple)
{
- double *t;
- t = &triple[0];
+ double *t;
- *t = x;
- *(t+1) = y;
- *(t+2) = z;
+ t = &triple[0];
- return;
+ *t = x;
+ *(t + 1) = y;
+ *(t + 2) = z;
+
+ return;
}
// compute coordinate differences
void coord_diff(int i, int j, double *r, double *dr)
{
- int k = 3*i, l = 3*j;
- double *rk, *rl, *drt;
- rk = &r[k];
- rl = &r[l];
- drt = &dr[0];
+ int k = 3 * i, l = 3 * j;
+ double *rk, *rl, *drt;
- if ((*rk == 0.) && (*(rk+1) == 0.) && (*(rk+2) == 0.)) {
- G_fatal_error(_("Coordinates of point no. %d are zeros."),i);
- }
+ rk = &r[k];
+ rl = &r[l];
+ drt = &dr[0];
- *drt = *rl - *rk; // dx
- *(drt+1) = *(rl+1) - *(rk+1); // dy
- *(drt+2) = *(rl+2) - *(rk+2); // dz
+ if ((*rk == 0.) && (*(rk + 1) == 0.) && (*(rk + 2) == 0.)) {
+ G_fatal_error(_("Coordinates of point no. %d are zeros."), i);
+ }
- return;
+ *drt = *rl - *rk; // dx
+ *(drt + 1) = *(rl + 1) - *(rk + 1); // dy
+ *(drt + 2) = *(rl + 2) - *(rk + 2); // dz
+
+ return;
}
// compute horizontal radius from coordinate differences
double radius_hz_diff(double *dr)
{
- double rds;
- rds = SQUARE(dr[0]) + SQUARE(dr[1]);
+ double rds;
- return rds;
+ rds = SQUARE(dr[0]) + SQUARE(dr[1]);
+
+ return rds;
}
// compute zenith angle
double zenith_angle(double dr[3])
{
- double zenith;
- zenith = atan2(dr[2], sqrt(SQUARE(dr[0]) + SQUARE(dr[1])));
+ double zenith;
- return zenith;
+ zenith = atan2(dr[2], sqrt(SQUARE(dr[0]) + SQUARE(dr[1])));
+
+ return zenith;
}
// compute size of the lag
-double lag_size(int direction, struct points *pnts, struct parameters *var_pars, struct write *report)
+double lag_size(int direction, struct points *pnts,
+ struct parameters *var_pars, struct write *report)
{
- // local variables
- int n = pnts->n; // # of input points
- double *r = pnts->r; // xyz of input points
- int type = var_pars->type;
- double max_dist = type == 2 ? var_pars->horizontal.max_dist : var_pars->max_dist; // maximum horizontal distance (or vertical for vertical variogram)
- double max_dist_vert = type == 2 ? var_pars->vertical.max_dist : var_pars->max_dist; // maximum vertical distance (just for bivariate variogram)
+ // local variables
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz of input points
+ int type = var_pars->type;
+ double max_dist = type == 2 ? var_pars->horizontal.max_dist : var_pars->max_dist; // maximum horizontal distance (or vertical for vertical variogram)
+ double max_dist_vert = type == 2 ? var_pars->vertical.max_dist : var_pars->max_dist; // maximum vertical distance (just for bivariate variogram)
- int i, j; // loop indices
- int n_vals; // number of the nearest neighbors (NN)
- int *j_vals; // indices of NN (NOT sorted by distance)
- struct ilist *list; // list of NN
-
- double dr[3]; // coordinate differences
- double d_min = SQUARE(max_dist) + SQUARE(max_dist_vert); // the shortest distance of the 5%-th NN
- double lagNN; // square root of d_min -> lag
- double dist2; // square of the distances between search point and NN
+ int i, j; // loop indices
+ int n_vals; // number of the nearest neighbors (NN)
+ int *j_vals; // indices of NN (NOT sorted by distance)
+ struct ilist *list; // list of NN
- double *search; // search point
- int perc5; // 5% from total number of NN
- int add_ident; // number of identical points to add to perc5
-
- for (i=0; i < n; i++) { // for each input point...
-
- add_ident = 0; // number of points identical to search point equals to zero
- search = &pnts->r[3 * i];
+ double dr[3]; // coordinate differences
+ double d_min = SQUARE(max_dist) + SQUARE(max_dist_vert); // the shortest distance of the 5%-th NN
+ double lagNN; // square root of d_min -> lag
+ double dist2; // square of the distances between search point and NN
- switch (direction) {
- case 12: // horizontal variogram
- list = find_NNs_within(2, search, pnts, max_dist, -1);
- break;
- case 3: // vertical variogram
- list = find_NNs_within(1, search, pnts, max_dist, max_dist_vert);
- break;
- default: // anisotropic and bivariate variogram
- list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
- break;
- }
- n_vals = list->n_values; // number of overlapping rectangles
-
- // find number of points identical to the search point
- if (n_vals > 0) {
- j_vals = &list->value[0]; // indices of overlapping rectangles (note: increased by 1)
+ double *search; // search point
+ int perc5; // 5% from total number of NN
+ int add_ident; // number of identical points to add to perc5
- for (j=0; j < n_vals; j++) { // for each overlapping rectangle:
- coord_diff(i, (*j_vals-1), pnts->r, dr); // compute coordinate differences
- switch (direction) { // and count those which equal to zero
- case 12: // horizontal
- if (dr[0] == 0. && dr[1] == 0.) {
- add_ident++;
- }
- break;
- case 3: // vertical
- if (dr[2] == 0.) {
- add_ident++;
- }
- break;
- default: // aniso / bivar
- if(dr[0] == 0. && dr[1] == 0. && dr[2] == 0.) {
- add_ident++;
- }
- break;
- } // end switch: test direction
- j_vals++; // index of the next NN
- } // end j for loop
- } //end if: n_vals > 0
-
- else {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("Error: There are no nearest neighbours of point %d..."), i);
- } // end error
+ for (i = 0; i < n; i++) { // for each input point...
- perc5 = (int) round(0.05 * n) + add_ident; // set up 5% increased by number of identical points
+ add_ident = 0; // number of points identical to search point equals to zero
+ search = &pnts->r[3 * i];
- coord_diff(i, perc5, pnts->r, dr); // coordinate differences of search point and 5%th NN
- switch (direction) { // compute squared distance for particular direction
- case 12: // horizontal
- dist2 = SQUARE(dr[0]) + SQUARE(dr[1]);
- break;
- case 3: // vertical
- dist2 = SQUARE(dr[2]);
- break;
- default: // aniso / bivar
- dist2 = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
- break;
- }
+ switch (direction) {
+ case 12: // horizontal variogram
+ list = find_NNs_within(2, search, pnts, max_dist, -1);
+ break;
+ case 3: // vertical variogram
+ list = find_NNs_within(1, search, pnts, max_dist, max_dist_vert);
+ break;
+ default: // anisotropic and bivariate variogram
+ list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
+ break;
+ }
+ n_vals = list->n_values; // number of overlapping rectangles
- if (dist2 != 0.) { // and compare it with minimum distance between search point and 5% NN
- d_min = MIN(d_min, dist2);
- }
-
- r += 3; // go to next search point
- } // end i for loop: distance of i-th search pt and 5%-th NN
+ // find number of points identical to the search point
+ if (n_vals > 0) {
+ j_vals = &list->value[0]; // indices of overlapping rectangles (note: increased by 1)
- lagNN = sqrt(d_min);
+ for (j = 0; j < n_vals; j++) { // for each overlapping rectangle:
+ coord_diff(i, (*j_vals - 1), pnts->r, dr); // compute coordinate differences
+ switch (direction) { // and count those which equal to zero
+ case 12: // horizontal
+ if (dr[0] == 0. && dr[1] == 0.) {
+ add_ident++;
+ }
+ break;
+ case 3: // vertical
+ if (dr[2] == 0.) {
+ add_ident++;
+ }
+ break;
+ default: // aniso / bivar
+ if (dr[0] == 0. && dr[1] == 0. && dr[2] == 0.) {
+ add_ident++;
+ }
+ break;
+ } // end switch: test direction
+ j_vals++; // index of the next NN
+ } // end j for loop
+ } //end if: n_vals > 0
- G_free_ilist(list); // free list memory
-
- return lagNN;
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Error: There are no nearest neighbours of point %d..."),
+ i);
+ } // end error
+
+ perc5 = (int)round(0.05 * n) + add_ident; // set up 5% increased by number of identical points
+
+ coord_diff(i, perc5, pnts->r, dr); // coordinate differences of search point and 5%th NN
+ switch (direction) { // compute squared distance for particular direction
+ case 12: // horizontal
+ dist2 = SQUARE(dr[0]) + SQUARE(dr[1]);
+ break;
+ case 3: // vertical
+ dist2 = SQUARE(dr[2]);
+ break;
+ default: // aniso / bivar
+ dist2 = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]);
+ break;
+ }
+
+ if (dist2 != 0.) { // and compare it with minimum distance between search point and 5% NN
+ d_min = MIN(d_min, dist2);
+ }
+
+ r += 3; // go to next search point
+ } // end i for loop: distance of i-th search pt and 5%-th NN
+
+ lagNN = sqrt(d_min);
+
+ G_free_ilist(list); // free list memory
+
+ return lagNN;
}
double linear(double x, double a, double b)
{
- double y;
- y = a * x + b;
+ double y;
- return 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
+ double y;
- return 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;
+ double y;
+ double ratio;
- return y;
+ 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));
+ double y;
+ double c2 = SQUARE(c);
+ double ratio;
- return y;
+ 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;
+ 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;
}
void optimize(double *lag, int *nLag, double max)
{
- if (*nLag > 20) {
- *nLag = 20;
- *lag = max / *nLag;
- }
+ if (*nLag > 20) {
+ *nLag = 20;
+ *lag = max / *nLag;
+ }
- return;
+ return;
}
// maximal horizontal and vertical distance to restrict variogram computing
-void variogram_restricts(struct int_par *xD, struct points *pnts, struct parameters *var_pars)
+void variogram_restricts(struct int_par *xD, struct points *pnts,
+ struct parameters *var_pars)
{
- struct write *report = &xD->report;
-
- double *min, *max; // extend
- double dr[3]; // coordinate differences
+ struct write *report = &xD->report;
- int dimension; // dimension: hz / vert / aniso
- char type[12];
- variogram_type(var_pars->type, type); // set up: 0 - "hz", 1 - "vert", 2 - "bivar", 3 - "aniso"
+ double *min, *max; // extend
+ double dr[3]; // coordinate differences
- // 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
+ int dimension; // dimension: hz / vert / aniso
+ char type[12];
- dr[0] = *max - *min; // x range
- dr[1] = *(max+1) - *(min+1); // y range
- dr[2] = *(max+2) - *(min+2); // z range
+ variogram_type(var_pars->type, type); // set up: 0 - "hz", 1 - "vert", 2 - "bivar", 3 - "aniso"
- switch (var_pars->type) {
- case 1: // vertical variogram
- var_pars->max_dist = dr[2]; // zmax - zmin
- var_pars->radius = SQUARE(var_pars->max_dist); // anisotropic distance (todo: try also dz)
- break;
- default:
- var_pars->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
- var_pars->max_dist = sqrt(var_pars->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
- break;
- }
-
- if (report->write2file == TRUE) { // report name available:
- write2file_varSetsIntro(var_pars->type, report); // describe properties
- }
+ // 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
- // set up code according to type
- switch (var_pars->type) {
- case 0: // horizontal
- dimension = 12;
- break;
- case 1: // vertical
- dimension = 3;
- break;
- case 3: // anisotropic
- dimension = 0;
- break;
- }
+ dr[0] = *max - *min; // x range
+ dr[1] = *(max + 1) - *(min + 1); // y range
+ dr[2] = *(max + 2) - *(min + 2); // z range
- if (var_pars->type == 2) { // bivariate variogram
- // horizontal direction:
- var_pars->lag = lag_size(12, pnts, var_pars, report); // lag distance
- var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // number of lags
- optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
- var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ switch (var_pars->type) {
+ case 1: // vertical variogram
+ var_pars->max_dist = dr[2]; // zmax - zmin
+ var_pars->radius = SQUARE(var_pars->max_dist); // anisotropic distance (todo: try also dz)
+ break;
+ default:
+ var_pars->radius = radius_hz_diff(dr) / 9.; // horizontal radius (1/9)
+ var_pars->max_dist = sqrt(var_pars->radius); // 1/3 max horizontal dist (Surfer, Golden SW)
+ break;
+ }
- // vertical direction
- var_pars->lag_vert = lag_size(3, pnts, var_pars, report); // lag distance
- var_pars->nLag_vert = lag_number(var_pars->lag_vert, &var_pars->max_dist_vert); // # of lags
- optimize(&var_pars->lag_vert, &var_pars->nLag_vert, var_pars->max_dist_vert);
- var_pars->max_dist_vert = var_pars->nLag_vert * var_pars->lag_vert; // max distance
- }
+ if (report->write2file == TRUE) { // report name available:
+ write2file_varSetsIntro(var_pars->type, report); // describe properties
+ }
- else { // univariate variograms (hz / vert / aniso)
- var_pars->lag = lag_size(dimension, pnts, var_pars, report); // lag size
- var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // # of lags
- optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
- var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
- }
+ // set up code according to type
+ switch (var_pars->type) {
+ case 0: // horizontal
+ dimension = 12;
+ break;
+ case 1: // vertical
+ dimension = 3;
+ break;
+ case 3: // anisotropic
+ dimension = 0;
+ break;
+ }
- if (report->write2file == TRUE) { // report name available:
- write2file_varSets(report, var_pars); // describe properties
- }
+ if (var_pars->type == 2) { // bivariate variogram
+ // horizontal direction:
+ var_pars->lag = lag_size(12, pnts, var_pars, report); // lag distance
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // number of lags
+ optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+
+ // vertical direction
+ var_pars->lag_vert = lag_size(3, pnts, var_pars, report); // lag distance
+ var_pars->nLag_vert = lag_number(var_pars->lag_vert, &var_pars->max_dist_vert); // # of lags
+ optimize(&var_pars->lag_vert, &var_pars->nLag_vert,
+ var_pars->max_dist_vert);
+ var_pars->max_dist_vert = var_pars->nLag_vert * var_pars->lag_vert; // max distance
+ }
+
+ else { // univariate variograms (hz / vert / aniso)
+ var_pars->lag = lag_size(dimension, pnts, var_pars, report); // lag size
+ var_pars->nLag = lag_number(var_pars->lag, &var_pars->max_dist); // # of lags
+ optimize(&var_pars->lag, &var_pars->nLag, var_pars->max_dist);
+ var_pars->max_dist = var_pars->nLag * var_pars->lag; // maximum distance
+ }
+
+ if (report->write2file == TRUE) { // report name available:
+ write2file_varSets(report, var_pars); // describe properties
+ }
}
void geometric_anisotropy(struct int_par *xD, struct points *pnts)
{
- int i;
- double ratio = xD->aniso_ratio;
+ int i;
+ double ratio = xD->aniso_ratio;
- if (ratio <= 0.) { // ratio is negative or zero:
- if (xD->report.name) { // close report file
- fprintf(xD->report.fp, "Error (see standard output). Process killed...");
- fclose(xD->report.fp);
+ if (ratio <= 0.) { // ratio is negative or zero:
+ if (xD->report.name) { // close report file
+ fprintf(xD->report.fp,
+ "Error (see standard output). Process killed...");
+ fclose(xD->report.fp);
+ }
+ G_fatal_error(_("Anisotropy ratio must be greater than zero..."));
}
- G_fatal_error(_("Anisotropy ratio must be greater than zero..."));
- }
- double *r;
- r = &pnts->r[0];
+ double *r;
- struct RTree *R_tree; // spatial index for input
- struct RTree_Rect *rect;
+ r = &pnts->r[0];
- R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
+ struct RTree *R_tree; // spatial index for input
+ struct RTree_Rect *rect;
- for (i=0; i < pnts->n; i++) {
- *(r+2) = ratio * *(r+2);
+ R_tree = RTreeCreateTree(-1, 0, 3); // create 3D spatial index
- rect = RTreeAllocRect(R_tree);
- RTreeSetRect3D(rect, R_tree, *r, *r, *(r+1), *(r+1), *(r+2), *(r+2));
- RTreeInsertRect(rect, i+1, R_tree);
+ for (i = 0; i < pnts->n; i++) {
+ *(r + 2) = ratio * *(r + 2);
- r += 3;
- }
+ rect = RTreeAllocRect(R_tree);
+ RTreeSetRect3D(rect, R_tree, *r, *r, *(r + 1), *(r + 1), *(r + 2),
+ *(r + 2));
+ RTreeInsertRect(rect, i + 1, R_tree);
- pnts->R_tree = R_tree;
+ r += 3;
+ }
+
+ pnts->R_tree = R_tree;
}
// Least Squares Method
-mat_struct *LSM(mat_struct *A, mat_struct *x)
+mat_struct *LSM(mat_struct * A, mat_struct * x)
{
- mat_struct *AT, *ATA, *ATA_Inv, *ATx, *T;
+ mat_struct *AT, *ATA, *ATA_Inv, *ATx, *T;
- /* 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 */
+ /* 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;
+ return T;
}
// set up type of function according to GUI (user)
int set_function(char *variogram, struct write *report)
{
- int function;
- if (strcmp(variogram, "linear") == 0) {
- function = 0;
- }
- else if (strcmp(variogram, "parabolic") == 0) {
- function = 1;
- }
- else if (strcmp(variogram, "exponential") == 0) {
- function = 2;
- }
- else if (strcmp(variogram, "spherical") == 0) {
- function = 3;
- }
- else if (strcmp(variogram, "gaussian") == 0) {
- function = 4;
- }
- else if (strcmp(variogram, "bivariate") == 0) {
- function = 5;
- }
- else {
- if (report->write2file == TRUE) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
+ int function;
+
+ if (strcmp(variogram, "linear") == 0) {
+ function = 0;
}
- G_fatal_error(_("Set up correct name of variogram function..."));
- }
-
- return function;
+ else if (strcmp(variogram, "parabolic") == 0) {
+ function = 1;
+ }
+ else if (strcmp(variogram, "exponential") == 0) {
+ function = 2;
+ }
+ else if (strcmp(variogram, "spherical") == 0) {
+ function = 3;
+ }
+ else if (strcmp(variogram, "gaussian") == 0) {
+ function = 4;
+ }
+ else if (strcmp(variogram, "bivariate") == 0) {
+ function = 5;
+ }
+ else {
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Set up correct name of variogram function..."));
+ }
+
+ return function;
}
// set up terminal and format of output plot
void set_gnuplot(char *fileformat, struct parameters *var_pars)
{
- if (strcmp(fileformat,"cdr") == 0) {
- strcpy(var_pars->term, "corel");
- strcpy(var_pars->ext, "cdr");
- }
- if (strcmp(fileformat,"dxf") == 0) {
- strcpy(var_pars->term, "dxf");
- strcpy(var_pars->ext, "dxf");
- }
- if (strcmp(fileformat,"eps") == 0) {
- strcpy(var_pars->term, "postscript");
- strcpy(var_pars->ext, "eps");
- }
- if (strcmp(fileformat,"pdf") == 0) {
- strcpy(var_pars->term, "pdfcairo");
- strcpy(var_pars->ext, "pdf");
- }
- if (strcmp(fileformat,"png") == 0) {
- strcpy(var_pars->term, "png");
- strcpy(var_pars->ext, "png");
- }
- if (strcmp(fileformat,"svg") == 0) {
- strcpy(var_pars->term, "svg");
- strcpy(var_pars->ext, "svg");
- }
+ if (strcmp(fileformat, "cdr") == 0) {
+ strcpy(var_pars->term, "corel");
+ strcpy(var_pars->ext, "cdr");
+ }
+ if (strcmp(fileformat, "dxf") == 0) {
+ strcpy(var_pars->term, "dxf");
+ strcpy(var_pars->ext, "dxf");
+ }
+ if (strcmp(fileformat, "eps") == 0) {
+ strcpy(var_pars->term, "postscript");
+ strcpy(var_pars->ext, "eps");
+ }
+ if (strcmp(fileformat, "pdf") == 0) {
+ strcpy(var_pars->term, "pdfcairo");
+ strcpy(var_pars->ext, "pdf");
+ }
+ if (strcmp(fileformat, "png") == 0) {
+ strcpy(var_pars->term, "png");
+ strcpy(var_pars->ext, "png");
+ }
+ if (strcmp(fileformat, "svg") == 0) {
+ strcpy(var_pars->term, "svg");
+ strcpy(var_pars->ext, "svg");
+ }
}
// plot experimental variogram
-void plot_experimental_variogram(struct int_par *xD, struct parameters *var_par)
+void plot_experimental_variogram(struct int_par *xD,
+ struct parameters *var_par)
{
- int bivar = var_par->type == 2 ? TRUE : FALSE;
-
- int i, j; // indices
- int nr, nc; // # of rows, cols
- double *h; // pointer to horizontal or anisotropic bins
- double *vert;
- mat_struct *gamma_M = var_par->gamma; // pointer to gamma matrix
- FILE *gp; // pointer to file
+ int bivar = var_par->type == 2 ? TRUE : FALSE;
- nr = gamma_M->rows; // # of rows of gamma matrix
- nc = gamma_M->cols; // # of cols of gamma matrix
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert;
+ mat_struct *gamma_M = var_par->gamma; // pointer to gamma matrix
+ FILE *gp; // pointer to file
- double *gamma;
- gamma = &gamma_M->vals[0]; // values of gamma matrix
-
- gp = fopen("dataE.dat","w"); // open file to write experimental variogram
- if (access("dataE.dat", W_OK) < 0) {
- G_fatal_error(_("Something went wrong opening tmp file..."));
- }
+ nr = gamma_M->rows; // # of rows of gamma matrix
+ nc = gamma_M->cols; // # of cols of gamma matrix
- /* 3D experimental variogram */
- if (bivar == TRUE) {
- vert = &var_par->vertical.h[0];
- for (i=0; i < nc+1; i++) { // for each row (nZ)
- h = &var_par->horizontal.h[0];
- for (j=0; j < nr+1; j++) { // for each col (nL)
- if (i == 0 && j == 0) {
- fprintf(gp, "%d ", nr); // write to file count of columns
- }
- else if (i == 0 && j != 0) { // 1st row
- fprintf(gp, "%f", *h); // write h to 1st row of the file
- h++;
- if (j < nr) {
- fprintf(gp, " "); // spaces between elements
- }
- } // end 1st row setting
+ double *gamma;
- else {
- if (j == 0 && i != 0) // 1st column
- fprintf(gp, "%f ", *vert);// write vert to 1st column of the file
- else {
- if (isnan(*gamma)) {
- fprintf(gp, "NaN"); // write other elements
- }
- else {
- fprintf(gp, "%f", *gamma); // write other elements
- }
-
- if (j != nr) {
- fprintf(gp, " "); // spaces between elements
- }
- gamma++;
- }
- } // end columns settings
- } // end j
- fprintf(gp, "\n");
- if (i!=0) { // vert must not increase in the 1st loop
- vert++;
- }
- } // end i
- } // end 3D
+ gamma = &gamma_M->vals[0]; // values of gamma matrix
- /* 2D experimental variogram */
- else {
- h = &var_par->h[0];
- for (i=0; i < nr; ++i) {
- if (isnan(*gamma)) {
- fprintf(gp, "%f NaN\n", *h);
- }
- else {
- fprintf(gp, "%f %f\n", *h, *gamma);
- }
- h++;
- gamma++;
+ gp = fopen("dataE.dat", "w"); // open file to write experimental variogram
+ if (access("dataE.dat", W_OK) < 0) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
}
- }
- fclose(gp);
- gp = popen(GNUPLOT, "w"); /* open file to create plots */
- if (gp == NULL)
- G_message("Unable to plot variogram");
+ /* 3D experimental variogram */
+ if (bivar == TRUE) {
+ vert = &var_par->vertical.h[0];
+ for (i = 0; i < nc + 1; i++) { // for each row (nZ)
+ h = &var_par->horizontal.h[0];
+ for (j = 0; j < nr + 1; j++) { // for each col (nL)
+ if (i == 0 && j == 0) {
+ fprintf(gp, "%d ", nr); // write to file count of columns
+ }
+ else if (i == 0 && j != 0) { // 1st row
+ fprintf(gp, "%f", *h); // write h to 1st row of the file
+ h++;
+ if (j < nr) {
+ fprintf(gp, " "); // spaces between elements
+ }
+ } // end 1st row setting
- 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 {
+ if (j == 0 && i != 0) // 1st column
+ fprintf(gp, "%f ", *vert); // write vert to 1st column of the file
+ else {
+ if (isnan(*gamma)) {
+ fprintf(gp, "NaN"); // write other elements
+ }
+ else {
+ fprintf(gp, "%f", *gamma); // write other elements
+ }
- else { // univariate variogram
- char dim[6];
- 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);
+ if (j != nr) {
+ fprintf(gp, " "); // spaces between elements
+ }
+ gamma++;
+ }
+ } // end columns settings
+ } // end j
+ fprintf(gp, "\n");
+ if (i != 0) { // vert must not increase in the 1st loop
+ vert++;
+ }
+ } // end i
+ } // end 3D
+
+ /* 2D experimental variogram */
+ else {
+ h = &var_par->h[0];
+ for (i = 0; i < nr; ++i) {
+ if (isnan(*gamma)) {
+ fprintf(gp, "%f NaN\n", *h);
+ }
+ else {
+ fprintf(gp, "%f %f\n", *h, *gamma);
+ }
+ h++;
+ gamma++;
+ }
}
- else { // 2D
- fprintf(gp, "set title \"Experimental variogram (2D) of <%s>\"\n", var_par->name);
+ 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");
}
- 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);
+ else { // univariate variogram
+ char dim[6];
+
+ if (xD->i3 == TRUE) { // 3D
+ if (var_par->type == 0) {
+ strcpy(dim, "hz");
+ }
+ if (var_par->type == 1) {
+ strcpy(dim, "vert");
+ }
+ if (var_par->type == 3) {
+ strcpy(dim, "aniso");
+ }
+ fprintf(gp,
+ "set title \"Experimental variogram (3D %s) of <%s>\"\n",
+ dim, var_par->name);
+ }
+ else { // 2D
+ fprintf(gp, "set title \"Experimental variogram (2D) of <%s>\"\n",
+ var_par->name);
+ }
+
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp,
+ "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5\n");
+ }
+ fclose(gp);
}
// plot experimental and theoretical variogram
void plot_var(int i3, int bivar, struct parameters *var_pars)
{
- int function;
- double nugget, nugget_h, nugget_v;
- double sill;
- double h_range;
- double *T;
+ int function;
+ double nugget, nugget_h, nugget_v;
+ double sill;
+ double h_range;
+ double *T;
- // setup theoretical variogram parameters
- if (bivar == FALSE) { // univariate variogram
- function = var_pars->function;
- if (function > 1) { // nonlinear and not parabolic variogram
- nugget = var_pars->nugget;
- sill = var_pars->sill - nugget;
- h_range = var_pars->h_range;
+ // setup theoretical variogram parameters
+ if (bivar == FALSE) { // univariate variogram
+ function = var_pars->function;
+ if (function > 1) { // nonlinear and not parabolic variogram
+ nugget = var_pars->nugget;
+ sill = var_pars->sill - nugget;
+ h_range = var_pars->h_range;
+ }
}
- }
- if (var_pars->function < 2) { // linear or parabolic variogram
- T = &var_pars->T->vals[0];
- }
+ if (var_pars->function < 2) { // linear or parabolic variogram
+ T = &var_pars->T->vals[0];
+ }
- int i, j; // indices
- int nr, nc; // # of rows, cols
- int type = var_pars->type;
- double *h; // pointer to horizontal or anisotropic bins
- double *vert; // pointer to vertical bins
- double h_ratio;
- double g_teor;
- double *gamma_var; // pointer to gamma matrix
- FILE *gp; // pointer to file
+ int i, j; // indices
+ int nr, nc; // # of rows, cols
+ int type = var_pars->type;
+ double *h; // pointer to horizontal or anisotropic bins
+ double *vert; // pointer to vertical bins
+ double h_ratio;
+ double g_teor;
+ double *gamma_var; // pointer to gamma matrix
+ FILE *gp; // pointer to file
- nr = var_pars->gamma->rows; // # of rows of gamma matrix
- nc = var_pars->gamma->cols; // # of cols of gamma matrix
+ nr = var_pars->gamma->rows; // # of rows of gamma matrix
+ nc = var_pars->gamma->cols; // # of cols of gamma matrix
- gamma_var = &var_pars->gamma->vals[0]; // values of gamma matrix
+ gamma_var = &var_pars->gamma->vals[0]; // values of gamma matrix
- if (type != 2) { // univariate variogram
- gp = fopen("dataE.dat","w"); // open file to write experimental variogram
- if (access("dataE.dat", W_OK) < 0) {
- G_fatal_error(_("Something went wrong opening tmp file..."));
- }
- if (gp == NULL) {
- G_fatal_error(_("Error writing file"));
+ if (type != 2) { // univariate variogram
+ gp = fopen("dataE.dat", "w"); // open file to write experimental variogram
+ if (access("dataE.dat", W_OK) < 0) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
+ }
+ if (gp == NULL) {
+ G_fatal_error(_("Error writing file"));
+ }
+
+ h = &var_pars->h[0];
+ for (i = 0; i < nr; i++) {
+ if (isnan(*gamma_var)) {
+ fprintf(gp, "%f NaN\n", *h); // write other elements
+ }
+ else {
+ fprintf(gp, "%f %f\n", *h, *gamma_var);
+ }
+ h++;
+ gamma_var++;
+ }
}
- h = &var_pars->h[0];
- for (i=0; i < nr; i++) {
- if (isnan(*gamma_var)) {
- fprintf(gp, "%f NaN\n", *h); // write other elements
- }
- else {
- fprintf(gp, "%f %f\n", *h, *gamma_var);
- }
- h++;
- gamma_var++;
- }
- }
+ else { // bivariate variogram
+ gp = fopen("dataE.dat", "r"); // open file to read experimental variogram
+ if (access("dataE.dat", F_OK) < 0) {
+ G_fatal_error(_("You have probably deleted dataE.dat - process middle phase again, please."));
+ }
+ }
- else { // bivariate variogram
- gp = fopen("dataE.dat","r"); // open file to read experimental variogram
- if (access("dataE.dat", F_OK) < 0) {
- G_fatal_error(_("You have probably deleted dataE.dat - process middle phase again, please."));
+ if (fclose(gp) != 0) {
+ G_fatal_error(_("Error closing file..."));
}
- }
-
- if (fclose(gp) != 0) {
- G_fatal_error(_("Error closing file..."));
- }
- /* Theoretical variogram */
- gp = fopen("dataT.dat","w"); // open file to write theoretical variogram
- if (access("dataT.dat", W_OK) < 0) {
- G_fatal_error(_("Something went wrong opening tmp file..."));
- }
+ /* Theoretical variogram */
+ gp = fopen("dataT.dat", "w"); // open file to write theoretical variogram
+ if (access("dataT.dat", W_OK) < 0) {
+ G_fatal_error(_("Something went wrong opening tmp file..."));
+ }
- double dd;
- double hh;
- double hv[2];
- h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
-
- // Compute theoretical variogram
- switch (bivar) {
- case 0: // Univariate variogram:
- for (i=0; i <= nr; i++) {
- hh = i==0 ? 0. : *h;
+ double dd;
+ double hh;
+ double hv[2];
- switch (function) {
- case 0: // linear
- dd = *T * hh; // todo: add nugget effect
- 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);
+ h = type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
- if (i > 0) {
- h++;
- }
- } // end i for cycle
- break;
+ // Compute theoretical variogram
+ switch (bivar) {
+ case 0: // Univariate variogram:
+ for (i = 0; i <= nr; i++) {
+ hh = i == 0 ? 0. : *h;
- case 1: // bivariate (3D)
- // Coefficients of theoretical variogram
- T = &var_pars->T->vals[0]; // values
- vert = &var_pars->vertical.h[0];
- for (i=0; i < nc+1; i++) { // for each row
- h = &var_pars->horizontal.h[0];
- for (j=0; j < nr+1; j++) { // for each col
- if (i == 0 && j == 0) { // the 0-th element...
- fprintf(gp, "%d ", nr); // ... is number of columns
- }
- else if (i == 0 && j != 0) { // for the other elements in 1st row...
- fprintf(gp, "%f", *h); // ... write h
- if (j < nr) {
- fprintf(gp, " "); // ... write a space
- }
- h++; // go to next h
- }
- else { // other rows
- if (j == 0 && i != 0) {// write vert to 1st column
- fprintf(gp, "%f ", *vert);
- }
- else { // fill the other columns with g_teor
- hv[0] = *h;
- hv[1] = *vert;
- g_teor = variogram_fction(var_pars, hv);
- fprintf(gp, "%f", g_teor);
- if (j != nr) {
- fprintf(gp, " ");
- }
- h++;
- } // end else: fill the other columns
- } // end fill the other rows
- } //end j
- if (i != 0 && j != 0) {
- vert++;
- }
- fprintf(gp, "\n");
- } //end i
- break;
- }
+ switch (function) {
+ case 0: // linear
+ dd = *T * hh; // todo: add nugget effect
+ 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 (fclose(gp) == EOF) {
- G_fatal_error(_("Error closing file..."));
- }
+ if (i > 0) {
+ h++;
+ }
+ } // end i for cycle
+ break;
- gp = popen(GNUPLOT, "w"); // open file to create plots
- if (gp == NULL) {
- G_message(_("Unable to plot variogram"));
- }
-
- if (strcmp(var_pars->term," ") != 0) {
- fprintf(gp, "set terminal %s size 750,450\n", var_pars->term);
-
- switch (var_pars->type) {
- case 0: // horizontal component of bivariate variogram
- if (var_pars->function == 5) { // component of 3D bivariate
- fprintf(gp, "set output \"variogram_bivar.%s\"\n", var_pars->ext);
- }
- else { // horizontal variogram
- fprintf(gp, "set output \"variogram_hz.%s\"\n", var_pars->ext);
- }
- break;
- case 1: // vertical variogram
- fprintf(gp, "set output \"variogram_vertical.%s\"\n", var_pars->ext);
- break;
- case 2: // bivariate variogram
- if (bivar == TRUE) {
- fprintf(gp, "set output \"variogram_bivariate.%s\"\n", var_pars->ext);
- }
- break;
- case 3: // anisotropic variogram
- fprintf(gp, "set output \"variogram_final.%s\"\n", var_pars->ext);
- break;
+ case 1: // bivariate (3D)
+ // Coefficients of theoretical variogram
+ T = &var_pars->T->vals[0]; // values
+ vert = &var_pars->vertical.h[0];
+ for (i = 0; i < nc + 1; i++) { // for each row
+ h = &var_pars->horizontal.h[0];
+ for (j = 0; j < nr + 1; j++) { // for each col
+ if (i == 0 && j == 0) { // the 0-th element...
+ fprintf(gp, "%d ", nr); // ... is number of columns
+ }
+ else if (i == 0 && j != 0) { // for the other elements in 1st row...
+ fprintf(gp, "%f", *h); // ... write h
+ if (j < nr) {
+ fprintf(gp, " "); // ... write a space
+ }
+ h++; // go to next h
+ }
+ else { // other rows
+ if (j == 0 && i != 0) { // write vert to 1st column
+ fprintf(gp, "%f ", *vert);
+ }
+ else { // fill the other columns with g_teor
+ hv[0] = *h;
+ hv[1] = *vert;
+ g_teor = variogram_fction(var_pars, hv);
+ fprintf(gp, "%f", g_teor);
+ if (j != nr) {
+ fprintf(gp, " ");
+ }
+ h++;
+ } // end else: fill the other columns
+ } // end fill the other rows
+ } //end j
+ if (i != 0 && j != 0) {
+ vert++;
+ }
+ fprintf(gp, "\n");
+ } //end i
+ break;
}
- }
- else {
- G_warning("\nVariogram output> If you wish some special file format, please set it at the beginning.\n");
- fprintf(gp, "set terminal wxt %d size 800,450\n", var_pars->type);
- }
- if (bivar == TRUE) { // bivariate variogram
- fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_pars->name);
-
- fprintf(gp, "set xlabel \"h interval No\"\n");
- fprintf(gp, "set ylabel \"dz interval No\"\n");
- fprintf(gp, "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
- fprintf(gp, "set key below\n");
- fprintf(gp, "set key box\n");
- fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
- fprintf(gp, "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\", 'dataT.dat' every ::1:1 matrix with lines title \"theoretical variogram\" palette\n");
- }
+ if (fclose(gp) == EOF) {
+ G_fatal_error(_("Error closing file..."));
+ }
- else { // univariate variogram
- if (i3 == TRUE) { // 3D variogram
- if (var_pars->type == 0) { // horizontal variogram
- if (i3 == TRUE) {
- fprintf(gp, "set title \"Experimental and theoretical variogram (3D hz) of <%s>\"\n", var_pars->name);
- }
- else {
- fprintf(gp, "set title \"Experimental and theoretical variogram (2D hz) of <%s>\"\n", var_pars->name);
- }
- }
- else if (var_pars->type == 1) { // vertical variogram
- fprintf(gp, "set title \"Experimental and theoretical variogram (3D vert) of <%s>\"\n", var_pars->name);
- }
- else if (var_pars->type == 3) {// anisotropic variogram
- fprintf(gp, "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n", var_pars->name);
- }
+ gp = popen(GNUPLOT, "w"); // open file to create plots
+ if (gp == NULL) {
+ G_message(_("Unable to plot variogram"));
}
- else { // 2D
- fprintf(gp, "set title \"Experimental and theoretical variogram (2D) of <%s>\"\n", var_pars->name);
+ if (strcmp(var_pars->term, " ") != 0) {
+ fprintf(gp, "set terminal %s size 750,450\n", var_pars->term);
+
+ switch (var_pars->type) {
+ case 0: // horizontal component of bivariate variogram
+ if (var_pars->function == 5) { // component of 3D bivariate
+ fprintf(gp, "set output \"variogram_bivar.%s\"\n",
+ var_pars->ext);
+ }
+ else { // horizontal variogram
+ fprintf(gp, "set output \"variogram_hz.%s\"\n",
+ var_pars->ext);
+ }
+ break;
+ case 1: // vertical variogram
+ fprintf(gp, "set output \"variogram_vertical.%s\"\n",
+ var_pars->ext);
+ break;
+ case 2: // bivariate variogram
+ if (bivar == TRUE) {
+ fprintf(gp, "set output \"variogram_bivariate.%s\"\n",
+ var_pars->ext);
+ }
+ break;
+ case 3: // anisotropic variogram
+ fprintf(gp, "set output \"variogram_final.%s\"\n", var_pars->ext);
+ break;
+ }
}
+ else {
+ G_warning
+ ("\nVariogram output> If you wish some special file format, please set it at the beginning.\n");
+ fprintf(gp, "set terminal wxt %d size 800,450\n", var_pars->type);
+ }
- switch (var_pars->function) {
- case 0: // linear
- fprintf(gp, "set label \"linear: gamma(h) = %f*h\" at screen 0.30,0.90\n", var_pars->T->vals[0]);
- break;
- case 1: // parabolic
- fprintf(gp, "set label \"parabolic: gamma(h) = %f*h^2\" at screen 0.25,0.90\n", var_pars->T->vals[0]);
- break;
- case 2: // exponential
- fprintf(gp, "set rmargin 5\n");
- fprintf(gp, "set label \"exponential: gamma(h) = %f + %f * (1 - exp(-3*h / %f))\" at screen 0.10,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range);
- break;
- case 3: // spherical
- fprintf(gp, "set rmargin 5\n");
- fprintf(gp, "set label \"spherical: gamma(h) = %f+%f*(1.5*h/%f-0.5*(h/%f)^3)\" at screen 0.05,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range, var_pars->h_range);
- break;
- case 4: // gaussian
- fprintf(gp, "set label \"gaussian: gamma(h) = %f + %f * (1 - exp(-3*(h / %f)^2))\" at screen 0.10,0.90\n", var_pars->nugget, var_pars->sill - var_pars->nugget, var_pars->h_range);
- break;
+ if (bivar == TRUE) { // bivariate variogram
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n",
+ var_pars->name);
+
+ fprintf(gp, "set xlabel \"h interval No\"\n");
+ fprintf(gp, "set ylabel \"dz interval No\"\n");
+ fprintf(gp,
+ "set zlabel \"gamma(h,dz) [units^2]\" rotate by 90 left\n");
+ fprintf(gp, "set key below\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp, "set dgrid3d %d,%d\n", nc, nr);
+ fprintf(gp,
+ "splot 'dataE.dat' every ::1:1 matrix title \"experimental variogram\", 'dataT.dat' every ::1:1 matrix with lines title \"theoretical variogram\" palette\n");
}
-
- 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");
+ else { // univariate variogram
+ if (i3 == TRUE) { // 3D variogram
+ if (var_pars->type == 0) { // horizontal variogram
+ if (i3 == TRUE) {
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (3D hz) of <%s>\"\n",
+ var_pars->name);
+ }
+ else {
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (2D hz) of <%s>\"\n",
+ var_pars->name);
+ }
+ }
+ else if (var_pars->type == 1) { // vertical variogram
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (3D vert) of <%s>\"\n",
+ var_pars->name);
+ }
+ else if (var_pars->type == 3) { // anisotropic variogram
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (3D) of <%s>\"\n",
+ var_pars->name);
+ }
+ }
+
+ else { // 2D
+ fprintf(gp,
+ "set title \"Experimental and theoretical variogram (2D) of <%s>\"\n",
+ var_pars->name);
+ }
+
+ switch (var_pars->function) {
+ case 0: // linear
+ fprintf(gp,
+ "set label \"linear: gamma(h) = %f*h\" at screen 0.30,0.90\n",
+ var_pars->T->vals[0]);
+ break;
+ case 1: // parabolic
+ fprintf(gp,
+ "set label \"parabolic: gamma(h) = %f*h^2\" at screen 0.25,0.90\n",
+ var_pars->T->vals[0]);
+ break;
+ case 2: // exponential
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp,
+ "set label \"exponential: gamma(h) = %f + %f * (1 - exp(-3*h / %f))\" at screen 0.10,0.90\n",
+ var_pars->nugget, var_pars->sill - var_pars->nugget,
+ var_pars->h_range);
+ break;
+ case 3: // spherical
+ fprintf(gp, "set rmargin 5\n");
+ fprintf(gp,
+ "set label \"spherical: gamma(h) = %f+%f*(1.5*h/%f-0.5*(h/%f)^3)\" at screen 0.05,0.90\n",
+ var_pars->nugget, var_pars->sill - var_pars->nugget,
+ var_pars->h_range, var_pars->h_range);
+ break;
+ case 4: // gaussian
+ fprintf(gp,
+ "set label \"gaussian: gamma(h) = %f + %f * (1 - exp(-3*(h / %f)^2))\" at screen 0.10,0.90\n",
+ var_pars->nugget, var_pars->sill - var_pars->nugget,
+ var_pars->h_range);
+ break;
+ }
+
+ fprintf(gp, "set xlabel \"h [m]\"\n");
+ fprintf(gp, "set ylabel \"gamma(h) [units^2]\"\n");
+ fprintf(gp, "set key bottom right\n");
+ fprintf(gp, "set key box\n");
+ fprintf(gp,
+ "plot 'dataE.dat' using 1:2 title \"experimental variogram\" pointtype 5, 'dataT.dat' using 1:2 title \"theoretical variogram\" with linespoints\n");
+ }
+ fclose(gp);
+
+ remove("dataE.dat");
+ remove("dataT.dat");
}
Modified: grass-addons/grass7/vector/v.kriging/utils_kriging.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_kriging.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/utils_kriging.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -2,675 +2,713 @@
void LMS_variogram(struct parameters *var_par, struct write *report)
{
- int nZ, nL; // # of gamma matrix rows and columns
- int nr, nc; // # of plan matrix rows and columns
- int i, j; // indices
- double *h, *vert, *gamma;
- mat_struct *gR;
-
- nL = var_par->gamma->rows; // # of cols (horizontal bins)
- nZ = var_par->gamma->cols; // # of rows (vertical bins)
- gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
-
- // Test length of design matrix
- nr = 0; // # of rows of design matrix A - depend on # of non-nan gammas
- for (i = 0; i < nZ; i++) {
- for (j = 0; j < nL; j++) {
- if (!isnan(*gamma)) { // todo: upgrade: !isnan(*c) L434 too
- nr++;
- }
- gamma++;
- } // end j
- } // end i
-
- // # of columns of design matrix A
- nc = var_par->function == 5 ? 3 : 1;
+ int nZ, nL; // # of gamma matrix rows and columns
+ int nr, nc; // # of plan matrix rows and columns
+ int i, j; // indices
+ double *h, *vert, *gamma;
+ mat_struct *gR;
- 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
+ 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
- if (var_par->function == 5) { // in case of bivariate variogram
- vert = &var_par->vertical.h[0]; // vertical direction
- }
+ // Test length of design matrix
+ nr = 0; // # of rows of design matrix A - depend on # of non-nan gammas
+ for (i = 0; i < nZ; i++) {
+ for (j = 0; j < nL; j++) {
+ if (!isnan(*gamma)) { // todo: upgrade: !isnan(*c) L434 too
+ nr++;
+ }
+ gamma++;
+ } // end j
+ } // end i
- for (i = 0; i < nZ; i++) {
- h = var_par->function == 5 ? &var_par->horizontal.h[0] : &var_par->h[0];
- for (j = 0; j < nL; j++) {
- if (!isnan(*gamma)) { // write to matrix - each valid element of gamma
- switch (var_par->function) { // function of theoretical variogram
- case 0: // linear variogram
- G_matrix_set_element(var_par->A, nr, 0, *h);
- //G_matrix_set_element(var_par->A, nr, 1, 1.);
- break;
- case 1: // parabolic variogram
- G_matrix_set_element(var_par->A, nr, 0, SQUARE(*h));
- G_matrix_set_element(var_par->A, nr, 1, 1.);
- break;
- 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.);
- break;
- } // end switch variogram fuction
- G_matrix_set_element(gR, nr, 0, *gamma);
- nr++; // length of vector of valid elements (not null)
- } // end test if !isnan(*gamma)
- h++;
- gamma++;
- } // end j
- if (var_par->function == 5) {
- vert++;
- }
- } // end i
+ // # of columns of design matrix A
+ nc = var_par->function == 5 ? 3 : 1;
- // 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);
+ var_par->A = G_matrix_init(nr, nc, nr); // initialise design matrix
+ gR = G_matrix_init(nr, 1, nr); // initialise vector of observations
+
+ // Set design matrix A
+ nr = 0; // index of matrix row
+ gamma = &var_par->gamma->vals[0]; // pointer to experimental variogram
+
+ if (var_par->function == 5) { // in case of bivariate variogram
+ vert = &var_par->vertical.h[0]; // vertical direction
}
- G_fatal_error("Unable to compute 3D theoretical variogram...");
- } // error
- // constant raster
- if (var_par->T->vals[0] == 0. && var_par->T->vals[1] == 0.) {
- var_par->const_val = 1;
- G_message(_("Input values to be interpolated are constant."));
- } // todo: constant raster for exponential etc.
+ for (i = 0; i < nZ; i++) {
+ h = var_par->function ==
+ 5 ? &var_par->horizontal.h[0] : &var_par->h[0];
+ for (j = 0; j < nL; j++) {
+ if (!isnan(*gamma)) { // write to matrix - each valid element of gamma
+ switch (var_par->function) { // function of theoretical variogram
+ case 0: // linear variogram
+ G_matrix_set_element(var_par->A, nr, 0, *h);
+ //G_matrix_set_element(var_par->A, nr, 1, 1.);
+ break;
+ case 1: // parabolic variogram
+ G_matrix_set_element(var_par->A, nr, 0, SQUARE(*h));
+ G_matrix_set_element(var_par->A, nr, 1, 1.);
+ break;
+ 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.);
+ break;
+ } // end switch variogram fuction
+ G_matrix_set_element(gR, nr, 0, *gamma);
+ nr++; // length of vector of valid elements (not null)
+ } // end test if !isnan(*gamma)
+ h++;
+ gamma++;
+ } // end j
+ if (var_par->function == 5) {
+ vert++;
+ }
+ } // end i
- // coefficients of theoretical variogram (linear)
- if (report->name) {
- fprintf(report->fp, "Parameters of bivariate variogram:\n");
- fprintf(report->fp, "Function: linear\n\n");
- fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + %f\n", var_par->T->vals[0], var_par->T->vals[1], var_par->T->vals[2]);
- } // end if: report
+ // Estimate theoretical variogram coefficients
+ var_par->T = LSM(var_par->A, gR); // Least Square Method
+
+ // Test of theoretical variogram estimation
+ if (var_par->T->vals == NULL) { // NULL values of theoretical variogram
+ if (report->write2file == TRUE) { // close report file
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error("Unable to compute 3D theoretical variogram...");
+ } // error
+
+ // constant raster
+ if (var_par->T->vals[0] == 0. && var_par->T->vals[1] == 0.) {
+ var_par->const_val = 1;
+ G_message(_("Input values to be interpolated are constant."));
+ } // todo: constant raster for exponential etc.
+
+ // coefficients of theoretical variogram (linear)
+ if (report->name) {
+ fprintf(report->fp, "Parameters of bivariate variogram:\n");
+ fprintf(report->fp, "Function: linear\n\n");
+ fprintf(report->fp, "gamma(h, vert) = %f * h + %f * vert + %f\n",
+ var_par->T->vals[0], var_par->T->vals[1],
+ var_par->T->vals[2]);
+ } // end if: report
}
-double bivar_sill(int direction, mat_struct *gamma)
+double bivar_sill(int direction, mat_struct * gamma)
{
- int i, j, k;
- int n; // # of lags
- int n_gamma = 0; // # of real gammas
- double gamma_i; // each real gamma
- double sum_gamma = 0.; // sum of real gammas
- double sill; // variogram sill
+ int i, j, k;
+ int n; // # of lags
+ int n_gamma = 0; // # of real gammas
+ double gamma_i; // each real gamma
+ double sum_gamma = 0.; // sum of real gammas
+ double sill; // variogram sill
- switch (direction) { // direction is:
- case 12: // horizontal
- n = gamma->rows; // use 1st column of experimental variogram matrix
- break;
- case 3: // vertical
- n = gamma->cols; // use 1st row of experimental variogram matrix
- break;
- }
-
- for (i=0; i<n; i++) {
- switch (direction) { // direction is:
- case 12: // horizontal
- j = i; // row index is increasing
- k = 0; // column index is constant
- break;
- case 3: // vertical
- j = 0; // row index is constant
- k = i; // column index is increasing
- break;
+ switch (direction) { // direction is:
+ case 12: // horizontal
+ n = gamma->rows; // use 1st column of experimental variogram matrix
+ break;
+ case 3: // vertical
+ n = gamma->cols; // use 1st row of experimental variogram matrix
+ break;
}
- gamma_i = G_matrix_get_element(gamma, j, k);
- if (!isnan(gamma_i)) { // gamma is real:
- sum_gamma += gamma_i; // sum all real elements of the matrix
- n_gamma++; // count them
- } // end if
- } // end for
+ for (i = 0; i < n; i++) {
+ switch (direction) { // direction is:
+ case 12: // horizontal
+ j = i; // row index is increasing
+ k = 0; // column index is constant
+ break;
+ case 3: // vertical
+ j = 0; // row index is constant
+ k = i; // column index is increasing
+ break;
+ }
- sill = sum_gamma / n_gamma;
+ gamma_i = G_matrix_get_element(gamma, j, k);
+ if (!isnan(gamma_i)) { // gamma is real:
+ sum_gamma += gamma_i; // sum all real elements of the matrix
+ n_gamma++; // count them
+ } // end if
+ } // end for
- return sill;
+ sill = sum_gamma / n_gamma;
+
+ return sill;
}
// Compute sill
void sill(struct parameters *var_par)
{
- // Local variables
- int type = var_par->type;
- mat_struct *gamma = var_par->gamma;
+ // Local variables
+ int type = var_par->type;
+ mat_struct *gamma = var_par->gamma;
- char var_type[12];
+ char var_type[12];
- switch (type) {
- case 2: // bivariate variogram
- var_par->horizontal.sill = bivar_sill(12, gamma);
- var_par->vertical.sill = bivar_sill(3, gamma);
- break;
+ switch (type) {
+ case 2: // bivariate variogram
+ var_par->horizontal.sill = bivar_sill(12, gamma);
+ var_par->vertical.sill = bivar_sill(3, gamma);
+ break;
- default: // hz / vert / aniso variogram
- var_par->sill = var_par->gamma_sum / var_par->gamma_n; // mean of gamma elements
-
- variogram_type(var_par->type, var_type); // describe code by string
- G_message(_("Default sill of %s variogram: %f"), var_type, var_par->sill);
- break;
- } // end switch: variogram type
+ default: // hz / vert / aniso variogram
+ var_par->sill = var_par->gamma_sum / var_par->gamma_n; // mean of gamma elements
+
+ variogram_type(var_par->type, var_type); // describe code by string
+ G_message(_("Default sill of %s variogram: %f"), var_type,
+ var_par->sill);
+ break;
+ } // end switch: variogram type
}
// compare sills
-void sill_compare(struct int_par *xD, struct flgs *flg, struct var_par *var_par, struct points *pnts)
+void sill_compare(struct int_par *xD, struct flgs *flg,
+ struct var_par *var_par, struct points *pnts)
{
- // local variables
- double sill_hz = var_par->hz.sill; // sill in horizontal direction
- double sill_vert = var_par->vert.sill; // sill in vertical direction
+ // local variables
+ double sill_hz = var_par->hz.sill; // sill in horizontal direction
+ double sill_vert = var_par->vert.sill; // sill in vertical direction
- double diff_sill_05, diff_sill;
+ double diff_sill_05, diff_sill;
- diff_sill_05 = sill_hz > sill_vert ? 0.05 * sill_hz : 0.05 * sill_vert; // critical value as 5% from bigger sill
- diff_sill = fabsf(sill_hz - sill_vert); // difference between the sills
-
- if (xD->bivar == TRUE || (!flg->univariate->answer && diff_sill > diff_sill_05)) { // zonal anisotropy
- var_par->fin.type = 2; // code for bivariate variogram
- var_par->fin.td = var_par->hz.td; // azimuth limit
- var_par->fin.horizontal.max_dist = var_par->hz.max_dist; // maximum dist in hz direction
- var_par->fin.horizontal.nLag = var_par->hz.nLag;
- var_par->fin.horizontal.lag = var_par->hz.lag;
+ diff_sill_05 = sill_hz > sill_vert ? 0.05 * sill_hz : 0.05 * sill_vert; // critical value as 5% from bigger sill
+ diff_sill = fabsf(sill_hz - sill_vert); // difference between the sills
- var_par->fin.vertical.max_dist = var_par->vert.max_dist; // maximum dist in vert direction
- var_par->fin.vertical.nLag = var_par->vert.nLag;
- var_par->fin.vertical.lag = var_par->vert.lag;
- }
-
- else if (xD->univar == TRUE || (!flg->bivariate->answer && diff_sill <= diff_sill_05)) { // geometric anisotropy
- var_par->fin.type = 3; // variogram type: anisotropic
- var_par->fin.max_dist = var_par->hz.max_dist; // maximum distance: hz
- var_par->fin.td = var_par->hz.td; // azimuth: hz
- xD->aniso_ratio = var_par->hz.h_range / var_par->vert.h_range; // anisotropic ratio
- geometric_anisotropy(xD, pnts); // exaggerate z coords and build a new spatial index
- }
+ if (xD->bivar == TRUE || (!flg->univariate->answer && diff_sill > diff_sill_05)) { // zonal anisotropy
+ var_par->fin.type = 2; // code for bivariate variogram
+ var_par->fin.td = var_par->hz.td; // azimuth limit
+ var_par->fin.horizontal.max_dist = var_par->hz.max_dist; // maximum dist in hz direction
+ var_par->fin.horizontal.nLag = var_par->hz.nLag;
+ var_par->fin.horizontal.lag = var_par->hz.lag;
- return 0;
+ var_par->fin.vertical.max_dist = var_par->vert.max_dist; // maximum dist in vert direction
+ var_par->fin.vertical.nLag = var_par->vert.nLag;
+ var_par->fin.vertical.lag = var_par->vert.lag;
+ }
+
+ else if (xD->univar == TRUE || (!flg->bivariate->answer && diff_sill <= diff_sill_05)) { // geometric anisotropy
+ var_par->fin.type = 3; // variogram type: anisotropic
+ var_par->fin.max_dist = var_par->hz.max_dist; // maximum distance: hz
+ var_par->fin.td = var_par->hz.td; // azimuth: hz
+ xD->aniso_ratio = var_par->hz.h_range / var_par->vert.h_range; // anisotropic ratio
+ geometric_anisotropy(xD, pnts); // exaggerate z coords and build a new spatial index
+ }
+
+ return 0;
}
// formulation of variogram functions
-double variogram_fction(struct parameters *var_par, double *dr)
+double variogram_fction(struct parameters *var_par, double *dr)
{
- // Local variables
- int i;
- int variogram = var_par->function; // function of theoretical variogram
- int type = var_par->type; // hz / vert / bivar / aniso
- double nugget;
- double part_sill;
- double h_range;
- double *T;
+ // Local variables
+ int i;
+ int variogram = var_par->function; // function of theoretical variogram
+ int type = var_par->type; // hz / vert / bivar / aniso
+ double nugget;
+ double part_sill;
+ double h_range;
+ double *T;
- if (type == 2 || var_par->function == 0) { // bivariate variogram:
- T = &var_par->T->vals[0]; // coefficients of theoretical variogram
- }
-
- double radius; // square distance of the point couple
- double h;
- double vert;
- double teor_var, result = 0.; // estimated value of the theoretical variogram
+ if (type == 2 || var_par->function == 0) { // bivariate variogram:
+ T = &var_par->T->vals[0]; // coefficients of theoretical variogram
+ }
- int n_cycles = (type == 2 && var_par->function != 5) ? 2 : 1; // # of cycles (bivariate (not linear) or univariate)
+ double radius; // square distance of the point couple
+ double h;
+ double vert;
+ double teor_var, result = 0.; // estimated value of the theoretical variogram
- for (i=0; i < n_cycles; i++) {
- if (n_cycles == 2) { // bivariate variogram
- variogram = i==0 ? var_par->horizontal.function : var_par->vertical.function;
- h = i==0 ? dr[0] : dr[1];
- radius = SQUARE(h);
- nugget = i==0 ? var_par->horizontal.nugget : var_par->vertical.nugget;
- part_sill = i==0 ? var_par->horizontal.sill - nugget : var_par->vertical.sill - nugget;
- h_range = i==0 ? var_par->horizontal.h_range : var_par->vertical.h_range;
+ int n_cycles = (type == 2 && var_par->function != 5) ? 2 : 1; // # of cycles (bivariate (not linear) or univariate)
+
+ for (i = 0; i < n_cycles; i++) {
+ if (n_cycles == 2) { // bivariate variogram
+ variogram =
+ i ==
+ 0 ? var_par->horizontal.function : var_par->vertical.function;
+ h = i == 0 ? dr[0] : dr[1];
+ radius = SQUARE(h);
+ nugget =
+ i ==
+ 0 ? var_par->horizontal.nugget : var_par->vertical.nugget;
+ part_sill =
+ i ==
+ 0 ? var_par->horizontal.sill -
+ nugget : var_par->vertical.sill - nugget;
+ h_range =
+ i ==
+ 0 ? var_par->horizontal.h_range : var_par->vertical.h_range;
+ }
+ else { // univariate variogram or linear bivariate
+ variogram = var_par->function;
+ if (variogram == 5) {
+ h = dr[0];
+ vert = dr[1];
+ }
+ else { // univariate variogram
+ radius = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]); // weighted dz
+ h = sqrt(radius);
+
+ nugget = var_par->nugget;
+ part_sill = var_par->sill - nugget;
+ h_range = var_par->h_range;
+ }
+ }
+
+ switch (variogram) {
+ case 0: // linear function
+ teor_var = linear(h, *T, *(T + 1));
+ break;
+ case 1: // parabolic function TODO: delete
+ teor_var = *T * radius + *(T + 1);
+ break;
+ case 2: // exponential function
+ teor_var = exponential(h, nugget, part_sill, h_range);
+ break;
+ case 3: // spherical function
+ teor_var = spherical(h, nugget, part_sill, h_range);
+ break;
+ case 4: // Gaussian function
+ teor_var = gaussian(radius, nugget, part_sill, h_range);
+ break;
+ case 5:
+ teor_var = *T * h + *(T + 1) * vert + *(T + 2);
+ break;
+ } // end switch (variogram)
+
+ result += teor_var;
}
- else { // univariate variogram or linear bivariate
- variogram = var_par->function;
- if (variogram == 5) {
- h = dr[0];
- vert = dr[1];
- }
- else { // univariate variogram
- radius = SQUARE(dr[0]) + SQUARE(dr[1]) + SQUARE(dr[2]); // weighted dz
- h = sqrt(radius);
+ if (type == 2 && var_par->function != 5) {
+ result -=
+ 0.5 * (var_par->horizontal.sill - nugget +
+ var_par->vertical.sill - nugget);
+ }
- nugget = var_par->nugget;
- part_sill = var_par->sill - nugget;
- h_range = var_par->h_range;
- }
- }
-
- switch (variogram) {
- case 0: // linear function
- teor_var = linear(h, *T, *(T+1));
- break;
- case 1: // parabolic function TODO: delete
- teor_var = *T * radius + *(T+1);
- break;
- case 2: // exponential function
- teor_var = exponential(h, nugget, part_sill, h_range);
- break;
- case 3: // spherical function
- teor_var = spherical(h, nugget, part_sill, h_range);
- break;
- case 4: // Gaussian function
- teor_var = gaussian(radius, nugget, part_sill, h_range);
- break;
- case 5:
- teor_var = *T * h + *(T+1) * vert + *(T+2);
- break;
- } // end switch (variogram)
-
- result += teor_var;
- }
- if (type == 2 && var_par->function != 5) {
- result -= 0.5 * (var_par->horizontal.sill - nugget + var_par->vertical.sill - nugget);
- }
-
- return result;
+ 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)
+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;
+ // 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
+ r0[0] = reg->west + (col + 0.5) * reg->ew_res; // x0
+ r0[1] = reg->north - (row + 0.5) * reg->ns_res; // y0
- if (i3 == TRUE) { // 3D interpolation
- switch (var_par->function) {
- case 5: // bivariate variogram
- r0[2] = reg->bot + (dep + 0.5) * reg->bt_res; // z0
- break;
- default: // univariate (anisotropic) function
- r0[2] = xD->aniso_ratio * (reg->bot + (dep + 0.5) * reg->bt_res); // z0
- break;
- } // end switch
- } // end if
+ if (i3 == TRUE) { // 3D interpolation
+ switch (var_par->function) {
+ case 5: // bivariate variogram
+ r0[2] = reg->bot + (dep + 0.5) * reg->bt_res; // z0
+ break;
+ default: // univariate (anisotropic) function
+ r0[2] = xD->aniso_ratio * (reg->bot + (dep + 0.5) * reg->bt_res); // z0
+ break;
+ } // end switch
+ } // end if
- else { // 2D interpolation
- r0[2] = 0.;
- }
+ else { // 2D interpolation
+ r0[2] = 0.;
+ }
}
// set up elements of G matrix
-mat_struct *set_up_G(struct points *pnts, struct parameters *var_par, struct write *report)
+mat_struct *set_up_G(struct points *pnts, struct parameters *var_par,
+ struct write *report)
{
- // Local variables
- int n = pnts->n; // # of input points
- double *r = pnts->r; // xyz coordinates of input points
- int type = var_par->type; // hz / vert / aniso / bivar
+ // Local variables
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // xyz coordinates of input points
+ int type = var_par->type; // hz / vert / aniso / bivar
- int i, j; // indices of matrix rows/cols
- int n1 = n+1; // # of matrix rows/cols
- double theor_var; // GM element = theor_var(distance)
- double *dr; // dx, dy, dz between point couples
- mat_struct *GM; // GM matrix
+ int i, j; // indices of matrix rows/cols
+ int n1 = n + 1; // # of matrix rows/cols
+ double theor_var; // GM element = theor_var(distance)
+ double *dr; // dx, dy, dz between point couples
+ mat_struct *GM; // GM matrix
- dr = (double *) G_malloc(3 * sizeof(double));
- GM = G_matrix_init(n1, n1, n1); // G[n1][n1] matrix
+ dr = (double *)G_malloc(3 * sizeof(double));
+ GM = G_matrix_init(n1, n1, n1); // G[n1][n1] matrix
- doublereal *md, *mu, *ml, *dbu, *dbl, *m1r, *m1c;
- dbu = &GM->vals[0]; // upper matrix elements
- dbl = &GM->vals[0]; // lower matrix elements
- md = &GM->vals[0]; // diagonal matrix elements
- m1r = &GM->vals[n1-1]; // GM - last matrix row
- m1c = &GM->vals[n1*n]; // GM - last matrix col
-
- // G[n;n]
- for (i = 0; i < n1; i++) { // for elements in each row
- dbu += n1; // new row of the U
- dbl++; // new col of the L
- mu = dbu; // 1st element in the new U row
- ml = dbl; // 1st element in the new L col
- for (j = i; j < n1; j++) { // elements in cols of upper matrix
- if (i != j && i != n && j != n) { // non-diagonal elements except last row/col
- coord_diff(i, j, r, dr); // compute coordinate differences
- if (type == 2) { // bivariate variogram
- dr[0] = sqrt(radius_hz_diff(dr));
- dr[1] = dr[2];
- }
- theor_var = variogram_fction(var_par, dr); // compute GM element
-
- if (isnan(theor_var)) { // not a number:
- if (report->name) {
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
- }
- G_fatal_error(_("Theoretical variogram is NAN..."));
- }
-
- *mu = *ml = (doublereal) theor_var; // set the value to the matrix
- mu += n1; // go to next element in the U row
- ml++; // go to next element in the L col
- } // end non-diagonal elements condition
- } // end j loop
+ doublereal *md, *mu, *ml, *dbu, *dbl, *m1r, *m1c;
- // go to the diagonal element in the next row
- dbu++; // U
- dbl += n1; // L
- *md = 0.0; // set diagonal
- md += n1+1;
+ 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
- // 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
+ // G[n;n]
+ for (i = 0; i < n1; i++) { // for elements in each row
+ dbu += n1; // new row of the U
+ dbl++; // new col of the L
+ mu = dbu; // 1st element in the new U row
+ ml = dbl; // 1st element in the new L col
+ for (j = i; j < n1; j++) { // elements in cols of upper matrix
+ if (i != j && i != n && j != n) { // non-diagonal elements except last row/col
+ coord_diff(i, j, r, dr); // compute coordinate differences
+ if (type == 2) { // bivariate variogram
+ dr[0] = sqrt(radius_hz_diff(dr));
+ dr[1] = dr[2];
+ }
+ theor_var = variogram_fction(var_par, dr); // compute GM element
- free(dr);
- return GM;
+ if (isnan(theor_var)) { // not a number:
+ if (report->name) {
+ fprintf(report->fp,
+ "Error (see standard output). Process killed...");
+ fclose(report->fp);
+ }
+ G_fatal_error(_("Theoretical variogram is NAN..."));
+ }
+
+ *mu = *ml = (doublereal) theor_var; // set the value to the matrix
+ mu += n1; // go to next element in the U row
+ ml++; // go to next element in the L col
+ } // end non-diagonal elements condition
+ } // end j loop
+
+ // go to the diagonal element in the next row
+ dbu++; // U
+ dbl += n1; // L
+ *md = 0.0; // set diagonal
+ md += n1 + 1;
+
+ // Set ones
+ if (i < n1 - 1) { // all elements except the last one...
+ *m1r = *m1c = 1.0; // ... shall be 1
+ m1r += n1; // go to next col in last row
+ m1c++; // go to next row in last col
+ } // end "last 1" condition
+ } // end i loop
+
+ free(dr);
+ return GM;
}
// make G submatrix for rellevant points
-mat_struct *submatrix(struct ilist *index, mat_struct *GM_all, struct write *report)
+mat_struct *submatrix(struct ilist * index, mat_struct * GM_all,
+ struct write * report)
{
- // Local variables
- int n = index->n_values; // # of selected points
- mat_struct *GM = GM_all; // whole G matrix
+ // Local variables
+ int n = index->n_values; // # of selected points
+ mat_struct *GM = GM_all; // whole G matrix
- int i, j, N1 = GM->rows, n1 = n+1, *dinR, *dini, *dinj;
- doublereal *dbo, *dbx, *dbu, *dbl, *md, *mu, *ml, *m1r, *m1c;
+ int i, j, N1 = GM->rows, n1 = n + 1, *dinR, *dini, *dinj;
+ doublereal *dbo, *dbx, *dbu, *dbl, *md, *mu, *ml, *m1r, *m1c;
- mat_struct *sub; // new submatrix
- sub = G_matrix_init(n1, n1, n1);
+ mat_struct *sub; // new submatrix
- if (sub == NULL) {
- if (report->name) { // close report file
- fprintf(report->fp, "Error (see standard output). Process killed...");
- fclose(report->fp);
+ 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..."));
}
- G_fatal_error(_("Unable to initialize G-submatrix..."));
- }
-
- // Initialize indexes: GM[0;1]
- // - cols
- dini = &index->value[0]; // i - set processing column
-
- // initialize new col
- dinR = &index->value[1]; // return to first processing row - lower GM
- dbo = &GM->vals[*dini * N1 + *dinR]; // 1st value in GM_all
- md = &sub->vals[0]; // GM - diagonal
- dbu = &sub->vals[n1]; // GM - upper
- dbl = &sub->vals[1]; // GM - lower
- m1r = &sub->vals[n1-1]; // GM - last row
- m1c = &sub->vals[n1*n]; // GM - last col
+ // Initialize indexes: GM[0;1]
+ // - cols
+ dini = &index->value[0]; // i - set processing column
- for (i=0; i<=n; i++) { // set up cols of submatrix
- // original matrix
- dbx = dbo; // new col in GM_all
- dinj = dinR; // orig: inicialize element (row) in (col) - upper matrix
- // submatrix
- mu = dbu; // sub: start new column
- ml = dbl; // sub: start new row
- for (j=i+1; j<n; j++) { // for rows
- //submatrix
- *mu = *ml = *dbx; // general elements
- mu += n1; // go to element in next column
- ml++; // go to element in next row
- // Original matrix
- dinj++; // set next ind
- dbx += *dinj - *(dinj-1);
- } // end j
+ // initialize new col
+ dinR = &index->value[1]; // return to first processing row - lower GM
- // Original matrix
- dini++; // set next ind
-
- dbu += n1+1; // new col in GM
- dbl += n1+1; // new row in GM
-
- dinR++;
+ 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
- dbo += (*dini - *(dini-1)) * N1 + (*dinR - *(dinR-1));
- // Set ones
- *m1r = *m1c = 1.0;
- m1r += n1;
- m1c++;
- // Set diagonals
- *md = 0.0;
- md += n1+1;
- } // end i
+ for (i = 0; i <= n; i++) { // set up cols of submatrix
+ // original matrix
+ dbx = dbo; // new col in GM_all
+ dinj = dinR; // orig: inicialize element (row) in (col) - upper matrix
+ // submatrix
+ mu = dbu; // sub: start new column
+ ml = dbl; // sub: start new row
+ for (j = i + 1; j < n; j++) { // for rows
+ //submatrix
+ *mu = *ml = *dbx; // general elements
+ mu += n1; // go to element in next column
+ ml++; // go to element in next row
+ // Original matrix
+ dinj++; // set next ind
+ dbx += *dinj - *(dinj - 1);
+ } // end j
- return sub;
+ // Original matrix
+ dini++; // set next ind
+
+ dbu += n1 + 1; // new col in GM
+ dbl += n1 + 1; // new row in GM
+
+ dinR++;
+
+ dbo += (*dini - *(dini - 1)) * N1 + (*dinR - *(dinR - 1));
+ // Set ones
+ *m1r = *m1c = 1.0;
+ m1r += n1;
+ m1c++;
+ // Set diagonals
+ *md = 0.0;
+ md += n1 + 1;
+ } // end i
+
+ return sub;
}
// set up g0 vector
-mat_struct *set_up_g0(struct int_par *xD, struct points *pnts, struct ilist *index, double *r0, struct parameters *var_par)
+mat_struct *set_up_g0(struct int_par * xD, struct points * pnts,
+ struct ilist * index, double *r0,
+ struct parameters * var_par)
{
- // Local variables
- int i3 = xD->i3; // interpolation: 2D or 3D
- int type = var_par->type; // variogram: hz / vert / aniso / bivar
- double *r; // xyz coordinates of input points
+ // Local variables
+ int i3 = xD->i3; // interpolation: 2D or 3D
+ int type = var_par->type; // variogram: hz / vert / aniso / bivar
+ double *r; // xyz coordinates of input points
- int n; // # of points (all or selected)
- int n_ind = index->n_values; // # of selected points
- int *lind; // pointer to indices of selected points
+ int n; // # of points (all or selected)
+ int n_ind = index->n_values; // # of selected points
+ int *lind; // pointer to indices of selected points
- if (n_ind == 0) { // no selected points:
- n = pnts->n; // use all points
- r = &pnts->r[0];
- }
- else { // any selected points:
- n = n_ind; // reduce # of input points
- lind = &index->value[0];
- r = &pnts->r[*lind * 3];
- }
+ if (n_ind == 0) { // no selected points:
+ n = pnts->n; // use all points
+ r = &pnts->r[0];
+ }
+ else { // any selected points:
+ n = n_ind; // reduce # of input points
+ lind = &index->value[0];
+ r = &pnts->r[*lind * 3];
+ }
- int n1 = n + 1;
-
- int i; // index of elements and length of g0 vector
- double *dr; // coordinate differences dx, dy, dz between couples of points
- mat_struct *g0; // vector of estimated differences between known and unknown values
+ int n1 = n + 1;
- dr = (double *) G_malloc(3 * sizeof(double)); // Coordinate differences
- g0 = G_matrix_init(n1,1,n1);
+ int i; // index of elements and length of g0 vector
+ double *dr; // coordinate differences dx, dy, dz between couples of points
+ mat_struct *g0; // vector of estimated differences between known and unknown values
- doublereal *g = g0->vals;
-
- for (i=0; i<n; i++) { // count of input points
- // Coord diffs (input points and cell/voxel center)
- dr[0] = *r - *r0; // dx
- dr[1] = *(r+1) - *(r0+1); // dy
- switch (i3) {
- case FALSE:
- // Cell value diffs estimation using linear variogram
- dr[2] = 0.0; // !!! optimalize - not to be necessary to use this line
- break;
- case TRUE:
- dr[2] = *(r+2) - *(r0+2); // dz
- break;
- }
-
- if (type == 2) { // bivariate variogram
- dr[0] = sqrt(radius_hz_diff(dr));
- dr[1] = dr[2];
- }
+ dr = (double *)G_malloc(3 * sizeof(double)); // Coordinate differences
+ g0 = G_matrix_init(n1, 1, n1);
- *g = variogram_fction(var_par, dr); // compute GM element using theoretical variogram function
- g++;
+ doublereal *g = g0->vals;
- if (n_ind == 0) { // no selected points:
- r += 3; // use each point
- }
- else { // selected points:
- lind++; // go to the next index
- r += 3 * (*lind - *(lind-1)); // use next selected point
- }
- } // end i loop
+ for (i = 0; i < n; i++) { // count of input points
+ // Coord diffs (input points and cell/voxel center)
+ dr[0] = *r - *r0; // dx
+ dr[1] = *(r + 1) - *(r0 + 1); // dy
+ switch (i3) {
+ case FALSE:
+ // Cell value diffs estimation using linear variogram
+ dr[2] = 0.0; // !!! optimalize - not to be necessary to use this line
+ break;
+ case TRUE:
+ dr[2] = *(r + 2) - *(r0 + 2); // dz
+ break;
+ }
- // condition: sum of weights = 1
- *g = 1.0; // g0 [n1x1]
+ if (type == 2) { // bivariate variogram
+ dr[0] = sqrt(radius_hz_diff(dr));
+ dr[1] = dr[2];
+ }
- return g0;
+ *g = variogram_fction(var_par, dr); // compute GM element using theoretical variogram function
+ g++;
+
+ if (n_ind == 0) { // no selected points:
+ r += 3; // use each point
+ }
+ else { // selected points:
+ lind++; // go to the next index
+ r += 3 * (*lind - *(lind - 1)); // use next selected point
+ }
+ } // end i loop
+
+ // condition: sum of weights = 1
+ *g = 1.0; // g0 [n1x1]
+
+ return g0;
}
// compute result to be located on reference point
-double result(struct points *pnts, struct ilist *index, mat_struct *w0)
+double result(struct points *pnts, struct ilist *index, mat_struct * w0)
{
- // local variables
- int n; // # of points
- double *vo; // pointer to input values
+ // local variables
+ int n; // # of points
+ double *vo; // pointer to input values
- int n_ind = index->n_values; // # of selected points
- int *indo; // pointer to indices of selected pt
+ int n_ind = index->n_values; // # of selected points
+ int *indo; // pointer to indices of selected pt
- // Setup n depending on NN points selection
- n = n_ind == 0 ? pnts->n : n_ind;
+ // Setup n depending on NN points selection
+ n = n_ind == 0 ? pnts->n : n_ind;
- int i;
- mat_struct *ins, *w, *rslt_OK;
- doublereal *vt, *wo, *wt;
+ int i;
+ mat_struct *ins, *w, *rslt_OK;
+ doublereal *vt, *wo, *wt;
- ins = G_matrix_init(n, 1, n); // matrix of selected vals
- w = G_matrix_init(1, n, 1); // matrix of selected weights
+ ins = G_matrix_init(n, 1, n); // matrix of selected vals
+ w = G_matrix_init(1, n, 1); // matrix of selected weights
- if (n_ind == 0) { // there are no selected points:
- vo = &pnts->invals[0]; // use all input values
- }
- else { // there are selected points:
- indo = &index->value[0]; // original indexes
- vo = &pnts->invals[*indo]; // original input values
- }
-
- wo = &w0->vals[0]; // pointer to weights of input vals
- vt = &ins->vals[0]; // pointer to selected inputs values
- wt = &w->vals[0]; // pointer to weights without the last one
-
- for (i=0; i<n; i++) { // for each input point:
- *vt = *vo; // subval = selected original
- *wt = *wo; // sub weight = selected original
- //G_debug(0,"%f %f", *vt, *wt);
-
- // go to the next:
- if (n_ind == 0) { // use all points:
- vo++; // input value
+ if (n_ind == 0) { // there are no selected points:
+ vo = &pnts->invals[0]; // use all input values
}
- else { // use selected points:
- indo++; // index of selected point
- vo += *indo - *(indo-1); // selected input value
+ else { // there are selected points:
+ indo = &index->value[0]; // original indexes
+ vo = &pnts->invals[*indo]; // original input values
}
- vt++; // element of value matrix
- wo++; // weight of the value
- wt++; // element of weight matrix
- } // end i for loop
- rslt_OK = G_matrix_product(w,ins); // interpolated value
-
- G_matrix_free(w);
- G_matrix_free(ins);
+ wo = &w0->vals[0]; // pointer to weights of input vals
+ vt = &ins->vals[0]; // pointer to selected inputs values
+ wt = &w->vals[0]; // pointer to weights without the last one
- return rslt_OK->vals[0];
+ for (i = 0; i < n; i++) { // for each input point:
+ *vt = *vo; // subval = selected original
+ *wt = *wo; // sub weight = selected original
+ //G_debug(0,"%f %f", *vt, *wt);
+
+ // go to the next:
+ if (n_ind == 0) { // use all points:
+ vo++; // input value
+ }
+ else { // use selected points:
+ indo++; // index of selected point
+ vo += *indo - *(indo - 1); // selected input value
+ }
+ vt++; // element of value matrix
+ wo++; // weight of the value
+ wt++; // element of weight matrix
+ } // end i for loop
+
+ rslt_OK = G_matrix_product(w, ins); // interpolated value
+
+ G_matrix_free(w);
+ G_matrix_free(ins);
+
+ return rslt_OK->vals[0];
}
// validation
-void crossvalidation(struct int_par *xD, struct points *pnts, struct parameters *var_par)
+void crossvalidation(struct int_par *xD, struct points *pnts,
+ struct parameters *var_par)
{
- int n = pnts->n; // # of input points
- double *r = pnts->r; // coordinates of points
- int i3 = xD->i3;
- double *vals = pnts->invals; // input values
+ int n = pnts->n; // # of input points
+ double *r = pnts->r; // coordinates of points
+ int i3 = xD->i3;
+ double *vals = pnts->invals; // input values
- int type = var_par->type;
- double ratio = type == 3 ? xD->aniso_ratio : 1.; // anisotropic ratio
- mat_struct *GM = var_par->GM; // GM = theor_var(distances)
-
- int i;
- int n_vals;
- double max_dist = type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
- double max_dist_vert = type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
-
- double *search; // coordinates of the search point
- struct ilist *list;
- mat_struct *GM_sub;
- mat_struct *GM_Inv, *g0, *w0;
- double rslt_OK;
+ int type = var_par->type;
+ double ratio = type == 3 ? xD->aniso_ratio : 1.; // anisotropic ratio
+ mat_struct *GM = var_par->GM; // GM = theor_var(distances)
- struct write *crossvalid = &xD->crossvalid;
- struct write *report = &xD->report;
- double *normal, *absval, *norm, *av;
+ int i;
+ int n_vals;
+ double max_dist =
+ type == 2 ? var_par->horizontal.max_dist : var_par->max_dist;
+ double max_dist_vert =
+ type == 2 ? var_par->vertical.max_dist : var_par->max_dist;
- search = (double *) G_malloc(3 * sizeof(double));
- normal = (double *) G_malloc(n * sizeof(double));
- absval = (double *) G_malloc(n * sizeof(double));
- norm = &normal[0];
- av = &absval[0];
+ double *search; // coordinates of the search point
+ struct ilist *list;
+ mat_struct *GM_sub;
+ mat_struct *GM_Inv, *g0, *w0;
+ double rslt_OK;
- FILE *fp;
- if (crossvalid->name) { // CV report file available:
- fp = fopen(crossvalid->name, "w");
- }
+ struct write *crossvalid = &xD->crossvalid;
+ struct write *report = &xD->report;
+ double *normal, *absval, *norm, *av;
- for (i=0; i<n; i++) { // for each input point [r0]:
- list = G_new_ilist(); // create list of overlapping rectangles
+ search = (double *)G_malloc(3 * sizeof(double));
+ normal = (double *)G_malloc(n * sizeof(double));
+ absval = (double *)G_malloc(n * sizeof(double));
+ norm = &normal[0];
+ av = &absval[0];
- search = &pnts->r[3 * i];
- if (i3 == TRUE) {
- list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
+ FILE *fp;
+
+ if (crossvalid->name) { // CV report file available:
+ fp = fopen(crossvalid->name, "w");
}
- else {
- list = find_NNs_within(2, search, pnts, max_dist, max_dist_vert);
- }
-
- n_vals = list->n_values; // # of overlapping rectangles
- if (n_vals > 0 ) { // if positive:
- correct_indices(list, r, pnts, var_par);
+ for (i = 0; i < n; i++) { // for each input point [r0]:
+ list = G_new_ilist(); // create list of overlapping rectangles
- GM_sub = submatrix(list, GM, &xD->report); // create submatrix using indices
- GM_Inv = G_matrix_inverse(GM_sub); // inverse matrix
- G_matrix_free(GM_sub);
-
- g0 = set_up_g0(xD, pnts, list, r, var_par); // Diffs inputs - unknowns (incl. cond. 1))
- w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
-
- G_matrix_free(g0);
- G_matrix_free(GM_Inv);
+ search = &pnts->r[3 * i];
+ if (i3 == TRUE) {
+ list = find_NNs_within(3, search, pnts, max_dist, max_dist_vert);
+ }
+ else {
+ list = find_NNs_within(2, search, pnts, max_dist, max_dist_vert);
+ }
- rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
- G_matrix_free(w0);
+ n_vals = list->n_values; // # of overlapping rectangles
- //Create output
- *norm = rslt_OK - *vals; // differences between input and interpolated values
- *av = fabsf(*norm); // absolute values of the differences (quantile computation)
+ if (n_vals > 0) { // if positive:
+ correct_indices(list, r, pnts, var_par);
- if (xD->i3 == TRUE) { // 3D interpolation:
- fprintf(fp,"%d %.3f %.3f %.2f %f %f %f\n", i, *r, *(r+1), *(r+2) / ratio, pnts->invals[i], rslt_OK, *norm);
- }
- else { // 2D interpolation:
- fprintf(fp,"%d %.3f %.3f %f %f %f\n", i, *r, *(r+1), *vals, rslt_OK, *norm);
- }
- } // end if: n_vals > 0
+ GM_sub = submatrix(list, GM, &xD->report); // create submatrix using indices
+ GM_Inv = G_matrix_inverse(GM_sub); // inverse matrix
+ G_matrix_free(GM_sub);
- else { // n_vals == 0:
- fprintf(fp,"%d. point does not have neighbours in given radius\n", i);
- }
+ g0 = set_up_g0(xD, pnts, list, r, var_par); // Diffs inputs - unknowns (incl. cond. 1))
+ w0 = G_matrix_product(GM_Inv, g0); // Vector of weights, condition SUM(w) = 1 in last row
- r += 3; // next rectangle
- vals++;
- norm++;
- av++;
+ G_matrix_free(g0);
+ G_matrix_free(GM_Inv);
- G_free_ilist(list); // free list memory
- } // end i for loop
+ rslt_OK = result(pnts, list, w0); // Estimated cell/voxel value rslt_OK = w x inputs
+ G_matrix_free(w0);
- fclose(fp);
- G_message(_("Cross validation results have been written into <%s>"), crossvalid->name);
+ //Create output
+ *norm = rslt_OK - *vals; // differences between input and interpolated values
+ *av = fabsf(*norm); // absolute values of the differences (quantile computation)
- if (report->name) {
- double quant95;
- fprintf(report->fp, "\n************************************************\n");
- fprintf(report->fp, "*** Cross validation results ***\n");
+ if (xD->i3 == TRUE) { // 3D interpolation:
+ fprintf(fp, "%d %.3f %.3f %.2f %f %f %f\n", i, *r, *(r + 1),
+ *(r + 2) / ratio, pnts->invals[i], rslt_OK, *norm);
+ }
+ else { // 2D interpolation:
+ fprintf(fp, "%d %.3f %.3f %f %f %f\n", i, *r, *(r + 1), *vals,
+ rslt_OK, *norm);
+ }
+ } // end if: n_vals > 0
- test_normality(n, normal, report);
-
- fprintf(report->fp, "Quantile of absolute values\n");
- quant95 = quantile(0.95, n, absval, report);
- }
+ else { // n_vals == 0:
+ fprintf(fp,
+ "%d. point does not have neighbours in given radius\n",
+ i);
+ }
+
+ r += 3; // next rectangle
+ vals++;
+ norm++;
+ av++;
+
+ G_free_ilist(list); // free list memory
+ } // end i for loop
+
+ fclose(fp);
+ G_message(_("Cross validation results have been written into <%s>"),
+ crossvalid->name);
+
+ if (report->name) {
+ double 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);
+ }
}
Modified: grass-addons/grass7/vector/v.kriging/utils_raster.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_raster.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/utils_raster.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -3,58 +3,65 @@
/* 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);
+ 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);
+ /* 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)
+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;
+ // 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;
+ int pass = 0; /* Number of processed cells */
- 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);
+ 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;
}
- pass++;
- break;
- }
- return pass;
+ return pass;
}
-
Modified: grass-addons/grass7/vector/v.kriging/utils_write.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/utils_write.c 2015-07-09 08:29:00 UTC (rev 65550)
+++ grass-addons/grass7/vector/v.kriging/utils_write.c 2015-07-09 09:17:29 UTC (rev 65551)
@@ -2,319 +2,335 @@
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;
- }
+ 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;
+ 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;
- }
+ 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)
+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;
- }
+ 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);
+ 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);
+ char type[12];
- fprintf(report->fp, "\n************************************************\n");
- fprintf(report->fp, "*** Variogram settings - %s ***\n\n", type);
+ variogram_type(code, type);
+
+ fprintf(report->fp,
+ "\n************************************************\n");
+ fprintf(report->fp, "*** Variogram settings - %s ***\n\n", type);
}
void write2file_varSets(struct write *report, struct parameters *var_pars)
{
- char type[12];
- variogram_type(var_pars->type, type);
+ char type[12];
- fprintf(report->fp, "- number of lags in %s direction: %d\n", type, var_pars->nLag);
- fprintf(report->fp, "- lag distance (%s): %f\n", type, var_pars->lag);
+ variogram_type(var_pars->type, type);
- if (var_pars->type == 2) {
- fprintf(report->fp, "- number of lags in vertical direction: %d\n", var_pars->nLag_vert);
- fprintf(report->fp, "- lag distance (vertical): %f\n", var_pars->lag);
- }
+ fprintf(report->fp, "- number of lags in %s direction: %d\n", type,
+ var_pars->nLag);
+ fprintf(report->fp, "- lag distance (%s): %f\n", type, var_pars->lag);
- fprintf(report->fp, "- azimuth: %f°\n", RAD2DEG(var_pars->dir));
- fprintf(report->fp, "- angular tolerance: %f°\n", RAD2DEG(var_pars->td));
-
- switch (var_pars->type) {
- case 2:
- fprintf(report->fp, "- maximum distance in horizontal direction (1/3): %f\n", var_pars->max_dist);
- break;
- default:
- fprintf(report->fp, "- maximum distance (1/3): %f\n", var_pars->max_dist);
- break;
- }
+ if (var_pars->type == 2) {
+ fprintf(report->fp, "- number of lags in vertical direction: %d\n",
+ var_pars->nLag_vert);
+ fprintf(report->fp, "- lag distance (vertical): %f\n", var_pars->lag);
+ }
+
+ fprintf(report->fp, "- azimuth: %f°\n", RAD2DEG(var_pars->dir));
+ fprintf(report->fp, "- angular tolerance: %f°\n", RAD2DEG(var_pars->td));
+
+ switch (var_pars->type) {
+ case 2:
+ fprintf(report->fp,
+ "- maximum distance in horizontal direction (1/3): %f\n",
+ var_pars->max_dist);
+ break;
+ default:
+ fprintf(report->fp, "- maximum distance (1/3): %f\n",
+ var_pars->max_dist);
+ break;
+ }
}
-void write2file_variogram_E(struct int_par *xD, struct parameters *var_pars, mat_struct *c_M)
+void write2file_variogram_E(struct int_par *xD, struct parameters *var_pars,
+ mat_struct * c_M)
{
- int i, j;
- int n_h;
- double *h, *vert;
- double *c;
- double *gamma;
- struct write *report = &xD->report;
+ int i, j;
+ int n_h;
+ double *h, *vert;
+ double *c;
+ double *gamma;
+ struct write *report = &xD->report;
- char type[12];
- variogram_type(var_pars->type, type);
+ char type[12];
- c = &c_M->vals[0];
- gamma = &var_pars->gamma->vals[0];
-
- fprintf(report->fp, "\n************************************************\n");
- fprintf(report->fp, "*** Experimental variogram - %s ***\n\n", type);
-
- if (var_pars->type == 2) { // bivariate variogram:
- vert = &var_pars->vertical.h[0];
-
- for (i=0; i<var_pars->vertical.nLag; i++) { // write header - verts
- if (i == 0) {
- fprintf(report->fp, " lagV ||"); // column for h
- }
- fprintf(report->fp, " %f ||", *vert);
- vert++;
- }
- fprintf(report->fp, "\n");
+ variogram_type(var_pars->type, type);
- for (i=0; i < var_pars->vertical.nLag; i++) { // write header - h c gamma
- if (i == 0) {
- fprintf(report->fp, " lagHZ ||"); // column for h
- }
- fprintf(report->fp, " c | ave ||");
+ c = &c_M->vals[0];
+ gamma = &var_pars->gamma->vals[0];
+
+ fprintf(report->fp,
+ "\n************************************************\n");
+ fprintf(report->fp, "*** Experimental variogram - %s ***\n\n", type);
+
+ if (var_pars->type == 2) { // bivariate variogram:
+ vert = &var_pars->vertical.h[0];
+
+ for (i = 0; i < var_pars->vertical.nLag; i++) { // write header - verts
+ if (i == 0) {
+ fprintf(report->fp, " lagV ||"); // column for h
+ }
+ fprintf(report->fp, " %f ||", *vert);
+ vert++;
+ }
+ fprintf(report->fp, "\n");
+
+ for (i = 0; i < var_pars->vertical.nLag; i++) { // write header - h c gamma
+ if (i == 0) {
+ fprintf(report->fp, " lagHZ ||"); // column for h
+ }
+ fprintf(report->fp, " c | ave ||");
+ }
+ fprintf(report->fp, "\n");
+
+ for (i = 0; i < var_pars->vertical.nLag; i++) { // write header - h c gamma
+ if (i == 0) {
+ fprintf(report->fp, "--------||"); // column for h
+ }
+ fprintf(report->fp, "-----------");
+ }
+ fprintf(report->fp, "\n");
}
- fprintf(report->fp, "\n");
+ else { // univariate variogram
+ fprintf(report->fp, " lag || # of pairs | average\n");
+ fprintf(report->fp, "------------------------------------\n");
+ }
- for (i=0; i<var_pars->vertical.nLag; i++) { // write header - h c gamma
- if (i == 0) {
- fprintf(report->fp, "--------||"); // column for h
- }
- fprintf(report->fp, "-----------");
+ // Write values
+ h = var_pars->type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
+ n_h = var_pars->type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
+ for (i = 0; i < n_h; i++) {
+ fprintf(report->fp, "%f ||", *h);
+ if (var_pars->type == 2) { // bivariate variogram
+ for (j = 0; j < var_pars->vertical.nLag; j++) {
+ fprintf(report->fp, " %d | %f ||", (int)*c, *gamma);
+ c++;
+ gamma++;
+ } // end for j
+ fprintf(report->fp, "\n");
+ for (j = 0; j < var_pars->vertical.nLag; j++) {
+ fprintf(report->fp, "-----------");
+ }
+ fprintf(report->fp, "\n");
+ }
+ else { // univariate variogram
+ fprintf(report->fp, " %d | %f\n", (int)*c, *gamma);
+ c++;
+ gamma++;
+ }
+ h++;
}
- fprintf(report->fp, "\n");
- }
- else { // univariate variogram
- fprintf(report->fp, " lag || # of pairs | average\n");
fprintf(report->fp, "------------------------------------\n");
- }
-
- // Write values
- h = var_pars->type == 2 ? &var_pars->horizontal.h[0] : &var_pars->h[0];
- n_h = var_pars->type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
- for (i=0; i < n_h; i++) {
- fprintf(report->fp, "%f ||", *h);
- if (var_pars->type == 2) { // bivariate variogram
- for (j=0; j<var_pars->vertical.nLag; j++) {
- fprintf(report->fp, " %d | %f ||", (int) *c, *gamma);
- c++;
- gamma++;
- } // end for j
- fprintf(report->fp, "\n");
- for (j=0; j<var_pars->vertical.nLag; j++) {
- fprintf(report->fp, "-----------");
- }
- fprintf(report->fp, "\n");
- }
- else { // univariate variogram
- fprintf(report->fp, " %d | %f\n", (int) *c, *gamma);
- c++;
- gamma++;
- }
- h++;
- }
- fprintf(report->fp, "------------------------------------\n");
}
void write2file_variogram_T(struct write *report)
{
- fprintf(report->fp, "\n");
- fprintf(report->fp, "************************************************\n");
- fprintf(report->fp, "*** Theoretical variogram ***\n\n");
+ fprintf(report->fp, "\n");
+ fprintf(report->fp, "************************************************\n");
+ fprintf(report->fp, "*** Theoretical variogram ***\n\n");
}
void write_temporary2file(struct int_par *xD, struct parameters *var_pars)
{
- // local variables
- int type = var_pars->type;
- int nLag = type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
- int nLag_vert = type == 2 ? var_pars->vertical.nLag : var_pars->nLag;
- double *h = type == 2 ? var_pars->horizontal.h: var_pars->h;
- double *vert = type == 2 ? var_pars->vertical.h: var_pars->h;;
- double *gamma = var_pars->gamma->vals;
- double sill = var_pars->sill;
+ // local variables
+ int type = var_pars->type;
+ int nLag = type == 2 ? var_pars->horizontal.nLag : var_pars->nLag;
+ int nLag_vert = type == 2 ? var_pars->vertical.nLag : var_pars->nLag;
+ double *h = type == 2 ? var_pars->horizontal.h : var_pars->h;
+ double *vert = type == 2 ? var_pars->vertical.h : var_pars->h;;
+ double *gamma = var_pars->gamma->vals;
+ double sill = var_pars->sill;
- int i; // index
- FILE *fp;
- int file_length;
-
- /* Introduction */
- switch (type) {
- case 0: // horizontal variogram
- fp = fopen("variogram_hz_tmp.txt", "w");
- if (xD->report.name) { // write name of report file
- file_length = strlen(xD->report.name);
- if (file_length < 4) { // 4 types of variogram
- G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
- }
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
- }
- fprintf(fp,"%d\n", var_pars->type); // write # of lags
- break;
+ int i; // index
+ FILE *fp;
+ int file_length;
- case 1: // vertical variogram
- fp = fopen("variogram_vert_tmp.txt", "w");
- if (xD->report.write2file == TRUE) { // write name of report file
- file_length = strlen(xD->report.name);
- if (file_length < 3) {
- G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
- }
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
- }
- fprintf(fp,"%d\n", var_pars->type); // write type
- break;
+ /* Introduction */
+ switch (type) {
+ case 0: // horizontal variogram
+ fp = fopen("variogram_hz_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) { // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ }
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp, "%d\n", var_pars->type); // write # of lags
+ break;
- case 2: // bivariate variogram
- fp = fopen("variogram_final_tmp.txt", "w");
+ case 1: // vertical variogram
+ fp = fopen("variogram_vert_tmp.txt", "w");
+ if (xD->report.write2file == TRUE) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 3) {
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ }
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp, "%d\n", var_pars->type); // write type
+ break;
- if (xD->report.write2file == TRUE) { // write name of report file
- file_length = strlen(xD->report.name);
- if (file_length < 4) { // 4 types of variogram
- G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
- }
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
- }
+ case 2: // bivariate variogram
+ fp = fopen("variogram_final_tmp.txt", "w");
- fprintf(fp,"%d\n", var_pars->type); // write type
- fprintf(fp,"%d\n", var_pars->vertical.nLag); // write # of lags
- fprintf(fp,"%f\n", var_pars->vertical.lag); // write size of lag
- fprintf(fp,"%f\n", var_pars->vertical.max_dist); // write maximum distance
+ if (xD->report.write2file == TRUE) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) { // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 2 characters...")); // todo: error
+ }
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
- fprintf(fp,"%d\n", var_pars->horizontal.nLag); // write # of lags
- fprintf(fp,"%f\n", var_pars->horizontal.lag); // write size of lag
- fprintf(fp,"%f\n", var_pars->horizontal.max_dist); // write maximum distance
- break;
+ fprintf(fp, "%d\n", var_pars->type); // write type
+ fprintf(fp, "%d\n", var_pars->vertical.nLag); // write # of lags
+ fprintf(fp, "%f\n", var_pars->vertical.lag); // write size of lag
+ fprintf(fp, "%f\n", var_pars->vertical.max_dist); // write maximum distance
- case 3: // anisotropic variogram
- fp = fopen("variogram_final_tmp.txt", "w");
- if (xD->report.name) { // write name of report file
- file_length = strlen(xD->report.name);
- if (file_length < 4) { // 4 types of variogram
- G_fatal_error(_("File name must contain more than 4 characters...")); // todo: error
- }
- fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ fprintf(fp, "%d\n", var_pars->horizontal.nLag); // write # of lags
+ fprintf(fp, "%f\n", var_pars->horizontal.lag); // write size of lag
+ fprintf(fp, "%f\n", var_pars->horizontal.max_dist); // write maximum distance
+ break;
+
+ case 3: // anisotropic variogram
+ fp = fopen("variogram_final_tmp.txt", "w");
+ if (xD->report.name) { // write name of report file
+ file_length = strlen(xD->report.name);
+ if (file_length < 4) { // 4 types of variogram
+ G_fatal_error(_("File name must contain more than 4 characters...")); // todo: error
+ }
+ fprintf(fp, "%d 9 %s\n", file_length, xD->report.name);
+ }
+ fprintf(fp, "%d\n", var_pars->type); // write type
+ fprintf(fp, "%f\n", xD->aniso_ratio); // write ratio of anisotropy
+ break;
}
- fprintf(fp,"%d\n", var_pars->type); // write type
- fprintf(fp,"%f\n", xD->aniso_ratio); // write ratio of anisotropy
- break;
- }
-
- if (type != 2) {
- fprintf(fp,"%d\n", nLag); // write # of lags
- fprintf(fp,"%f\n", var_pars->lag); // write size of lag
- fprintf(fp,"%f\n", var_pars->max_dist); // write maximum distance
- }
- if (type != 1) {
- fprintf(fp,"%f\n", var_pars->td); // write maximum distance
- }
- switch (type) {
- case 2: // bivariate variogram
- nLag = var_pars->horizontal.nLag;
- nLag_vert = var_pars->vertical.nLag;
- // write h
- for (i=0; i<nLag; i++) { // write experimental variogram
- fprintf(fp,"%f\n", *h);
- h++;
+ if (type != 2) {
+ fprintf(fp, "%d\n", nLag); // write # of lags
+ fprintf(fp, "%f\n", var_pars->lag); // write size of lag
+ fprintf(fp, "%f\n", var_pars->max_dist); // write maximum distance
}
- // write vert
- for (i=0; i<nLag_vert; i++) { // write experimental variogram
- fprintf(fp,"%f\n", *vert);
- vert++;
+ if (type != 1) {
+ fprintf(fp, "%f\n", var_pars->td); // write maximum distance
}
- // write gamma
- for (i=0; i<nLag * nLag_vert; i++) { // write experimental variogram
- fprintf(fp,"%f\n", *gamma);
- gamma++;
+
+ switch (type) {
+ case 2: // bivariate variogram
+ nLag = var_pars->horizontal.nLag;
+ nLag_vert = var_pars->vertical.nLag;
+ // write h
+ for (i = 0; i < nLag; i++) { // write experimental variogram
+ fprintf(fp, "%f\n", *h);
+ h++;
+ }
+ // write vert
+ for (i = 0; i < nLag_vert; i++) { // write experimental variogram
+ fprintf(fp, "%f\n", *vert);
+ vert++;
+ }
+ // write gamma
+ for (i = 0; i < nLag * nLag_vert; i++) { // write experimental variogram
+ fprintf(fp, "%f\n", *gamma);
+ gamma++;
+ }
+ fprintf(fp, "%f\n", var_pars->horizontal.sill); // write sill
+ fprintf(fp, "%f\n", var_pars->vertical.sill); // write sill
+ break;
+ default:
+ for (i = 0; i < nLag; i++) { // write experimental variogram
+ fprintf(fp, "%f %f\n", *h, *gamma);
+ h++;
+ gamma++;
+ }
+ fprintf(fp, "%f", sill); // write sill
+ break;
}
- fprintf(fp,"%f\n", var_pars->horizontal.sill);// write sill
- fprintf(fp,"%f\n", var_pars->vertical.sill);// write sill
- break;
- default:
- for (i=0; i<nLag; i++) { // write experimental variogram
- fprintf(fp,"%f %f\n", *h, *gamma);
- h++;
- gamma++;
- }
- fprintf(fp,"%f", sill);// write sill
- break;
- }
- fclose(fp);
+ fclose(fp);
}
More information about the grass-commit
mailing list