help with r.neighbors modification
Lars Schylberg
larss at fmi.kth.se
Wed Jul 12 08:00:00 EDT 1995
Hi, I was working with similar thing a couple of years ago.
What I did, was to write a new function for r.mapcalc. I wrote a mode
function that will be in Grass4.2. I guess you could modify it to give you
the second most frequent value and also the third most frequent value.
It would be fun to see more mapcalc funtions around.
I have included the new fucntion in the file xmode.c, then function.h and
Gmakefile have to be updated as I show below. Have fun !
Lars
Lars Schylberg Email: lasc at celsiustech.se
CelsiusTech IT or: larss at fmi.kth.se
S-175 88 Jarfalla Tel. +46 8 580 847 03
Sweden Fax. +46 8 580 123 20
===========================================================================
/* mode function contributed by
* Lars Schylberg
* Department of Geodesy and Photogrammetry,
* Royal Inst. of Technology, Stockholm, Sweden
*/
#include "gis.h"
static int icmp(i,j)
CELL *i, *j;
{
return(*i - *j);
}
static int dcmp(i,j)
double *i, *j;
{
if (*i < *j) return -1;
if (*i > *j) return 1;
return 0;
}
i_mode (argc, argv, cell, ncols)
CELL *argv[];
register CELL *cell;
register int ncols;
{
register int i, k, m, cols, max_freq_i ,max_value_i;
static int nargs = 0;
static CELL *arr = NULL;
static CELL *mfq_value = NULL;
static CELL *mfq_count = NULL;
if ( argc > nargs )
{
arr = (CELL *) G_realloc ( arr, argc*sizeof(CELL));
mfq_count = (CELL *) G_realloc ( mfq_count, argc*sizeof(CELL)+1);
mfq_value = (CELL *) G_realloc ( mfq_value, argc*sizeof(CELL)+1);
nargs = argc;
}
while (ncols-- > 0)
{
for ( i=0; i<argc; i++)
arr[i] = argv[i][ncols];
/* sort all values */
qsort( arr, argc, sizeof(CELL), icmp);
/* sort by frequency */
for ( m=0; m<argc; m++)
mfq_count[m] = 0;
mfq_value[0] = arr[0];
mfq_count[0] = 1;
i = 0;
for ( m=0; m<argc; m++)
{
if ( arr[m] == arr[m -1])
mfq_count[i]++;
else
{
i++;
mfq_value[i] = arr[m];
mfq_count[i]++;
}
}
/* determine which frequency that */
/* occurs the most */
/* > for lowest value when ties */
/* >= for higest value when ties */
max_freq_i = 0; /* don't return a zero value */
max_value_i = 0; /* unless all zeros */
for (k=0; k<=i; k++ )
{
if ( ((mfq_count[k]) > max_freq_i) && !(mfq_value[k] == 0) )
{
max_value_i = k;
max_freq_i = mfq_count[k];
}
}
cell[ncols] = mfq_value[max_value_i];
}
}
x_mode (argc, argv, cell, ncols)
double *argv[];
register double *cell;
register int ncols;
{
register int i, k, m, cols, max_freq_i ,max_value_i;
static int nargs = 0;
static double *arr = NULL;
static double *mfq_value = NULL;
static double *mfq_count = NULL;
if ( argc > nargs )
{
arr = (double *) G_realloc ( arr, argc*sizeof(double));
mfq_count = (double *) G_realloc ( mfq_count, argc*sizeof(double)+1);
mfq_value = (double *) G_realloc ( mfq_value, argc*sizeof(double)+1);
nargs = argc;
}
while (ncols-- > 0)
{
for ( i=0; i<argc; i++)
arr[i] = argv[i][ncols];
/* sort all values */
qsort( arr, argc, sizeof(double), dcmp);
/* sort by frequency */
for ( m=0; m<argc; m++)
mfq_count[m] = 0;
mfq_value[0] = arr[0];
mfq_count[0] = 1;
i = 0;
for ( m=0; m<argc; m++)
{
if ( arr[m] == arr[m -1])
mfq_count[i]++;
else
{
i++;
mfq_value[i] = arr[m];
mfq_count[i]++;
}
}
/* determine which frequency that */
/* occurs the most */
/* > for lowest value when ties */
/* >= for higest value when ties */
max_freq_i = 0; /* don't return a zero value */
max_value_i = 0; /* unless all zeros */
for (k=0; k<=i; k++ )
{
if ( ((mfq_count[k]) > max_freq_i) && !(mfq_value[k] == 0) )
{
max_value_i = k;
max_freq_i = mfq_count[k];
}
}
cell[ncols] = mfq_value[max_value_i];
}
}
n_mode (n,name) char *name;
{
if (n>1) return 1;
fprintf (stderr, "%s - ", name);
if (n==0)
fprintf (stderr, "no arguments ");
else
fprintf (stderr, "only one argument ");
fprintf (stderr, "specified. usage: %s(x,y[,z...])\n", name);
return 0;
}
/* end of xmode.c */
PGM = r.mapcalc
HOME=$(BIN_MAIN_CMD)
LIST = main.o\
compare.o\
convert.o\
evaluate.o\
execute.o\
expression.o\
function.o\
fpe.o\
input.o\
logical.o\
math.o\
maps.o\
min_max.o\
polish.o\
pool.o\
round.o\
support.o\
variables.o\
xexp.o\
x2float.o\
x2int.o\
xabs.o\
xcos.o\
xeval.o\
xif.o\
xlog.o\
xmax.o\
xmedian.o\
xmin.o\
xround.o\
xsqrt.o\
xsin.o\
xtan.o\
xatan.o\
xmode.o
LIBES = $(GISLIB) $(BTREELIB) $(ROWIOLIB)
$(HOME)/$(PGM): $(LIST) $(LIBES)
$(CC) $(LDFLAGS) $(LIST) $(LIBES) $(MATHLIB) -o $@
compare.o: glob.h
convert.o: glob.h
evaluate.o: glob.h
evaluate.o: function.h
execute.o: glob.h
expression.o: glob.h
fpe.o: glob.h
function.o: function.h
logical.o: glob.h
main.o: glob.h
main.o: function.h
math.o: glob.h
maps.o: glob.h
polish.o: glob.h
polish.o: function.h
round.o: glob.h
variables.o: glob.h
x2int.o: glob.h
xlog.o: glob.h
xsqrt.o: glob.h
xtan.o: glob.h
xatan.o: glob.h
xmode.o: function.h
$(GISLIB): #
$(BTREELIB): #
/**************************************************************************
the following function table describes the currently implemented functions
for the map calculator
the function table is defined as follows
name name of the function as specified by the user
icall routine to call if the input args are CELL buffers
if set to NOFUNC, all input args are converted to double
and the xcall routine is called.
xcall routine to call if the input args are doubles.
if set to NOFUNC, the result is a zeroed buffer.
nargs routine to call to check number of args
should return 1 if ok, 0 otherwise and print error msg.
double_result for xcall only. if true the result from the xcall routine
is double. else it is assumed to be CELL (ie, the
xcall routine produces CELL from doubles)
note: when writing a new function, make every effort to do something
reasonable with floating point exceptions. there is a mechanism
in place which makes this fairly simple to do. look at the
implementation of the log() function (xlog.c) to see how friendly
this function is.
***************************************************************************/
#define NOFUNC ( (int(*)())0 )
struct function_list
{
char *name;
int (*icall)(); /* cell arguments, cell result */
int (*xcall)(); /* double arguments */
int (*nargs)(); /* argument checker */
int double_result; /* xcall() produces double result or not */
};
#ifndef MAIN
extern struct function_list function_list[];
#else
/****************** FUNCTON DECLARATIONS *********************************/
int x_sqrt(), n_sqrt();
int x_exp(), n_exp();
int x_log(), n_log();
int x_cos(), n_cos();
int x_sin(), n_sin();
int x_tan(), n_tan();
int x_atan(), n_atan();
int i_abs(), x_abs(), n_abs();
int i_if(), x_if(), n_if();
int i_max(), x_max(), n_max();
int i_median(),x_median(), n_median();
int i_min(), x_min(), n_min();
int i_round(), x_round(), n_round();
int x_2float(), n_2float();
int i_2int(), x_2int(), n_2int();
int i_eval(), x_eval(), n_eval();
int i_mode(), x_mode(), n_mode();
struct function_list function_list[]=
{
"abs", i_abs, x_abs, n_abs, 1,
"eval", i_eval, x_eval, n_eval, 1,
"float", NOFUNC, x_2float, n_2float, 1,
"if", i_if, x_if, n_if, 1,
"int", i_2int, x_2int, n_2int, 0,
"exp", NOFUNC, x_exp, n_exp, 1,
"log", NOFUNC, x_log, n_log, 1,
"max", i_max, x_max, n_max, 1,
"min", i_min, x_min, n_min, 1,
"median",i_median,x_median, n_median, 1,
"round", i_round, x_round, n_round, 0,
"sqrt", NOFUNC, x_sqrt, n_sqrt, 1,
"sin", NOFUNC, x_sin, n_sin, 1,
"cos", NOFUNC, x_cos, n_cos, 1,
"tan", NOFUNC, x_tan, n_tan, 1,
"atan", NOFUNC, x_atan, n_atan, 1,
"mode", i_mode, x_mode, n_mode, 1,
/* add new functions above this line */
"", NOFUNC, NOFUNC, NOFUNC, 0
};
#endif
More information about the grass-dev
mailing list