[GRASS-SVN] r35963 - grass-addons/raster/r.univar2.zonal
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 19 18:20:08 EST 2009
Author: mmetz
Date: 2009-02-19 18:20:07 -0500 (Thu, 19 Feb 2009)
New Revision: 35963
Added:
grass-addons/raster/r.univar2.zonal/Makefile
grass-addons/raster/r.univar2.zonal/globals.h
grass-addons/raster/r.univar2.zonal/r.univar.html
grass-addons/raster/r.univar2.zonal/r.univar_main.c
grass-addons/raster/r.univar2.zonal/r3.univar.html
grass-addons/raster/r.univar2.zonal/sort.c
grass-addons/raster/r.univar2.zonal/stats.c
Log:
zonal statistics: here they are
Added: grass-addons/raster/r.univar2.zonal/Makefile
===================================================================
--- grass-addons/raster/r.univar2.zonal/Makefile (rev 0)
+++ grass-addons/raster/r.univar2.zonal/Makefile 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,23 @@
+
+MODULE_TOPDIR = ../..
+
+LIBES = $(G3DLIB) $(GISLIB)
+DEPENDENCIES = $(G3DDEP) $(GISDEP)
+
+#needed for htmlmulti
+PROGRAMS = r.univar.zonal
+
+include $(MODULE_TOPDIR)/include/Make/Multi.make
+
+#R3UNIVAR = $(BIN)/r3.univar.zonal$(EXE)
+RUNIVAR = $(BIN)/r.univar.zonal$(EXE)
+
+default: $(RUNIVAR)
+ $(MAKE) htmlmulti
+
+$(RUNIVAR): $(OBJDIR)/r.univar_main.o $(OBJDIR)/sort.o $(OBJDIR)/stats.o
+ $(CC) $(LDFLAGS) -o $@ $^ $(FMODE_OBJ) $(LIBES) $(XDRLIB) $(MATHLIB)
+
+#$(R3UNIVAR): $(OBJDIR)/r3.univar_main.o $(OBJDIR)/sort.o $(OBJDIR)/stats.o
+# $(CC) $(LDFLAGS) -o $@ $^ $(FMODE_OBJ) $(LIBES) $(XDRLIB) $(MATHLIB)
+
Added: grass-addons/raster/r.univar2.zonal/globals.h
===================================================================
--- grass-addons/raster/r.univar2.zonal/globals.h (rev 0)
+++ grass-addons/raster/r.univar2.zonal/globals.h 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,80 @@
+/*
+ * Calculates univariate statistics from the non-null cells
+ *
+ * Copyright (C) 2004-2007 by the GRASS Development Team
+ * Author(s): Soeren Gebbert
+ * Based on r.univar from Hamish Bowman, University of Otago, New Zealand
+ * and Martin Landa
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ */
+
+#ifndef _GLOBALS_H_
+#define _GLOBALS_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/G3d.h>
+#include <grass/glocale.h>
+
+/*- Parameters and global variables -----------------------------------------*/
+typedef struct
+{
+ double sum;
+ double sumsq;
+ double min;
+ double max;
+ unsigned int n_perc;
+ int *perc;
+ double sum_abs;
+ int n;
+ int size;
+ DCELL *dcell_array;
+ FCELL *fcell_array;
+ CELL *cell_array;
+ int map_type;
+} univar_stat;
+
+typedef struct
+{
+ CELL min, max, n_zones;
+ struct Categories cats;
+} zone_type;
+
+typedef struct
+{
+ int index;
+ void *nextp;
+} next_point;
+
+/* command line options are the same for raster and raster3d maps */
+typedef struct
+{
+ struct Option *inputfile, *zonefile, *percentile, *output_file;
+ struct Flag *shell_style, *extended;
+} param_type;
+
+#ifdef MAIN
+param_type param;
+zone_type zone_info;
+next_point *np;
+#else
+extern param_type param;
+extern zone_type zone_info;
+extern next_point *np;
+#endif
+
+/* fn prototypes */
+void heapsort_double(double *data, int n);
+void heapsort_float(float *data, int n);
+void heapsort_int(int *data, int n);
+int print_stats(univar_stat * stats);
+univar_stat *create_univar_stat_struct(int map_type, int size, int n_perc);
+void free_univar_stat_struct(univar_stat * stats);
+
+#endif
Added: grass-addons/raster/r.univar2.zonal/r.univar.html
===================================================================
--- grass-addons/raster/r.univar2.zonal/r.univar.html (rev 0)
+++ grass-addons/raster/r.univar2.zonal/r.univar.html 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,54 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.univar</em> calculates the univariate statistics of a raster map.
+This includes the number of cells counted, minimum and maximum cell values,
+range, arithmetic mean, population variance, standard deviation, and
+coefficient of variation.
+If the <b>-e</b> extended statistics flag is given the 1st quartile, median,
+3rd quartile, and given <b>percentile</b> are calculated.
+If the <b>-g</b> flag is given the results are presented in a format suitable
+for use in a shell script.
+
+<h2>NOTES</h2>
+
+As with most GRASS raster modules, <em>r.univar</em> operates on the cell
+array defined by the current region settings, not the original extent and
+resolution of the input map. See <em><a href="g.region.html">g.region</a></em>.
+<p>
+This module can use large amounts of system memory when the <b>-e</b>
+extended statistics flag is used with a very large region setting. If the
+region is too large the module should exit gracefully with a memory allocation
+error. Basic statistics can be calculated using any size input region.
+
+
+<h2>TODO</h2>
+
+<i>mode, skewness, kurtosis</i>
+
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a><br>
+<a href="r3.univar.html">r3.univar</a><br>
+<a href="r.univar.sh.html">r.univar.sh</a><br>
+<a href="r.average.html">r.average</a><br>
+<a href="r.median.html">r.median</a><br>
+<a href="r.mode.html">r.mode</a><br>
+<a href="r.sum.html">r.sum</a><br>
+<a href="r.series.html">r.series</a><br>
+<a href="r.stats.html">r.stats</a><br>
+<a href="r.statistics.html">r.statistics</a><br>
+<a href="v.univar.html">v.univar</a><br>
+<a href="v.univar.sh.html">v.univar.sh</a><br>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Hamish Bowman, Otago University, New Zealand<br>
+Extended statistics by Martin Landa
+
+<p>
+<i>Last changed: $Date: 2008-05-16 21:11:45 +0200 (Fri, 16 May 2008) $</i>
Added: grass-addons/raster/r.univar2.zonal/r.univar_main.c
===================================================================
--- grass-addons/raster/r.univar2.zonal/r.univar_main.c (rev 0)
+++ grass-addons/raster/r.univar2.zonal/r.univar_main.c 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,343 @@
+/*
+ * r.univar
+ *
+ * Calculates univariate statistics from the non-null cells of a GRASS raster map
+ *
+ * Copyright (C) 2004-2006 by the GRASS Development Team
+ * Author(s): Hamish Bowman, University of Otago, New Zealand
+ * Extended stats: Martin Landa
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ * This program is a replacement for the r.univar shell script
+ */
+
+#include <assert.h>
+#include <string.h>
+#define MAIN
+#include "globals.h"
+
+/* local proto */
+void set_params();
+
+/* ************************************************************************* */
+/* Set up the arguments we are expecting ********************************** */
+/* ************************************************************************* */
+void set_params()
+{
+ param.inputfile = G_define_standard_option(G_OPT_R_INPUT);
+
+ param.zonefile = G_define_standard_option(G_OPT_R_MAP);
+ param.zonefile->key = "zones";
+ param.zonefile->required = YES;
+ param.zonefile->description =
+ _("Raster map used for zoning, must be of type CELL");
+
+ param.output_file = G_define_standard_option(G_OPT_F_OUTPUT);
+ param.output_file->required = NO;
+ param.output_file->description =
+ _("Name for output file (if omitted or \"-\" output to stdout)");
+
+ param.percentile = G_define_option();
+ param.percentile->key = "percentile";
+ param.percentile->type = TYPE_INTEGER;
+ param.percentile->required = NO;
+ param.percentile->multiple = YES;
+ param.percentile->options = "0-100";
+ param.percentile->answer = "90";
+ param.percentile->description =
+ _("Percentile to calculate (requires extended statistics flag)");
+
+ param.shell_style = G_define_flag();
+ param.shell_style->key = 'g';
+ param.shell_style->description =
+ _("Print the stats in shell script style");
+
+ param.extended = G_define_flag();
+ param.extended->key = 'e';
+ param.extended->description = _("Calculate extended statistics");
+
+ return;
+}
+
+static int open_raster(const char *infile);
+static univar_stat *univar_stat_with_percentiles(int map_type, int size);
+static void process_raster(univar_stat * stats, int fd, int fdz,
+ const struct Cell_head *region);
+
+/* *************************************************************** */
+/* **** the main functions for r.univar ************************** */
+/* *************************************************************** */
+int main(int argc, char *argv[])
+{
+ unsigned int rows, cols; /* totals */
+ /* int rasters; */
+
+ struct Cell_head region;
+ struct GModule *module;
+ univar_stat *stats;
+ char *p, *z;
+ int fd, fdz, cell_type, min, max;
+ struct Range zone_range;
+ char *mapset, *name;
+
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster, statistics");
+ module->description =
+ _("Calculates univariate statistics from the non-null cells of a raster map.");
+
+ /* Define the different options */
+ set_params();
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ name = param.output_file->answer;
+ if (name != NULL && strcmp(name, "-") != 0) {
+ if (NULL == freopen(name, "w", stdout)) {
+ G_fatal_error(_("Unable to open file <%s> for writing"), name);
+ }
+ }
+
+ G_get_window(®ion);
+ rows = region.rows;
+ cols = region.cols;
+
+ zone_info.min = 0.0 / 0.0; /*set to nan as default */
+ zone_info.max = 0.0 / 0.0; /*set to nan as default */
+ zone_info.n_zones = 0;
+
+ /* open zoning raster */
+ G_message("find zoning raster");
+ z = param.zonefile->answer;
+ mapset = G_find_cell2(z, "");
+
+ G_message("open zoning raster");
+ fdz = open_raster(z);
+
+ cell_type = G_get_raster_map_type(fdz);
+ if (cell_type != CELL_TYPE)
+ G_fatal_error("Zoning raster must be of type CELL");
+
+ if (G_read_range(z, mapset, &zone_range) == -1)
+ G_fatal_error("Can not read range for zoning raster");
+ if (G_get_range_min_max(&zone_range, &min, &max))
+ G_fatal_error("Can not read range for zoning raster");
+ G_message("get cats");
+ if (G_read_raster_cats(z, mapset, &(zone_info.cats)))
+ G_warning("no category support for zoning raster");
+
+ zone_info.min = min;
+ zone_info.max = max;
+ zone_info.n_zones = max - min + 1;
+
+ /* process input raster */
+ G_message("process input raster");
+ size_t cells = cols * rows;
+ int map_type = param.extended->answer ? -2 : -1;
+
+ stats = ((map_type == -1)
+ ? create_univar_stat_struct(-1, cells, 0)
+ : 0);
+
+ p = param.inputfile->answer;
+ fd = open_raster(p);
+
+ if (map_type != -1) {
+ /* NB: map_type must match when doing extended stats */
+ int this_type = G_get_raster_map_type(fd);
+
+ assert(this_type > -1);
+ if (map_type < -1) {
+ assert(stats == 0);
+ map_type = this_type;
+ stats = univar_stat_with_percentiles(map_type, cells);
+ }
+ else if (this_type != map_type) {
+ G_fatal_error(_("Raster <%s> type mismatch"), p);
+ }
+ }
+
+ G_message("process_raster");
+ process_raster(stats, fd, fdz, ®ion);
+
+ if (!(param.shell_style->answer))
+ G_percent(rows, rows, 2); /* finish it off */
+
+ /* create the output */
+ G_message("print_stats");
+ print_stats(stats);
+
+ /* release memory */
+ free_univar_stat_struct(stats);
+
+ exit(EXIT_SUCCESS);
+}
+
+static int open_raster(const char *infile)
+{
+ char *mapset;
+ int fd;
+
+ mapset = G_find_cell2(infile, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Raster map <%s> not found"), infile);
+ }
+
+ fd = G_open_cell_old(infile, mapset);
+ if (fd < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), infile);
+
+ /* . */
+ return fd;
+}
+
+static univar_stat *univar_stat_with_percentiles(int map_type, int size)
+{
+ univar_stat *stats;
+ int i, j;
+
+ i = 0;
+ while (param.percentile->answers[i])
+ i++;
+ stats = create_univar_stat_struct(map_type, size, i);
+ for (i = 0; i < zone_info.n_zones; i++) {
+ for (j = 0; j < stats[i].n_perc; j++) {
+ sscanf(param.percentile->answers[j], "%i", &stats[i].perc[j]);
+ }
+ }
+
+ /* . */
+ return stats;
+}
+
+static void
+process_raster(univar_stat * stats, int fd, int fdz, const struct Cell_head *region)
+{
+ /* use G_window_rows(), G_window_cols() here? */
+ const int rows = region->rows;
+ const int cols = region->cols;
+ const int n_zones = zone_info.n_zones;
+ int *first, i;
+
+ const RASTER_MAP_TYPE map_type = G_get_raster_map_type(fd);
+ const size_t value_sz = G_raster_size(map_type);
+ unsigned int row;
+ void *raster_row;
+ CELL *zoneraster_row;
+
+ G_message("initialize processing");
+ np = (next_point *) G_calloc(n_zones, sizeof(next_point));
+
+ switch (map_type) {
+ case CELL_TYPE:
+ for (i = 0; i < n_zones; i++) {
+ np[i].nextp
+ = ((!param.extended->answer) ? 0 : (void *)stats[i].cell_array);
+ }
+ break;
+
+ case FCELL_TYPE:
+ for (i = 0; i < n_zones; i++) {
+ np[i].nextp
+ = ((!param.extended->answer) ? 0 : (void *)stats[i].fcell_array);
+ }
+ break;
+
+ case DCELL_TYPE:
+ for (i = 0; i < n_zones; i++) {
+ np[i].nextp
+ = ((!param.extended->answer) ? 0 : (void *)stats[i].dcell_array);
+ }
+ break;
+
+ default:
+ G_fatal_error("input raster has unknown type");
+ break;
+ }
+
+ G_message("start processing");
+
+ first = (int *) G_malloc((n_zones - 1) * sizeof(int));
+ {
+ int i;
+ for (i = 0; i < n_zones; i++)
+ first[i] = (stats[i].n < 1);
+ }
+
+
+ raster_row = G_calloc(cols, value_sz);
+ zoneraster_row = G_allocate_c_raster_buf();
+
+ for (row = 0; row < rows; row++) {
+ void *ptr;
+ CELL *zptr;
+ unsigned int col;
+ int zone;
+
+ if (G_get_raster_row(fd, raster_row, row, map_type) < 0)
+ G_fatal_error(_("Reading row %d"), row);
+ if (G_get_c_raster_row(fdz, zoneraster_row, row) < 0)
+ G_fatal_error(_("Reading row %d"), row);
+
+ ptr = raster_row;
+ zptr = zoneraster_row;
+
+ for (col = 0; col < cols; col++) {
+
+ if (G_is_null_value(ptr, map_type)) {
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ continue;
+ }
+ if (G_is_c_null_value(zptr)) {
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ continue;
+ }
+
+ zone = *zptr + zone_info.min;
+
+ if (np[zone].nextp) {
+ /* put the value into stats->XXXcell_array */
+ memcpy(np[zone].nextp, ptr, value_sz);
+ np[zone].nextp = G_incr_void_ptr(np[zone].nextp, value_sz);
+ }
+
+ {
+ double val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr)
+ : (map_type == FCELL_TYPE) ? *((FCELL *) ptr)
+ : *((CELL *) ptr));
+
+ stats[zone].sum += val;
+ stats[zone].sumsq += val * val;
+ stats[zone].sum_abs += fabs(val);
+
+ if (first[zone]) {
+ stats[zone].max = val;
+ stats[zone].min = val;
+ first[zone] = FALSE;
+ }
+ else {
+ if (val > stats[zone].max)
+ stats[zone].max = val;
+ if (val < stats[zone].min)
+ stats[zone].min = val;
+ }
+ }
+
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ stats[zone].n++;
+ }
+ if (!(param.shell_style->answer))
+ G_percent(row, rows, 2);
+ }
+ G_message("finished processing");
+
+}
Added: grass-addons/raster/r.univar2.zonal/r3.univar.html
===================================================================
--- grass-addons/raster/r.univar2.zonal/r3.univar.html (rev 0)
+++ grass-addons/raster/r.univar2.zonal/r3.univar.html 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,57 @@
+<h2>DESCRIPTION</h2>
+
+<em>r3.univar</em> calculates the univariate statistics for raster3d maps.
+This includes the number of cells counted, minimum and maximum cell values,
+range, arithmetic mean, population variance, standard deviation, and
+coefficient of variation.
+If the <b>-e</b> extended statistics flag is given the 1st quartile, median,
+3rd quartile, and given <b>percentile</b> are calculated.
+If the <b>-g</b> flag is given the results are presented in a format suitable
+for use in a shell script.
+
+<h2>NOTES</h2>
+
+As with most GRASS raster3d modules, <em>r3.univar</em> operates on the cell
+array defined by the current 3d region settings, not the original extent and
+resolution of the input map. See <em><a href="g.region.html">g.region</a></em>.
+<p>
+This module can use large amounts of system memory when the <b>-e</b>
+extended statistics flag is used with a very large region setting. If the
+region is too large the module should exit gracefully with a memory allocation
+error. Basic statistics can be calculated using any size input region.
+
+
+<h2>TODO</h2>
+
+<i>mode, skewness, kurtosis</i>
+
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a><br>
+<a href="r.univar.html">r.univar</a><br>
+<a href="r.univar.sh.html">r.univar.sh</a><br>
+<a href="r.average.html">r.average</a><br>
+<a href="r.median.html">r.median</a><br>
+<a href="r.mode.html">r.mode</a><br>
+<a href="r.sum.html">r.sum</a><br>
+<a href="r.series.html">r.series</a><br>
+<a href="r.stats.html">r.stats</a><br>
+<a href="r.statistics.html">r.statistics</a><br>
+<a href="v.univar.html">v.univar</a><br>
+<a href="v.univar.sh.html">v.univar.sh</a><br>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Soeren Gebbert<br>
+Code is based on r.univar from<br>
+Hamish Bowman, Otago University, New Zealand<br>
+and Martin Landa
+
+
+<p>
+<i>Last changed: $Date: 2008-05-16 21:11:45 +0200 (Fri, 16 May 2008) $</i>
Added: grass-addons/raster/r.univar2.zonal/sort.c
===================================================================
--- grass-addons/raster/r.univar2.zonal/sort.c (rev 0)
+++ grass-addons/raster/r.univar2.zonal/sort.c 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,152 @@
+/*
+ * Copyright (C) 2004-2007 by the GRASS Development Team
+ * Author(s): 1998, 1999 Sebastian Cyris
+ * 2007 modified by Soeren Gebbert
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ */
+
+#include "globals.h"
+static void downheap_int(int *array, int n, int k);
+static void downheap_float(float *array, int n, int k);
+static void downheap_double(double *array, int n, int k);
+
+/* *************************************************************** */
+/* *************************************************************** */
+/* *************************************************************** */
+void downheap_int(int *array, int n, int k)
+{
+ int j, v;
+
+ v = array[k];
+ while (k <= n / 2) {
+ j = k + k;
+
+ if (j < n && array[j] < array[j + 1])
+ j++;
+ if (v >= array[j])
+ break;
+
+ array[k] = array[j];
+ k = j;
+ }
+
+ array[k] = v;
+}
+
+/* *************************************************************** */
+/* *************************************************************** */
+/* *************************************************************** */
+void downheap_float(float *array, int n, int k)
+{
+ int j;
+ float v;
+
+ v = array[k];
+ while (k <= n / 2) {
+ j = k + k;
+
+ if (j < n && array[j] < array[j + 1])
+ j++;
+ if (v >= array[j])
+ break;
+
+ array[k] = array[j];
+ k = j;
+ }
+
+ array[k] = v;
+}
+
+/* *************************************************************** */
+/* *************************************************************** */
+/* *************************************************************** */
+void downheap_double(double *array, int n, int k)
+{
+ int j;
+ double v;
+
+ v = array[k];
+ while (k <= n / 2) {
+ j = k + k;
+
+ if (j < n && array[j] < array[j + 1])
+ j++;
+ if (v >= array[j])
+ break;
+
+ array[k] = array[j];
+ k = j;
+ }
+
+ array[k] = v;
+}
+
+/* *************************************************************** */
+/* ****** heapsort for int arrays of size n ********************** */
+/* *************************************************************** */
+void heapsort_int(int *array, int n)
+{
+ int k, t;
+
+ --n;
+
+ for (k = n / 2; k >= 0; k--)
+ downheap_int(array, n, k);
+
+ while (n > 0) {
+ t = array[0];
+ array[0] = array[n];
+ array[n] = t;
+ downheap_int(array, --n, 0);
+ }
+ return;
+}
+
+
+/* *************************************************************** */
+/* ****** heapsort for float arrays of size n ******************** */
+/* *************************************************************** */
+void heapsort_float(float *array, int n)
+{
+ int k;
+ float t;
+
+ --n;
+
+ for (k = n / 2; k >= 0; k--)
+ downheap_float(array, n, k);
+
+ while (n > 0) {
+ t = array[0];
+ array[0] = array[n];
+ array[n] = t;
+ downheap_float(array, --n, 0);
+ }
+ return;
+}
+
+/* *************************************************************** */
+/* ****** heapsort for double arrays of size n ******************* */
+/* *************************************************************** */
+void heapsort_double(double *array, int n)
+{
+ int k;
+ double t;
+
+ --n;
+
+ for (k = n / 2; k >= 0; k--)
+ downheap_double(array, n, k);
+
+ while (n > 0) {
+ t = array[0];
+ array[0] = array[n];
+ array[n] = t;
+ downheap_double(array, --n, 0);
+ }
+ return;
+}
Added: grass-addons/raster/r.univar2.zonal/stats.c
===================================================================
--- grass-addons/raster/r.univar2.zonal/stats.c (rev 0)
+++ grass-addons/raster/r.univar2.zonal/stats.c 2009-02-19 23:20:07 UTC (rev 35963)
@@ -0,0 +1,261 @@
+/*
+ * Calculates univariate statistics from the non-null cells
+ *
+ * Copyright (C) 2004-2007 by the GRASS Development Team
+ * Author(s): Hamish Bowman, University of Otago, New Zealand
+ * Martin Landa and Soeren Gebbert
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ */
+
+#include "globals.h"
+
+/* *************************************************************** */
+/* **** univar_stat constructor ********************************** */
+/* *************************************************************** */
+univar_stat *create_univar_stat_struct(int map_type, int size, int n_perc)
+{
+ univar_stat *stats;
+ int i;
+ int n_zones = zone_info.n_zones;
+
+ stats = (univar_stat *) G_calloc(n_zones, sizeof(univar_stat));
+
+ for (i = 0; i < n_zones; i++) {
+ stats[i].sum = 0.0;
+ stats[i].sumsq = 0.0;
+ stats[i].min = 0.0 / 0.0; /*set to nan as default */
+ stats[i].max = 0.0 / 0.0; /*set to nan as default */
+ stats[i].n_perc = n_perc;
+ if (n_perc > 0)
+ stats[i].perc = (int *)G_malloc(n_perc * sizeof(int));
+ else
+ stats[i].perc = NULL;
+ stats[i].sum_abs = 0.0;
+ stats[i].n = 0;
+ stats[i].size = size;
+ stats[i].dcell_array = NULL;
+ stats[i].fcell_array = NULL;
+ stats[i].cell_array = NULL;
+ stats[i].map_type = map_type;
+
+ /* alloacte memory for extended computation */
+ if (param.extended->answer) {
+ if (map_type == DCELL_TYPE)
+ stats[i].dcell_array =
+ (DCELL *) G_calloc(stats[i].size, sizeof(DCELL));
+ if (map_type == FCELL_TYPE)
+ stats[i].fcell_array =
+ (FCELL *) G_calloc(stats[i].size, sizeof(FCELL));
+ if (map_type == CELL_TYPE)
+ stats[i].cell_array = (CELL *) G_calloc(stats[i].size, sizeof(CELL));
+ }
+ }
+
+ return stats;
+}
+
+
+/* *************************************************************** */
+/* **** univar_stat destructor *********************************** */
+/* *************************************************************** */
+void free_univar_stat_struct(univar_stat * stats)
+{
+ int i;
+
+ for (i = 0; i < zone_info.n_zones; i++){
+ if (stats[i].perc)
+ G_free(stats[i].perc);
+ if (stats[i].dcell_array)
+ G_free(stats[i].dcell_array);
+ if (stats[i].fcell_array)
+ G_free(stats[i].fcell_array);
+ if (stats[i].cell_array)
+ G_free(stats[i].cell_array);
+ }
+
+ G_free(stats);
+
+ return;
+}
+
+
+/* *************************************************************** */
+/* **** compute and print univar statistics to stdout ************ */
+/* *************************************************************** */
+int print_stats(univar_stat * stats)
+{
+ int z;
+
+ for (z = 0; z < zone_info.n_zones; z++) {
+ char sum_str[100];
+ double mean, variance, stdev, var_coef;
+
+ /*for extendet stats */
+ double quartile_25 = 0.0, quartile_75 = 0.0, *quartile_perc;
+ double median = 0.0;
+ unsigned int i;
+ int qpos_25, qpos_75, *qpos_perc;
+
+ /* stats collected for this zone? */
+ if (stats[z].n == 0)
+ continue;
+
+
+ /* all these calculations get promoted to doubles, so any DIV0 becomes nan */
+ mean = stats[z].sum / stats[z].n;
+ variance = (stats[z].sumsq - stats[z].sum * stats[z].sum / stats[z].n) / stats[z].n;
+ if (variance < GRASS_EPSILON)
+ variance = 0.0;
+ stdev = sqrt(variance);
+ var_coef = (stdev / mean) * 100.; /* perhaps stdev/fabs(mean) ? */
+
+ sprintf(sum_str, "%.10f", stats[z].sum);
+ G_trim_decimal(sum_str);
+
+ fprintf(stdout, "\nzone %d %s\n\n", z - zone_info.min, G_get_c_raster_cat(&z + zone_info.min, &(zone_info.cats)));
+
+ if (!param.shell_style->answer) {
+ fprintf(stdout, "total null and non-null cells: %d\n", stats[z].size);
+ fprintf(stdout, "total null cells: %d\n\n", stats[z].size - stats[z].n);
+ fprintf(stdout, "Of the non-null cells:\n----------------------\n");
+ }
+
+ if (param.shell_style->answer) {
+ fprintf(stdout, "n=%d\n", stats[z].n);
+ fprintf(stdout, "null_cells=%d\n", stats[z].size - stats[z].n);
+ fprintf(stdout, "min=%.15g\n", stats[z].min);
+ fprintf(stdout, "max=%.15g\n", stats[z].max);
+ fprintf(stdout, "range=%.15g\n", stats[z].max - stats[z].min);
+ fprintf(stdout, "mean=%.15g\n", mean);
+ fprintf(stdout, "mean_of_abs=%.15g\n", stats[z].sum_abs / stats[z].n);
+ fprintf(stdout, "stddev=%.15g\n", stdev);
+ fprintf(stdout, "variance=%.15g\n", variance);
+ fprintf(stdout, "coeff_var=%.15g\n", var_coef);
+ fprintf(stdout, "sum=%s\n", sum_str);
+ }
+ else {
+ fprintf(stdout, "n: %d\n", stats[z].n);
+ fprintf(stdout, "minimum: %g\n", stats[z].min);
+ fprintf(stdout, "maximum: %g\n", stats[z].max);
+ fprintf(stdout, "range: %g\n", stats[z].max - stats[z].min);
+ fprintf(stdout, "mean: %g\n", mean);
+ fprintf(stdout, "mean of absolute values: %g\n",
+ stats[z].sum_abs / stats[z].n);
+ fprintf(stdout, "standard deviation: %g\n", stdev);
+ fprintf(stdout, "variance: %g\n", variance);
+ fprintf(stdout, "variation coefficient: %g %%\n", var_coef);
+ fprintf(stdout, "sum: %s\n", sum_str);
+ }
+
+
+ /* TODO: mode, skewness, kurtosis */
+ if (param.extended->answer) {
+ qpos_perc = (int *)G_calloc(stats[z].n_perc, sizeof(int));
+ quartile_perc = (double *)G_calloc(stats[z].n_perc, sizeof(double));
+ for (i = 0; i < stats[z].n_perc; i++) {
+ qpos_perc[i] = (int)(stats[z].n * 1e-2 * stats[z].perc[i] - 0.5);
+ }
+ qpos_25 = (int)(stats[z].n * 0.25 - 0.5);
+ qpos_75 = (int)(stats[z].n * 0.75 - 0.5);
+
+ switch (stats[z].map_type) {
+ case CELL_TYPE:
+ heapsort_int(stats[z].cell_array, stats[z].n);
+
+ quartile_25 = (double)stats[z].cell_array[qpos_25];
+ if (stats[z].n % 2) /* odd */
+ median = (double)stats[z].cell_array[(int)(stats[z].n / 2)];
+ else /* even */
+ median =
+ (double)(stats[z].cell_array[stats[z].n / 2 - 1] +
+ stats[z].cell_array[stats[z].n / 2]) / 2.0;
+ quartile_75 = (double)stats[z].cell_array[qpos_75];
+ for (i = 0; i < stats[z].n_perc; i++) {
+ quartile_perc[i] = (double)stats[z].cell_array[qpos_perc[i]];
+ }
+ break;
+
+ case FCELL_TYPE:
+ heapsort_float(stats[z].fcell_array, stats[z].n);
+
+ quartile_25 = (double)stats[z].fcell_array[qpos_25];
+ if (stats[z].n % 2) /* odd */
+ median = (double)stats[z].fcell_array[(int)(stats[z].n / 2)];
+ else /* even */
+ median =
+ (double)(stats[z].fcell_array[stats[z].n / 2 - 1] +
+ stats[z].fcell_array[stats[z].n / 2]) / 2.0;
+ quartile_75 = (double)stats[z].fcell_array[qpos_75];
+ for (i = 0; i < stats[z].n_perc; i++) {
+ quartile_perc[i] = (double)stats[z].fcell_array[qpos_perc[i]];
+ }
+ break;
+
+ case DCELL_TYPE:
+ heapsort_double(stats[z].dcell_array, stats[z].n);
+
+ quartile_25 = stats[z].dcell_array[qpos_25];
+ if (stats[z].n % 2) /* odd */
+ median = stats[z].dcell_array[(int)(stats[z].n / 2)];
+ else /* even */
+ median =
+ (stats[z].dcell_array[stats[z].n / 2 - 1] +
+ stats[z].dcell_array[stats[z].n / 2]) / 2.0;
+ quartile_75 = stats[z].dcell_array[qpos_75];
+ for (i = 0; i < stats[z].n_perc; i++) {
+ quartile_perc[i] = stats[z].dcell_array[qpos_perc[i]];
+ }
+ break;
+
+ default:
+ break;
+ }
+
+ if (param.shell_style->answer) {
+ fprintf(stdout, "first_quartile=%g\n", quartile_25);
+ fprintf(stdout, "median=%g\n", median);
+ fprintf(stdout, "third_quartile=%g\n", quartile_75);
+ for (i = 0; i < stats[z].n_perc; i++) {
+ fprintf(stdout, "percentile_%d=%g\n", stats[z].perc[i],
+ quartile_perc[i]);
+ }
+ }
+ else {
+ fprintf(stdout, "1st quartile: %g\n", quartile_25);
+ if (stats[z].n % 2)
+ fprintf(stdout, "median (odd number of cells): %g\n", median);
+ else
+ fprintf(stdout, "median (even number of cells): %g\n",
+ median);
+ fprintf(stdout, "3rd quartile: %g\n", quartile_75);
+
+
+ for (i = 0; i < stats[z].n_perc; i++) {
+ if (stats[z].perc[i] % 10 == 1 && stats[z].perc[i] != 11)
+ fprintf(stdout, "%dst percentile: %g\n", stats[z].perc[i],
+ quartile_perc[i]);
+ else if (stats[z].perc[i] % 10 == 2 && stats[z].perc[i] != 12)
+ fprintf(stdout, "%dnd percentile: %g\n", stats[z].perc[i],
+ quartile_perc[i]);
+ else if (stats[z].perc[i] % 10 == 3 && stats[z].perc[i] != 13)
+ fprintf(stdout, "%drd percentile: %g\n", stats[z].perc[i],
+ quartile_perc[i]);
+ else
+ fprintf(stdout, "%dth percentile: %g\n", stats[z].perc[i],
+ quartile_perc[i]);
+ }
+ }
+ G_free((void *)quartile_perc);
+ G_free((void *)qpos_perc);
+ }
+
+ if (!(param.shell_style->answer))
+ G_message("\n");
+ }
+
+ return 1;
+}
More information about the grass-commit
mailing list