[GRASS-SVN] r65501 - grass-addons/grass7/vector/v.kriging
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 20 04:23:47 PDT 2015
Author: evas
Date: 2015-06-20 04:23:47 -0700 (Sat, 20 Jun 2015)
New Revision: 65501
Added:
grass-addons/grass7/vector/v.kriging/main.c
Modified:
grass-addons/grass7/vector/v.kriging/Makefile
grass-addons/grass7/vector/v.kriging/local_proto.h
grass-addons/grass7/vector/v.kriging/v.kriging.html
Log:
v.kriging: kd-tree (PCL library) replaced by R-tree (GRASS GIS) and 2D kriging enabled (currently: inaccurate results for sparse and spatially heterogeneous data)
Modified: grass-addons/grass7/vector/v.kriging/Makefile
===================================================================
--- grass-addons/grass7/vector/v.kriging/Makefile 2015-06-20 09:32:39 UTC (rev 65500)
+++ grass-addons/grass7/vector/v.kriging/Makefile 2015-06-20 11:23:47 UTC (rev 65501)
@@ -1,12 +1,12 @@
-MODULE_TOPDIR = ../..
+MODULE_TOPDIR = /home/evergreen/grass7/dist.x86_64-unknown-linux-gnu
PGM = v.kriging
# other PCL version:
#-lpcl_apps -lpcl_common -lpcl_features -lpcl_filters -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_visualization
-LIBES = $(RASTER3DLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GMATHLIB) $(GISLIB) $(IOSTREAMLIB) -lpcl_common -lpcl_features -lpcl_filters -lpcl_io_ply -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_outofcore -lpcl_people -lpcl_recognition -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_tracking -lpcl_visualization
+LIBES = $(RASTER3DLIB) $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GMATHLIB) $(GISLIB) $(IOSTREAMLIB) $(RTREELIB) -lpcl_common -lpcl_features -lpcl_filters -lpcl_io_ply -lpcl_io -lpcl_kdtree -lpcl_keypoints -lpcl_octree -lpcl_outofcore -lpcl_people -lpcl_recognition -lpcl_registration -lpcl_sample_consensus -lpcl_search -lpcl_segmentation -lpcl_surface -lpcl_tracking -lpcl_visualization
-DEPENDENCIES= $(RASTER3DDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GMATHDEP) $(GISDEP) $(IOSTREAMDEP)
+DEPENDENCIES= $(RASTER3DDEP) $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GMATHDEP) $(GISDEP) $(IOSTREAMDEP) $(RTREEDEP)
EXTRA_INC = $(VECT_INC) -I/usr/include/pcl-1.7 -I/usr/include/eigen3 $(PCL_INCLUDE_DIRS) $(PROJINC)
EXTRA_CFLAGS = $(VECT_CFLAGS)
Modified: grass-addons/grass7/vector/v.kriging/local_proto.h
===================================================================
--- grass-addons/grass7/vector/v.kriging/local_proto.h 2015-06-20 09:32:39 UTC (rev 65500)
+++ grass-addons/grass7/vector/v.kriging/local_proto.h 2015-06-20 11:23:47 UTC (rev 65501)
@@ -3,16 +3,11 @@
#include <stdio.h>
#include <stdlib.h>
+#include <unistd.h>
#include <string.h>
#include <math.h>
#include <time.h>
-#include <pcl/point_types.h>
-#include <pcl/io/pcd_io.h>
-#include <pcl/point_cloud.h>
-#include <pcl/kdtree/impl/kdtree_flann.hpp>
-
-extern "C" {
#include <grass/vector.h>
#include <grass/dbmi.h>
#include <grass/config.h>
@@ -21,10 +16,12 @@
#include <grass/la.h>
#include <grass/raster3d.h>
#include <grass/glocale.h>
-}
+#include <grass/rtree.h>
#ifndef PI
#define PI M_PI
+#define DEG2RAD(ang) (ang / 180. * PI)
+#define RAD2DEG(ang) (ang / PI * 180.)
#endif
#ifndef SQUARE
@@ -35,7 +32,7 @@
struct opts
{
- struct Option *input, *output, *phase, *report, *crossvalid, *function_var_hz, *function_var_vert, *function_var_final, *form_file, *field, *intpl, *zcol, *var_dir_hz, *var_dir_vert, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert, *nugget_final, *sill_hz, *sill_vert, *sill_final, *range_hz, *range_vert, *range_final;
+ struct Option *input, *output, *phase, *report, *crossvalid, *function_var_hz, *function_var_vert, *function_var_final, *function_var_final_vert, *form_file, *field, *intpl, *zcol, *var_dir_hz, *var_dir_vert, *nL, *nZ, *td_hz, *td_vert, *nugget_hz, *nugget_vert, *nugget_final, *nugget_final_vert, *sill_hz, *sill_vert, *sill_final, *sill_final_vert, *range_hz, *range_vert, *range_final, *range_final_vert;
};
struct flgs
@@ -48,40 +45,38 @@
struct select
{
- int n; // # of selected
- int out; // # of the others
- int total; // total number
+ int n; // # of selected
+ int out; // # of the others
+ int total; // total number
int *indices; // indices of selected
};
-struct points // inputs
+struct points // inputs
{
- int n; // number of points
- double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
- double *r_min; // min coords
- double *r_max; // max coords
- double center[3]; // center
- double *invals; // values to be interpolated
- struct select in_reg; // points in region
+ int n; // number of points
+ double *r; // triples of coordinates (e.g. x0 y0 z0... xn yn zn)
+ struct RTree *R_tree; // spatial index
+ struct RTree *Rtree_hz; // spatial index
+ struct RTree *Rtree_vert; // spatial index
+ double *r_min; // min coords
+ double *r_max; // max coords
+ double max_dist; // maximum distance
+ double center[3]; // center
+ double *invals; // values to be interpolated
+ struct select in_reg; // points in region
mat_struct *trend;
};
-struct pcl_utils
-{
- pcl::PointCloud<pcl::PointXYZ>::Ptr pnts;
- pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree;
-
- pcl::PointCloud<pcl::PointXYZ>::Ptr pnts_hz;
- pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree_A;
-
- pcl::PointCloud<pcl::PointXYZ>::Ptr pnts_vert;
- pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr kd_tree_xy;
-};
-
struct bivar
{
int vert;
char *variogram;
+
+ double *h;
+ double max_dist;
+ int nLag;
+ double lag;
+
int function;
double sill;
double nugget;
@@ -90,41 +85,44 @@
struct parameters
{
- int function;
- int type;
- int const_val;
- double dir;
- double td; // range of directions
+ int function; // variogram function: lin, exp, spher, Gauss, bivar
+ int type; // variogram type: hz / vert / aniso / bivar
+ int const_val;
+ double dir; // azimuth for variogram computing
+ double td; // maximum azimuth
- double radius; // radius
- double max_dist; // maximum distance
- double max_dist_vert;
+ double radius; // radius (squared maximum distance)
+ double max_dist; // maximum distance - hz / aniso
+ double max_dist_vert; // maximum distance - vert
- int nLag; // number of length pieces
- double lag; // maximum distance between nearest neighbours (variogram lag)
+ int nLag; // number of lags - hz / aniso
+ double lag; // lag size
- int nLag_vert;
- double lag_vert;
+ int nLag_vert; // number of lags - vert
+ double lag_vert; // lag size
- double *h; // horizontal intervals used to estimate experimental variogram
- double *vert; // vertical intervals used to estimate experimental variogram
- mat_struct *c; // number of elements in each lag (final variogram only)
- mat_struct *gamma; // value of experimental variogram
+ double *h; // lag distance from search point - hz / aniso
+ double *vert; // lag distance from search point - vert
+
+ int gamma_n; // # of dissimilarities between input points
+ mat_struct *gamma; // experimental variogram matrix
+ double gamma_sum; // sum of gamma values
+
double nugget;
double sill;
double part_sill;
double h_range;
- struct bivar horizontal;
- struct bivar vertical;
+ struct bivar horizontal; // horizontal variogram properties
+ struct bivar vertical; // vertical variogram properties
- mat_struct *A;
- mat_struct *T; // coefficients of theoretical variogram
- mat_struct *GM; // matrix of diffences between point values based on the distances and the theoretical variogram
+ mat_struct *A; // plan matrix
+ mat_struct *T; // coefficients of theoretical variogram
+ mat_struct *GM; // GM = theor_var(dist: input, output points)
- char *name; // name of input vector layer
- char term[12]; // output format - gnuplot terminal
- char ext[4]; // output format - file format
+ char *name; // name of input vector layer
+ char term[12]; // output format - gnuplot terminal
+ char ext[4]; // output format - file format
};
struct var_par // parameters of experimental variogram
@@ -159,58 +157,70 @@
struct write crossvalid;
};
-extern "C" {
- struct reg_par // Region settings -> output extent and resolution
- {
- struct Cell_head reg_2d; // region for 2D interpolation
- RASTER3D_Region reg_3d; // region for 3D interpolation
- double west; // region.west
- double east;
- double north; // region.north
- double south;
- double bot; // region.bottom
- double top;
- double ew_res; // east-west resolution
- double ns_res; // north-south resolution
- double bt_res; // bottom-top resolution
- int nrows; // number of rows
- int ncols; // number of cols
- int ndeps; // number of deps
- int nrcd;
- };
+struct reg_par // Region settings -> output extent and resolution
+{
+ struct Cell_head reg_2d; // region for 2D interpolation
+ RASTER3D_Region reg_3d; // region for 3D interpolation
+ double west; // region.west
+ double east; // region.east
+ double north; // region.north
+ double south; // region.south
+ double bot; // region.bottom
+ double top; // region.top
+ double ew_res; // east-west resolution
+ double ns_res; // north-south resolution
+ double bt_res; // bottom-top resolution
+ int nrows; // # of rows
+ int ncols; // # of cols
+ int ndeps; // # of deps
+ int nrcd; // # of cells
+};
- struct output
- {
- DCELL *dcell;
- int fd_2d;
- RASTER3D_Map *fd_3d;
- char *name; // name of output 2D/3D raster
- };
-}
+struct output
+{
+ DCELL *dcell;
+ int fd_2d;
+ RASTER3D_Map *fd_3d;
+ char *name; // name of output 2D/3D raster
+};
double *get_col_values(struct Map_info *, struct int_par *, struct points *, int, const char *, int);
-void test_normality(int , double *, struct write *);
-void read_points(struct Map_info *, struct reg_par *, struct points *, struct pcl_utils *, struct int_par, const char *, int, struct write *);
+void test_normality(int, double *, struct write *);
+void read_points(struct Map_info *, struct reg_par *, struct points *, struct int_par *, const char *, int, struct write *);
double min(double *, struct points *);
double max(double *, struct points *);
+struct RTree *create_spatial_index(struct int_par *);
+void insert_rectangle(int, int, struct points *);
+struct ilist *find_NNs_within(int, int, struct points *, double, double);
+struct ilist *find_n_NNs(int, int, struct points *, int);
+double sum_NN(int, int, struct ilist *, struct points *);
+
+void correct_indices(int, struct ilist *, double *, struct points *, struct parameters *);
+int cmpVals(const void *, const void *);
void coord_diff(int, int, double *, double *);
double distance_diff(double *);
double radius_hz_diff(double *);
-double *triple(double, double, double);
-double lag_distance(int, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *, struct write *);
+double zenith_angle(double *);
+void triple(double, double, double, double *);
+double lag_size(int, struct int_par *, struct points *, struct parameters *, struct write *);
int lag_number(double, double *);
-void variogram_restricts(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *);
-void geometric_anisotropy(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr);
+void variogram_restricts(struct int_par *, struct points *, struct parameters *);
+void geometric_anisotropy(struct int_par *, struct points *);
double find_intersect_x(double *, double *, double *, double *, struct write *);
double find_intersect_y(double *, double *, double *, double *, double , struct write *);
mat_struct *LSM(mat_struct *, mat_struct *);
mat_struct *nonlin_LMS(int , double *, double *);
-void E_variogram(int, struct int_par *, struct points *, struct pcl_utils *, struct var_par *);
-void T_variogram(int, struct opts, struct parameters *, struct write *);
-void ordinary_kriging(struct int_par *, struct reg_par *, struct points *, struct pcl_utils *, struct var_par *, struct output *);
-int set_function(char *, struct parameters *, struct write *);
+void E_variogram(int, struct int_par *, struct points *, struct reg_par *, struct var_par *);
+void T_variogram(int, int, struct opts, struct parameters *, struct write *);
+void ordinary_kriging(struct int_par *, struct reg_par *, struct points *, struct var_par *, struct output *);
+
+void linear_variogram(struct parameters *, struct write *);
+double bivar_sill(int, mat_struct *);
+void sill(struct parameters *);
+int sill_compare(struct int_par *, struct flgs *, struct var_par *, struct points *);
+int set_function(char *, struct write *);
double RBF(double *);
double linear(double, double, double);
double exponential(double, double, double, double);
@@ -218,32 +228,30 @@
double gaussian(double, double, double, double);
double variogram_fction(struct parameters *, double *);
void set_gnuplot(char *, struct parameters *);
-void plot_experimental_variogram(struct int_par *, double *, mat_struct *, struct parameters *);
-void plot_var(int, double *, mat_struct *, struct parameters *);
+void plot_experimental_variogram(struct int_par *, struct parameters *);
+void plot_var(int, int, struct parameters *);
void variogram_type(int, char *);
void write2file_basics(struct int_par *, struct opts *);
void write2file_vector(struct int_par *, struct points *);
void write2file_values(struct write *, const char *);
void write2file_varSetsIntro(int, struct write *);
void write2file_varSets(struct write *, struct parameters *);
-void write2file_variogram_E(struct int_par *, struct parameters *);
+void write2file_variogram_E(struct int_par *, struct parameters *, mat_struct *);
void write2file_variogram_T(struct write *);
-void write_temporary2file(struct int_par *, struct parameters *, mat_struct *);
+void write_temporary2file(struct int_par *, struct parameters *);
void read_tmp_vals(const char *, struct parameters *, struct int_par *);
mat_struct *set_up_G(struct points *, struct parameters *, struct write *);
-mat_struct *set_up_g0(struct int_par *, std::vector<int>, struct points *, double *, struct parameters *);
-mat_struct *submatrix(std::vector<int> , mat_struct *, struct write *);
-double result(std::vector<int>, struct points *, mat_struct *);
+mat_struct *set_up_g0(struct int_par *, struct points *, struct ilist *, double *, struct parameters *);
+mat_struct *submatrix(struct ilist *, mat_struct *, struct write *);
+double result(struct points *, struct ilist *, mat_struct *);
-void crossvalidation(struct int_par *, struct points *, pcl::PointCloud<pcl::PointXYZ>::Ptr, struct parameters *);
+void crossvalidation(struct int_par *, struct points *, struct parameters *);
void cell_centre(unsigned int, unsigned int, unsigned int, struct int_par *, struct reg_par *, double *, struct parameters *);
-extern "C" {
void get_region_pars(struct int_par *, struct reg_par *);
void open_layer(struct int_par *, struct reg_par *, struct output *);
int write2layer(struct int_par *, struct reg_par *, struct output *, unsigned int, unsigned int, unsigned int, double);
-}
static inline double get_quantile(int);
static void get_slot_counts(int, double *);
@@ -251,6 +259,6 @@
static void fill_bins(int, double *);
static int compare(const void *, const void *);
static void sort_bins(void);
-static void compute_quantiles(int, double, struct write *);
+static void compute_quantiles(int, double *, struct write *);
double quantile(double, int, double *, struct write *);
#endif
Added: grass-addons/grass7/vector/v.kriging/main.c
===================================================================
--- grass-addons/grass7/vector/v.kriging/main.c (rev 0)
+++ grass-addons/grass7/vector/v.kriging/main.c 2015-06-20 11:23:47 UTC (rev 65501)
@@ -0,0 +1,594 @@
+
+/****************************************************************
+ *
+ * MODULE: v.kriging
+ * AUTHOR: Eva Stopková
+ * in case of functions taken from another modules,
+ * they are cited above the function
+ * or at the beginning of the file (e.g. quantile.cpp
+ * that uses slightly modified functions
+ * taken from the module r.quantile (Clemens, G.))
+ * PURPOSE: Module interpolates the values to two- or three-dimensional grid using input values
+ * located on 2D/3D point vector layer. Interpolation method
+ * Ordinary kriging has been extended for 3D points (v = f(x,y) -> v = f(x,y,z)).
+ *
+ * COPYRIGHT: (C) 2012-2014 Eva Stopková and by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ **************************************************************/
+#include "local_proto.h"
+
+
+int main(int argc, char *argv[])
+{
+ // Vector layer and module
+ struct Map_info map; // Input vector map
+ struct GModule *module; // Module
+
+ struct reg_par reg; // Region parameters
+ struct points pnts; // Input points (coordinates, extent, values, etc.)
+
+ // Geostatistical parameters
+ struct int_par xD; // 2D/3D interpolation for 2D/3D vector layer
+ struct var_par var_pars; // Variogram (experimental and theoretical)
+
+ // Outputs
+ struct output out; // Output layer properties
+ FILE *fp;
+
+ // Settings
+ int field;
+ struct opts opt;
+ struct flgs flg;
+
+ /* ------- Module creation ------- */
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("3D raster"));
+ G_add_keyword(_("ordinary kriging - for 2D and 3D data"));
+ module->description =
+ _("Interpolates 2D or 3D raster based on input values located on 2D or 3D point vector layer (method ordinary kriging extended to 3D).");
+
+ // Setting options
+ opt.input = G_define_standard_option(G_OPT_V_INPUT); // Vector input layer
+ opt.input->label = _("Name of input vector points map");
+
+ flg.d23 = G_define_flag(); // Option to process 2D or 3D interpolation
+ flg.d23->key = '2';
+ flg.d23->description = _("Force 2D interpolation even if input is 3D");
+ flg.d23->guisection = _("3D");
+
+ opt.field = G_define_standard_option(G_OPT_V_FIELD);
+
+ opt.phase = G_define_option();
+ opt.phase->key = "phase";
+ opt.phase->options = "initial, middle, final";
+ opt.phase->description = _("Phase of interpolation. 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.output = G_define_option(); // Output layer
+ opt.output->key = "output";
+ opt.output->description =
+ _("Name for output 2D/3D raster map");
+ opt.output->guisection = _("Final");
+
+ opt.report = G_define_standard_option(G_OPT_F_OUTPUT); // Report file
+ opt.report->key = "report";
+ opt.report->description = _("File to write the report");
+ opt.report->required = NO;
+ opt.report->guisection = _("Initial");
+
+ 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");
+
+ flg.bivariate = G_define_flag();
+ flg.bivariate->key = 'b';
+ flg.bivariate->description = _("Compute bivariate variogram (3D interpolation only)");
+ flg.bivariate->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_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_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");
+
+ 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_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");
+
+ flg.detrend = G_define_flag();
+ flg.detrend->key = 't';
+ flg.detrend->description = _("Eliminate trend if variogram is parabolic");
+ flg.detrend->guisection = _("Initial");
+
+ opt.form_file = G_define_option(); // Variogram plot - various output formats
+ opt.form_file->key = "fileformat";
+ opt.form_file->options = "cdr,dxf,eps,tex,pdf,png,svg";
+ opt.form_file->description = _("File format to save variogram plot (empty: preview in Gnuplot terminal)");
+ opt.form_file->guisection = _("Middle");
+
+ opt.intpl = G_define_standard_option(G_OPT_DB_COLUMN); // Input values for interpolation
+ opt.intpl->key = "icolumn";
+ opt.intpl->description =
+ _("Attribute column containing input values for interpolation");
+ opt.intpl->required = YES;
+
+ opt.zcol = G_define_standard_option(G_OPT_DB_COLUMN); // Column with z coord (2D points)
+ opt.zcol->key = "zcolumn";
+ opt.zcol->description =
+ _("Attribute column containing z coordinates (only for 3D interpolation based on 2D point layer)");
+ opt.zcol->required = NO;
+ opt.zcol->guisection = _("3D");
+
+ opt.var_dir_hz = G_define_option();
+ opt.var_dir_hz->key = "azimuth";
+ opt.var_dir_hz->type = TYPE_DOUBLE;
+ opt.var_dir_hz->required = NO;
+ opt.var_dir_hz->answer = "0.0";
+ opt.var_dir_hz->description =
+ _("Azimuth of variogram computing (isotrophic)");
+ opt.var_dir_hz->guisection = _("Initial");
+
+ opt.var_dir_vert = G_define_option();
+ opt.var_dir_vert->key = "zenith_angle";
+ opt.var_dir_vert->type = TYPE_DOUBLE;
+ opt.var_dir_vert->required = NO;
+ opt.var_dir_vert->answer = "0.0";
+ opt.var_dir_vert->description =
+ _("Zenith angle of variogram computing (isotrophic)");
+ opt.var_dir_vert->guisection = _("Initial");
+
+ opt.nL = G_define_option();
+ opt.nL->key = "lpieces";
+ opt.nL->type = TYPE_INTEGER;
+ opt.nL->required = NO;
+ opt.nL->description = _("Number of horizontal lags");
+ opt.nL->guisection = _("Initial");
+
+ opt.nZ = G_define_option();
+ opt.nZ->key = "vpieces";
+ opt.nZ->type = TYPE_INTEGER;
+ opt.nZ->required = NO;
+ opt.nZ->description =
+ _("Number of vertical lags (only for 3D variogram)");
+ opt.nZ->guisection = _("Initial");
+
+ opt.td_hz = G_define_option();
+ opt.td_hz->key = "td";
+ opt.td_hz->type = TYPE_DOUBLE;
+ opt.td_hz->answer = "90.0";
+ opt.td_hz->description = _("Angle of variogram processing");
+ opt.td_hz->required = NO;
+ opt.td_hz->guisection = _("Initial");
+
+ opt.nugget_hz = G_define_option();
+ opt.nugget_hz->key = "hz_nugget";
+ opt.nugget_hz->type = TYPE_DOUBLE;
+ opt.nugget_hz->description = _("Nugget effect of horizontal variogram");
+ opt.nugget_hz->required = NO;
+ opt.nugget_hz->guisection = _("Middle");
+
+ opt.nugget_vert = G_define_option();
+ opt.nugget_vert->key = "vert_nugget";
+ opt.nugget_vert->type = TYPE_DOUBLE;
+ opt.nugget_vert->description = _("Nugget effect of vertical variogram");
+ opt.nugget_vert->required = NO;
+ opt.nugget_vert->guisection = _("Middle");
+
+ opt.nugget_final = G_define_option();
+ opt.nugget_final->key = "final_nugget";
+ opt.nugget_final->type = TYPE_DOUBLE;
+ opt.nugget_final->multiple = TRUE;
+ 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_final_vert = G_define_option();
+ opt.nugget_final_vert->key = "final_vert_nugget";
+ opt.nugget_final_vert->type = TYPE_DOUBLE;
+ 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_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_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.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_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_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_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.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_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");
+ /* --------------------------------------------------------- */
+
+ G_gisinit(argv[0]);
+
+ if (G_parser(argc, argv)) {
+ exit(EXIT_FAILURE);
+ }
+
+ /* Get parameters from the parser */
+ if (strcmp(opt.phase->answer, "initial") == 0) {
+ xD.phase = 0; // estimate range
+ }
+ else if (strcmp(opt.phase->answer, "middle") == 0) {
+ xD.phase = 1; // estimate anisotropic variogram
+ }
+ else if (strcmp(opt.phase->answer, "final") == 0) {
+ xD.phase = 2; // compute kriging
+ }
+
+ // Open report file if desired
+ if (opt.report->answer) {
+ xD.report.write2file = TRUE;
+ xD.report.name = opt.report->answer;
+ xD.report.fp = fopen(xD.report.name, "w");
+ time(&xD.report.now);
+ fprintf(xD.report.fp, "v.kriging started on %s\n\n", ctime(&xD.report.now));
+ G_message(_("Report is being written to %s..."), xD.report.name);
+ }
+ else
+ xD.report.write2file = FALSE;
+
+ if (opt.crossvalid->answer) {
+ xD.crossvalid.write2file = TRUE;
+ xD.crossvalid.name = opt.crossvalid->answer;
+ xD.crossvalid.fp = fopen(xD.crossvalid.name, "w");
+ }
+ else {
+ xD.crossvalid.write2file = FALSE;
+ }
+
+ var_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);
+ }
+ G_fatal_error(_("You should mark either univariate, or bivariate variogram, not both of them..."));
+ } // error
+
+ /* ---------------------------------------------------------- */
+ Vect_set_open_level(2); // Open input vector map
+
+ if (0 > Vect_open_old2(&map, opt.input->answer, "", opt.field->answer)) {
+ if (xD.report.write2file == TRUE) {
+ fclose(xD.report.fp);
+ remove(xD.report.name);
+ }
+ G_fatal_error(_("Unable to open vector map <%s>"), opt.input->answer);
+ } // error
+ Vect_set_error_handler_io(&map, NULL);
+ /* ---------------------------------------------------------- */
+
+ /* Perform 2D or 3D interpolation ? */
+ xD.i3 = flg.d23->answer ? FALSE : TRUE; // 2D/3D interpolation
+ xD.v3 = Vect_is_3d(&map); // 2D/3D layer
+
+ /* What could happen:
+ *-------------------
+ * 3D interpolation + 3D points = 3D GIS
+ * 3D interpolation + 2D points = 2,5D -> 3D GIS (needs attribute column with z and 3D region)
+ * 2D interpolation + 3D points = 3D -> 2,5D GIS
+ * 2D interpolation + 2D points = 2,5D GIS */
+
+ // 3D interpolation
+ if (xD.i3 == TRUE) { // 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..."));
+ }
+ }
+
+ field = Vect_get_field_number(&map, opt.field->answer);
+ if (xD.report.write2file == TRUE)
+ write2file_basics(&xD, &opt);
+ /* ---------------------------------------------------------- */
+
+ /* Get... */
+ get_region_pars(&xD, ®); // ... region parameters
+ read_points(&map, ®, &pnts, &xD, opt.zcol->answer, field, &xD.report); // ... coordinates of points
+ pnts.invals = get_col_values(&map, &xD, &pnts, field, opt.intpl->answer, flg.detrend->answer); // ... values for interpolation
+
+ /* Estimate 2D/3D variogram */
+ switch (xD.phase) {
+ case 0: // initial phase
+ var_pars.hz.type = 0; // horizontal variogram
+ if (xD.i3 == TRUE) { // 3D interpolation:
+ var_pars.vert.type = 1; // vertical variogram
+ }
+
+ 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
+ }
+ }
+ 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;
+
+ 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 (xD.report.name) { // report file available:
+ xD.report.write2file = TRUE;
+ xD.report.fp = fopen(xD.report.name, "a");
+ }
+
+ 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
+
+ /* 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) { // anisotropic 3D:
+ if (opt.function_var_final_vert->answer || opt.nugget_final_vert->answer || 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) { // nonlinear:
+ 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 of horizontal component of bivariate variogram..."));
+ } // 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 of horizontal component of bivariate variogram..."));
+ } // end if range missing
+ }
+ } // 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
+ break;
+ }
+
+ /* Ordinary kriging (including 2D/3D raster creation) */
+ ordinary_kriging(&xD, ®, &pnts, &var_pars, &out);
+ /* ---------------------------------------------------------- */
+
+ end:
+
+ G_free(pnts.r);
+ Vect_close(&map); // Close vector map
+ exit(EXIT_SUCCESS);
+}
+
Modified: grass-addons/grass7/vector/v.kriging/v.kriging.html
===================================================================
--- grass-addons/grass7/vector/v.kriging/v.kriging.html 2015-06-20 09:32:39 UTC (rev 65500)
+++ grass-addons/grass7/vector/v.kriging/v.kriging.html 2015-06-20 11:23:47 UTC (rev 65501)
@@ -1,7 +1,7 @@
<h2>DESCRIPTION</h2>
<em>v.kriging</em> interpolates unknown values using method <i>ordinary
-kriging</i>. Output can be 2D (temporarily denied) or 3D.
+kriging</i>. Output can be 2D or 3D.
<h2>EXAMPLES</h2> Input layer should contain 3D coordinates (xyz) and
values to be interpolated (in attribute table). The commands can look
@@ -19,6 +19,21 @@
out=name crossval=crossval_file.txt
</pre></div>
+<h2>TODO</h2>
+<ul>
+<li><b>anisotropy</b> in horizontal direction missing
+<li>current version is suitable just for <b>metric coordinate systems</b>
+<li>enable <b>mask usage</b>
+<li><b>bivariate variogram</b> needs to be rebuilt (theory)
+<li><b>2D interpolation from 3D input layer</b> needs to be rebuilt (especially in case that there are too many points located on identical horizontal coordinates with different elevation)
+</ul>
+
+<h2>BUGS</h2>
+<ul>
+<li>2D variogram gives quite inaccurate results for <b>sparse (or spatially heterogeneous) data</b>
+<li><b>bivariate variogram</b> - too inaccurate results
+</ul>
+
<h2>SEE ALSO</h2>
<em>
@@ -33,7 +48,6 @@
GMATH library (GRASS Numerical Library)<br>
<a href="http://www.netlib.org/lapack">http://www.netlib.org/lapack</a>
(usually available on Linux distros)
-<li><a href="http://pointclouds.org">Point Cloud Library</a>
</ul>
<h2>AUTHOR</h2>
More information about the grass-commit
mailing list