/*
Gstat, a program for geostatistical modelling, prediction and simulation
Copyright 1992, 2011 (C) Edzer Pebesma
Edzer Pebesma, edzer.pebesma@uni-muenster.de
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. As a special exception, linking
this program with the Qt library is permitted.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
(read also the files COPYING and Copyright)
*/
/*
* plot.c: kriging weights/variogram plot drivers to gnuplot and jgraph
*/
#include <stdio.h>
#include <stdlib.h> /* getenv() */
#include <string.h>
#include <math.h>
#include "defs.h"
#include "userio.h"
#include "debug.h"
#include "data.h"
#include "utils.h"
#include "vario.h"
#include "glvars.h"
#include "sem.h"
#include "version.h"
#include "plot.h"
/* jgraph defaults: */
#define SIZE 0.013 /* rel distance (in) of number from cross */
#define N_PTS 200 /* number of semivariogram model points for jgraph */
#define PLOT_WINDOW_NR(v) \
(v->id1==v->id2?v->id1:(get_n_vars()+LTI2(v->id1,v->id2)))
typedef struct {
PLOT_TYPE p;
int pt, lt;
char *term_name, *term_opt, *gamma_str, *distance_str;
} GNUPLOT_TERM;
GNUPLOT_TERM gnuplot_terms[] = {
{ GNUPLOT, 2, 3, NULL, NULL, NULL, NULL }, /* use default option */
{ PSLATEX, 1, 1, "pslatex", "", "$\\gamma(h)$", "$h$" },
{ CGM, 1, 1, "cgm", "", NULL },
{ GIF, 2, 2, "gif",
"transparent size 640, 480 xffffff x000000 x404040 xff0000 x0000ff x00ff00",
NULL, NULL },
{ EPS, 1, 1, "postscript", "eps solid 17", NULL, NULL },
{ PNG, 2, 3,
"png", "small transparent xffffff x000000 x404040 xff0000",
NULL, NULL },
{ EEPIC, 2, 3, "eepic", "", "{$\\gamma(h)$}", "{$h$}" },
{ UNKNOWN, 0, 0, NULL, NULL, NULL, NULL }
};
static int set_key(FILE *f, const VARIOGRAM *v, double min, double max);
static void get_minmax_y(const VARIOGRAM *v, double max_x, double *min_y, double *max_y);
#ifndef USING_R
int fprint_gnuplot_variogram(FILE *stream, const VARIOGRAM *v,
char *fname, PLOT_TYPE p, int window_nr) {
double min_y = 0.0, max_y = 0.0, max_x = 0.0, range;
int i, key_outside = 1;
GNUPLOT_TERM t = gnuplot_terms[0];
const char *est = NULL;
int inline_data = 0;
assert(v->ev != NULL);
assert(stream != NULL);
/* plot sample vgm and/or model: */
assert(v->ev->n_est > 0 || v->n_models > 0);
for (i = 0; gnuplot_terms[i].p != UNKNOWN; i++) /* find the entry */
if (gnuplot_terms[i].p == p) {
t = gnuplot_terms[i];
break;
}
assert(t.p != UNKNOWN);
if (gl_gpterm != NULL) {
t.term_name = gl_gpterm;
t.term_opt = "";
}
fprintf(stream, "#\n# gnuplot file, created by gstat %s\n#\n\n", VERSION);
fprintf(stream, "set tics out\n");
if (fname) {
fprintf(stream, "## for %s term:\n", t.term_name);
if (t.term_name != NULL)
fprintf(stream, "set term %s %s\n", t.term_name, t.term_opt);
if (p == PSLATEX || p == EPS || p == EEPIC) {
fprintf(stream, "set border 3\n");
fprintf(stream, "set xtics nomirror\nset ytics nomirror\n");
}
/* output file: */
if (*fname != '\0')
fprintf(stream, "%sset output '%s'\n\n",
p == GNUPLOT ? "# " : "", fname);
} else { /* plot to screen terminal */
#ifdef WIN32
fprintf(stream, "set term windows\n");
#else
if (gl_gnuplot35 == NULL && getenv("DISPLAY") != NULL
&& getenv("GNUTERM") == NULL)
/* will default to x11; set nr for 3.7: */
fprintf(stream, "set term x11 %d\n",
window_nr >= 0 ? window_nr : PLOT_WINDOW_NR(v));
#endif
}
/* print basic models: */
for (i = 1; v_models[i].name != NULL; i++)
fprintf(stream, "%s\n",
is_variogram(v) ? v_models[i].v_gnuplot : v_models[i].c_gnuplot);
fprintf(stream, "\n# erase previous label settings:\nset nolabel\n");
if (p != EEPIC)
fprintf(stream, "\n# set dottet lines at y=0:\nset zeroaxis\n");
fprintf(stream, "\n# set x and y labels:\n");
fprintf(stream, "set xlabel '%s'\n",
t.distance_str == NULL ? "distance" : t.distance_str);
fprintf(stream, "set ylabel '");
switch (v->ev->evt) {
case SEMIVARIOGRAM:
fprintf(stream, "%s",
t.gamma_str == NULL ? "semivariance" : t.gamma_str);
break;
case CROSSVARIOGRAM:
if (v->ev->pseudo)
fprintf(stream, "%s", "pseudo cross semivariance");
else
fprintf(stream, "%s", "cross semivariance");
break;
case COVARIOGRAM:
fprintf(stream, "%s", "C(h)");
break;
case CROSSCOVARIOGRAM:
fprintf(stream, "%s", "cross covariance");
break;
case NOTSPECIFIED:
ErrMsg(ER_IMPOSVAL, "evt unspecified");
break;
}
if (v->ev->cloud)
fprintf(stream, " cloud");
if (v->ev->is_directional) {
fprintf(stream, " (dir. ");
if (gl_tol_hor < 90) {
fprintf(stream, "<x,y> %g +/- %g", gl_alpha, gl_tol_hor);
if (gl_tol_ver < 90)
fprintf(stream, "; ");
}
if (gl_tol_ver < 90)
fprintf(stream, "<z> %g +/- %g", gl_beta, gl_tol_ver);
fprintf(stream, ")");
}
fprintf(stream, "'\n\n");
if (v->ev->plot_numbers)
fprintf(stream, "\n# set number of point pairs as labels:\n");
for (i = 0; i < (v->ev->zero == ZERO_AVOID ?
v->ev->n_est - 1 : v->ev->n_est); i++) {
if (v->ev->nh[i] > 0) {
if (v->ev->plot_numbers && v->ev->n_est <= gl_dots) {
if (v->ev->cloud)
fprintf(stream, "set label '%ld,%ld' at %e,%e left\n",
HIGH_NH(v->ev->nh[i]) + 1, LOW_NH(v->ev->nh[i]) + 1,
v->ev->dist[i] + SIZE * v->ev->cutoff, v->ev->gamma[i]);
else
fprintf(stream, "set label '%ld' at %e,%e left\n",
v->ev->nh[i], v->ev->dist[i] + SIZE * v->ev->cutoff,
v->ev->gamma[i]);
}
}
}
if (v->ev->n_est == 0)
max_x = 2.5 * effective_range(v);
else
max_x = 1.05 * v->ev->cutoff;
get_minmax_y(v, max_x, &min_y, &max_y);
if (gl_gnuplot35 == NULL)
key_outside = set_key(stream, v, min_y, max_y);
fprintf(stream, "\n# set axis limits:\n");
fprintf(stream, "xmin = 0; xmax = %e\n", max_x);
fprintf(stream, "set xrange [xmin:xmax]\n");
if (max_y > min_y) {
range = max_y - min_y;
fprintf(stream, "ymin = %e; ymax = %e\n",
min_y < 0.0 ? min_y - 0.05 * range : 0.0,
max_y + 0.05 * range + 0.15 * key_outside * range);
fprintf(stream, "set yrange [ymin:ymax]\n");
}
/* plot file: */
fprintf(stream, "\n# do the data & function plotting:\n");
fprintf(stream, "plot ");
if (v->ev->n_est > 0) {
if (v->fname == NULL && o_filename == NULL) {
inline_data = 1;
est = "-";
} else
est = v->fname ? v->fname : o_filename;
fprintf(stream, "'%s'", est);
/* using: */
if (v->ev->cloud)
fprintf(stream, "%s", " using 3:4 ");
else
fprintf(stream, "%s", " using 4:5 ");
/* title: */
if (v->id1 == v->id2)
fprintf(stream, "title '%s' ", name_identifier(v->id1));
else
fprintf(stream, "title '%s x %s' ", name_identifier(v->id1),
name_identifier(v->id2));
if (v->ev->n_est > gl_dots)
fprintf(stream, "with dots 1");
else
fprintf(stream, "with points pt %d", t.pt);
}
if (v->n_models > 0) {
fprintf(stream, "%s\\\n ", v->ev->n_est > 0 ? "," : "");
fprint_gnuplot_model(stream, v, 0);
fprintf(stream, "\\\n title '%s' with lines lt %d", v->descr, t.lt);
/*
if (p == PSLATEX)
fprintf(stream, " with lines 1 1");
*/
}
fprintf(stream, "\n"); /* do it ! */
if (inline_data) {
fprint_sample_vgm(stream, v->ev);
fprintf(stream, "e\n");
}
#ifdef WIN32
/* prevent wgnuplot from disappearing; kms */
fprintf(stream, "pause -1 \"Hit Enter to continue\"\n");
#endif
fflush(stream);
return 0;
}
void fprint_gnuplot_model(FILE *f, const VARIOGRAM *vgm, int fit) {
int i;
char a[50], c[50];
for (i = 0; i < vgm->n_models; i++) {
if (fit && vgm->part[i].fit_sill) /* print parameter to string: */
sprintf(c, "c%d", i);
else
sprintf(c, "%.2e", vgm->part[i].sill);
if (fit && vgm->part[i].fit_range) /* print parameter to string: */
sprintf(a, "a%d", i);
else
sprintf(a, "%.2e", vgm->part[i].range[0] *
relative_norm(vgm->part[i].tm_range, vgm->ev->direction.x,
vgm->ev->direction.y, vgm->ev->direction.z));
fprintf(f, "%s * %s(%s,x)", c,
v_models[vgm->part[i].model].name, a);
if (i < vgm->n_models - 1)
fprintf(f, " + ");
}
}
int fprint_jgraph_variogram(FILE *f, const VARIOGRAM *v) {
/*
* print variogram to a jgraph input file
* (jgraph is a program that converts such a file to PS/EPS)
* URL: file://ftp.princeton.edu/pub/jgraph.Z
*
* don't play tricks: f might be a pipe.
*/
int i, sill_reached = 0;
double max_y = 0.0, min_y = 0.0, max_x = 0.0, x, y, z, s;
if (v->ev->n_est <= 0)
return 1;
fprintf(f, "(*\n jgraph file, created by gstat %s\n*)\n",
VERSION);
fprintf(f, "newgraph\n");
fprintf(f, "(* draw estimates as crosses: *)\n");
fprintf(f, "newcurve marktype cross label : ");
if (v->id1 == v->id2)
fprintf(f, "%s\n pts\n", name_identifier(v->id1));
else
fprintf(f, "%s x %s\n pts\n", name_identifier(v->id1),
name_identifier(v->id2));
for (i = 0; i < ((v->ev->zero == ZERO_AVOID) ?
v->ev->n_est - 1 : v->ev->n_est); i++)
if (v->ev->nh[i] > 0)
fprintf(f, "%g %g\n", v->ev->dist[i], v->ev->gamma[i]);
if (v->ev->plot_numbers) {
fprintf(f, "(* draw numbers *)\n");
for (i = 0; i < v->ev->n_est; i++) {
if (v->ev->nh[i] > 0)
fprintf(f, "newstring x %g y %g hjl vjc fontsize 8 : %ld\n",
v->ev->dist[i] + SIZE * v->ev->cutoff,
v->ev->gamma[i], v->ev->nh[i]);
} /* for */
max_x = (1.0 + 2.0 * SIZE) * v->ev->cutoff;
} else
max_x = v->ev->cutoff;
get_minmax_y(v, max_x, &min_y, &max_y);
/*
* LEGEND, TITLES:
*/
fprintf(f, "(* draw legend, title and axis *)\nlegend top\n");
fprintf(f, "(* title : (fill in title) *)\n");
/*
* XAXIS:
*/
fprintf(f, "xaxis min 0 max %g size 4.0 label : Distance\n", max_x);
if (v->ev->cutoff > 99999.0 && v->ev->cutoff < 100001.0)
fprintf(f, "hash 25000 mhash 4\n");
if (v->ev->plot_numbers)
fprintf(f, "(* Note: remove `max ...' if not correct *)\n");
/*
* YAXIS:
*/
min_y *= 1.05;
max_y *= 1.05;
if (min_y == max_y) {
min_y = -1.0;
max_y = 1.0;
}
fprintf(f, "yaxis min %g ", min_y);
if (max_y > 0)
fprintf(f, "max %g ", max_y);
else
fprintf(f, "(* max 0 *) hash_format g ");
fprintf(f, "size 2.5 label : ");
switch (v->ev->evt) {
case SEMIVARIOGRAM:
fprintf(f, "Semivariance\n");
break;
case CROSSVARIOGRAM:
if (v->ev->pseudo)
fprintf(f, "Pseudo Cross Semivariance\n");
else
fprintf(f, "Cross Semivariance\n");
break;
case COVARIOGRAM:
fprintf(f, "Covariance\n");
break;
case CROSSCOVARIOGRAM:
fprintf(f, "Cross Covariance\n");
break;
case NOTSPECIFIED:
ErrMsg(ER_IMPOSVAL, "evt unspecified");
break;
}
/*
* process model:
*/
if (v->n_models > 0) {
fprintf(f, "(* draw variogram model as line: *)\n");
fprintf(f, "newline label : %s\npts\n", v->descr);
for (i = 0, x = 0; i < N_PTS; i++) {
if (sill_reached)
i = N_PTS - 1; /* jump to last one */
x = (i * v->ev->direction.x) / (N_PTS-1.0) * max_x;
y = (i * v->ev->direction.y) / (N_PTS-1.0) * max_x;
z = (i * v->ev->direction.z) / (N_PTS-1.0) * max_x;
if (x < 1.0e-30 && y < 1.0e-30 && z < 1.0e-30)
x = 1.0e-30; /* avoid line from 0 to nugget */
s = (is_variogram(v) ?
get_semivariance(v, x, y, z) :
get_covariance(v, x, y, z));
fprintf(f, "\t%g %g\n", transform_norm(NULL, x, y, z), s);
if (x > v->max_range)
sill_reached = 1;
} /* for i */
}
return 0;
}
#endif
static void get_minmax_y(const VARIOGRAM *v, double max_x, double *min_y, double *max_y) {
int i;
double x, y, z, s;
for (i = 0; i < v->ev->n_est; i++) {
if (v->ev->nh[i] > 0) {
*max_y = MAX(*max_y, v->ev->gamma[i]); /* max_y >= 0 */
*min_y = MIN(*min_y, v->ev->gamma[i]); /* min_y <= 0 */
}
}
if (v->n_models > 0) { /* get max_y, min_y: */
for (i = 0, x = 0; i < N_PTS; i++) {
x = (i * v->ev->direction.x) / (N_PTS-1.0) * max_x;
y = (i * v->ev->direction.y) / (N_PTS-1.0) * max_x;
z = (i * v->ev->direction.z) / (N_PTS-1.0) * max_x;
if (x < 1.0e-30 && y < 1.0e-30 && z < 1.0e-30)
x = 1.0e-30; /* avoid line from 0 to nugget */
s = (is_variogram(v) ?
get_semivariance(v, x, y, z) :
get_covariance(v, x, y, z));
*max_y = MAX(*max_y, s);
*min_y = MIN(*min_y, s);
}
}
}
static int set_key(FILE *f, const VARIOGRAM *v, double min, double max) {
/*
* returns 1 if key should be set outside plot region
* returns 0 if it fits inside plot region (and sets key top or bottom)
*/
int i;
enum { POS_INIT, POS_TOP, POS_BOTTOM, POS_OUTSIDE } pos = POS_INIT;
min = min + 0.25 * (max - min);
max = max - 0.25 * (max - min);
for (i = 0; i < ((v->ev->zero == ZERO_AVOID) ?
v->ev->n_est - 1 : v->ev->n_est); i++) {
if (v->ev->nh[i] > 0 && v->ev->dist[i] > 0.5 * v->ev->cutoff) {
if (pos == POS_INIT) {
pos = POS_OUTSIDE;
if (v->ev->gamma[i] < max)
pos = POS_TOP;
if (v->ev->gamma[i] > min)
pos = POS_BOTTOM;
} else {
if ((pos == POS_TOP && v->ev->gamma[i] > max) ||
(pos == POS_BOTTOM && v->ev->gamma[i] < min))
pos = POS_OUTSIDE;
}
}
}
switch (pos) {
case POS_INIT:
case POS_OUTSIDE:
fprintf(f, "\n# set key outside plotting region:\nset key\n");
return 1;
case POS_TOP:
fprintf(f, "\n# set key in plotting region:\nset key top\n");
break;
case POS_BOTTOM:
fprintf(f, "\n# set key in plotting region:\nset key bottom\n");
break;
}
return 0;
}