[GRASS-SVN] r58439 - in grass-addons/grass7/raster: . r.univar2
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Dec 10 14:25:27 PST 2013
Author: zarch
Date: 2013-12-10 14:25:27 -0800 (Tue, 10 Dec 2013)
New Revision: 58439
Added:
grass-addons/grass7/raster/r.univar2/
grass-addons/grass7/raster/r.univar2/Makefile
grass-addons/grass7/raster/r.univar2/globals.h
grass-addons/grass7/raster/r.univar2/r.univar.html
grass-addons/grass7/raster/r.univar2/r.univar_main.c
grass-addons/grass7/raster/r.univar2/sort.c
grass-addons/grass7/raster/r.univar2/stats.c
Log:
Add a modified version of r.univar, at the moment works obly whit extended statistics and zones
Added: grass-addons/grass7/raster/r.univar2/Makefile
===================================================================
--- grass-addons/grass7/raster/r.univar2/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.univar2/Makefile 2013-12-10 22:25:27 UTC (rev 58439)
@@ -0,0 +1,15 @@
+
+MODULE_TOPDIR = ../..
+
+PGM = r.univar2
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+PROGRAMS = r.univar2
+
+r_univar_OBJS = r.univar_main.o sort.o stats.o
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.univar2/globals.h
===================================================================
--- grass-addons/grass7/raster/r.univar2/globals.h (rev 0)
+++ grass-addons/grass7/raster/r.univar2/globals.h 2013-12-10 22:25:27 UTC (rev 58439)
@@ -0,0 +1,115 @@
+/*
+ * Calculates univariate statistics from the non-null cells
+ *
+ * Copyright (C) 2004-2010 by the GRASS Development Team
+ * Author(s): Soeren Gebbert
+ * Based on r.univar from Hamish Bowman, University of Otago, New Zealand
+ * and Martin Landa
+ * zonal loop by Markus Metz
+ *
+ * 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 "../r.flow/mem.h"
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+/*- Parameters and global variables -----------------------------------------*/
+typedef struct
+{
+ int zone;
+ char *cat;
+ int null_cells;
+ double sum;
+ double sum2;
+ double sum3;
+ double sum4;
+ double sum_abs;
+ double min;
+ double max;
+
+ double range;
+ double mean;
+ double meandev;
+
+ double variance;
+ double stddev;
+ double var_coef;
+
+ double skewness;
+ double kurtosis;
+
+ double variance2;
+ double stddev2;
+ double var_coef2;
+
+ double skewness2;
+ double kurtosis2;
+
+ /* extend statistics */
+ double quartile_25;
+ double median;
+ double quartile_75;
+ double *perc;
+ double mode;
+ int occurrences;
+
+ /* need for processing */
+ int n;
+ int size;
+
+ RASTER_MAP_TYPE map_type;
+ DCELL *array;
+ void *nextp;
+ int n_alloc;
+} univar_stat;
+
+typedef struct
+{
+ CELL min, max, n_zones;
+ struct Categories cats;
+ char *sep;
+ int *len;
+ int n_alloc;
+} zone_type;
+
+/* command line options are the same for raster and raster3d maps */
+typedef struct
+{
+ struct Option *inputfile, *zonefile, *percentile, *tolerance, *output_file, *separator;
+ struct Flag *shell_style, *extended, *table;
+ int n_perc;
+ int *index_perc;
+ double *quant_perc;
+ double *perc;
+ double tol;
+} param_type;
+
+extern param_type param;
+extern zone_type zone_info;
+
+/* fn prototypes */
+void heapsort_double(double *data, int n);
+void heapsort_float(float *data, int n);
+void heapsort_int(int *data, int n);
+
+int compute_stats(univar_stat *, double, int);
+
+/* int print_stats(univar_stat *); */
+int print_stats_table(univar_stat *);
+
+univar_stat *create_univar_stat_struct();
+void free_univar_stat_struct(univar_stat * stats);
+
+
+#endif
Added: grass-addons/grass7/raster/r.univar2/r.univar.html
===================================================================
--- grass-addons/grass7/raster/r.univar2/r.univar.html (rev 0)
+++ grass-addons/grass7/raster/r.univar2/r.univar.html 2013-12-10 22:25:27 UTC (rev 58439)
@@ -0,0 +1,136 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.univar</em> calculates the univariate statistics of one or several raster
+map(s). This includes the number of cells counted, minimum and maximum cell
+values, range, arithmetic mean, population variance, standard deviation, and
+coefficient of variation. Statistics are calculated separately for every
+category/zone found in the <b>zones</b> input map if given.
+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.
+If the <b>-t</b> flag is given the results are presented in tabular format
+with the given field separator. The table can immediately be converted to a
+vector attribute table which can then be linked to a vector, e.g. the vector
+that was rasterized to create the <b>zones</b> input raster.
+<p>
+When multiple input maps are given to <em>r.univar</em>, the overall statistics
+are calculated. This is useful for a time series of the same variable, as well as
+for the case of a segmented/tiled dataset. Allowing multiple raster maps to be
+specified saves the user from using a temporary raster map for the result of
+<em>r.series</em> or <em>r.patch</em>.
+
+<h2>NOTES</h2>
+
+As with most GRASS raster modules, <em>r.univar</em> operates on the raster
+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.
+<p>
+Without a <b>zones</b> input raster, the <em>r.quantile</em> module will
+be significantly more efficient for calculating percentiles with large maps.
+
+<h2>EXAMPLE</h2>
+
+Calculate the raster statistics for zones within a vector map coverage
+and upload the results for mean, min and max back to the vector map:
+
+<div class="code"><pre>
+#### set the raster region to match the map
+g.region vect=fields res=10 -ap
+
+#### create rasterized version of vector map
+v.to.rast in=fields out=fields.10m use=cat type=area labelcolumn=label
+r.colors fields.10m color=random
+
+#### perform analysis
+r.univar -t map=elevation.10m zones=fields.10m | \
+ cut -f1,5,6,8 -d'|' > fields_stats.txt
+
+
+#### populate vector DB with stats
+
+# create working copy of vector map
+g.copy vect=fields,fields_stats
+
+# create new attribute columns to hold output
+v.db.addcol map=fields_stats \
+ columns='mean_elev DOUBLE PRECISION, min_elev DOUBLE PRECISION, max_elev DOUBLE PRECISION'
+
+# create SQL command file, and execute it
+sed -e '1d' fields_stats.txt | awk -F'|' \
+ '{print "UPDATE fields_stats SET min_elev = "$2", max_elev = "$3", \
+ mean_elev = "$4" WHERE cat = "$1";"}' \
+ > fields_stats_sqlcmd.txt
+
+db.execute input=fields_stats_sqlcmd.txt
+<!--
+
+###### alternate method with db.in.ogr: (needs work) ######
+
+#### convert text file table to a database table
+# not safe for commas in the label
+tr '|' ',' < fields_stats.txt > fields_stats.csv
+echo '"Integer","String","Real","Real","Real"' > fields_stats.csvt
+
+# import table
+db.in.ogr dsn=fields_stats.csv output=fields_data
+
+# view table
+db.select fields_data
+
+# remove temporary files
+rm fields_stats.csv fields_stats.csvt fields_stats.txt
+
+
+#### populate vector DB with stats
+
+# create working copy of vector map
+g.copy vect=fields,fields_stats
+
+# create new attribute columns to hold output
+v.db.addcol map=fields_stats \
+ columns='mean_elev DOUBLE PRECISION, min_elev DOUBLE PRECISION, max_elev DOUBLE PRECISION'
+
+# perform DB step (broken)
+## how to automatically collate by key column, ie copy between tables?
+## SELECT INTO? JOIN?
+echo "INSERT INTO fields_stats (mean_elev,min_elev,max_elev) SELECT mean,min,max FROM fields_data" | db.execute
+-->
+
+#### view completed table
+v.db.select fields_stats
+</pre></div>
+
+
+<h2>TODO</h2>
+
+<i>mode, skewness, kurtosis</i>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="g.region.html">g.region</a>,
+<a href="r3.univar.html">r3.univar</a>,
+<a href="r.mode.html">r.mode</a>,
+<a href="r.quantile.html">r.quantile</a>,
+<a href="r.series.html">r.series</a>,
+<a href="r.stats.html">r.stats</a>,
+<a href="r.statistics.html">r.statistics</a>,
+<a href="v.rast.stats.html">v.rast.stats</a>,
+<a href="v.univar.html">v.univar</a>
+</em>
+
+
+<h2>AUTHORS</h2>
+
+Hamish Bowman, Otago University, New Zealand<br>
+Extended statistics by Martin Landa<br>
+Multiple input map support by Ivan Shmakov<br>
+Zonal loop by Markus Metz
+
+<p><i>Last changed: $Date: 2013-02-07 09:56:17 +0000 (Thu, 07 Feb 2013) $</i>
Added: grass-addons/grass7/raster/r.univar2/r.univar_main.c
===================================================================
--- grass-addons/grass7/raster/r.univar2/r.univar_main.c (rev 0)
+++ grass-addons/grass7/raster/r.univar2/r.univar_main.c 2013-12-10 22:25:27 UTC (rev 58439)
@@ -0,0 +1,382 @@
+/*
+ * r.univar
+ *
+ * Calculates univariate statistics from the non-null cells of a GRASS raster map
+ *
+ * Copyright (C) 2004-2006, 2012 by the GRASS Development Team
+ * Author(s): Hamish Bowman, University of Otago, New Zealand
+ * Extended stats: Martin Landa
+ * Zonal stats: Markus Metz
+ *
+ * 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>
+#include "globals.h"
+#include "../../vector/v.vol.rst/userglobs.h"
+
+param_type param;
+zone_type zone_info;
+
+/* ************************************************************************* */
+/* Set up the arguments we are expecting ********************************** */
+/* ************************************************************************* */
+void set_params()
+{
+ param.inputfile = G_define_standard_option(G_OPT_R_MAPS);
+
+ param.zonefile = G_define_standard_option(G_OPT_R_MAP);
+ param.zonefile->key = "zones";
+ param.zonefile->required = NO;
+ 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.output_file->guisection = _("Output settings");
+
+ param.percentile = G_define_option();
+ param.percentile->key = "percentile";
+ param.percentile->type = TYPE_DOUBLE;
+ 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.percentile->guisection = _("Extended");
+
+ param.tolerance = G_define_option();
+ param.tolerance->key = "tolerance";
+ param.tolerance->type = TYPE_DOUBLE;
+ param.tolerance->required = NO;
+ param.tolerance->multiple = NO;
+ param.tolerance->answer = "0.5";
+ param.tolerance->description =
+ _("Tolerance to consider float number equal to another when computing the mode");
+ param.tolerance->guisection = _("Extended");
+
+ param.separator = G_define_standard_option(G_OPT_F_SEP);
+ param.separator->guisection = _("Formatting");
+
+ param.shell_style = G_define_flag();
+ param.shell_style->key = 'g';
+ param.shell_style->description =
+ _("Print the stats in shell script style");
+ param.shell_style->guisection = _("Formatting");
+
+ param.extended = G_define_flag();
+ param.extended->key = 'e';
+ param.extended->description = _("Calculate extended statistics");
+ param.extended->guisection = _("Extended");
+
+ param.table = G_define_flag();
+ param.table->key = 't';
+ param.table->description = _("Table output format instead of standard output format");
+ param.table->guisection = _("Formatting");
+
+ return;
+}
+
+static int open_raster(const char *infile);
+static univar_stat *univar_stat_with_percentiles();
+static void process_raster(univar_stat * stats,
+ int fd, int fdz,
+ const struct Cell_head *region);
+static void init_zones(int fdz, const struct Cell_head *region);
+
+
+void init_zones(int fdz, const struct Cell_head *region)
+{
+ const unsigned int rows = region->rows;
+ const unsigned int cols = region->cols;
+ unsigned int row;
+ unsigned int col;
+ int zone = 0, n_zones = zone_info.n_zones;
+
+ CELL *zoneraster_row;
+ CELL *zptr;
+
+ zone_info.n_alloc = 0;
+ zone_info.len = (int *) G_calloc(zone_info.n_zones, sizeof(int));
+
+ zoneraster_row = Rast_allocate_c_buf();
+ G_message("Reading the zones map.");
+ for (row = 0; row < rows; row++) {
+ Rast_get_c_row(fdz, zoneraster_row, row);
+ zptr = zoneraster_row;
+
+ for (col = 0; col < cols; col++) {
+ if (Rast_is_c_null_value(zptr)) {
+ zptr++;
+ continue;
+ }
+ zone = *zptr - zone_info.min;
+ zone_info.len[zone] += 1;
+ zptr++;
+ }
+ G_percent(row, rows, 2);
+ }
+ /* find the max lenght to understand how much memory nust be allocated */
+ for (zone = 0; zone < n_zones; zone++) {
+ if (zone_info.len[zone] > zone_info.n_alloc)
+ zone_info.n_alloc = zone_info.len[zone];
+ }
+ return;
+}
+
+
+
+
+/* *************************************************************** */
+/* **** the main functions for r.univar ************************** */
+/* *************************************************************** */
+int main(int argc, char *argv[])
+{
+ 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;
+ const char *mapset, *name;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("statistics"));
+ G_add_keyword(_("univariate statistics"));
+ G_add_keyword(_("zonal 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);
+
+ /* table field separator */
+ zone_info.sep = param.separator->answer;
+ if (strcmp(zone_info.sep, "\\t") == 0)
+ zone_info.sep = "\t";
+ if (strcmp(zone_info.sep, "tab") == 0)
+ zone_info.sep = "\t";
+ if (strcmp(zone_info.sep, "space") == 0)
+ zone_info.sep = " ";
+ if (strcmp(zone_info.sep, "comma") == 0)
+ zone_info.sep = ",";
+
+ 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;
+
+ fdz = -1;
+
+ /* open zoning raster */
+ z = param.zonefile->answer;
+ mapset = G_find_raster2(z, "");
+
+ fdz = open_raster(z);
+
+ cell_type = Rast_get_map_type(fdz);
+ if (cell_type != CELL_TYPE)
+ G_fatal_error("Zoning raster must be of type CELL");
+
+ if (Rast_read_range(z, mapset, &zone_range) == -1)
+ G_fatal_error("Can not read range for zoning raster");
+ Rast_get_range_min_max(&zone_range, &min, &max);
+ if (Rast_read_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;
+
+ init_zones(fdz, ®ion);
+
+ sscanf(param.tolerance->answers[0], "%lf", ¶m.tol);
+
+ /* count the input rasters given */
+ for (p = (char **)param.inputfile->answers, rasters = 0;
+ *p; p++, rasters++) ;
+
+ /* process all input rasters */
+ int map_type = param.extended->answer ? -2 : -1;
+
+ stats = ((map_type == -1)
+ ? create_univar_stat_struct(-1, 0)
+ : 0);
+
+ for (p = param.inputfile->answers; *p; p++) {
+ fd = open_raster(*p);
+
+ if (map_type != -1) {
+ /* NB: map_type must match when doing extended stats */
+ int this_type = Rast_get_map_type(fd);
+
+ assert(this_type > -1);
+ if (map_type < -1) {
+ /* extended stats */
+ assert(stats == 0);
+ map_type = this_type;
+ stats = univar_stat_with_percentiles(map_type);
+ }
+ else if (this_type != map_type) {
+ G_fatal_error(_("Raster <%s> type mismatch"), *p);
+ }
+ }
+ process_raster(stats, fd, fdz, ®ion);
+
+ /* close input raster */
+ Rast_close(fd);
+ }
+
+ /* close zoning raster */
+ if (z)
+ Rast_close(fdz);
+
+ /* create the output */
+ if (param.table->answer)
+ print_stats_table(stats);
+ /* else
+ print_stats(stats); */
+
+ /* release memory */
+ free_univar_stat_struct(stats);
+
+ exit(EXIT_SUCCESS);
+}
+
+static int open_raster(const char *infile)
+{
+ const char *mapset;
+ int fd;
+
+ mapset = G_find_raster2(infile, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Raster map <%s> not found"), infile);
+ }
+
+ fd = Rast_open_old(infile, mapset);
+
+ return fd;
+}
+
+static univar_stat *univar_stat_with_percentiles()
+{
+ univar_stat *stats;
+ unsigned int z, p, n_perc=0;
+ unsigned int n_zones = zone_info.n_zones;
+
+ while (param.percentile->answers[n_perc])
+ n_perc++;
+
+ /* allocate memory */
+ G_debug(3, "Allocate memory for percentile, n_perc=%d", n_perc);
+ param.n_perc = n_perc;
+ param.index_perc = (int *) G_calloc(n_perc + 1, sizeof(int));
+ param.quant_perc = (double *) G_calloc(n_perc + 1, sizeof(double));
+ param.perc = (double *) G_calloc(n_perc + 1, sizeof(double));
+ for (p = 0; p < n_perc; p++) {
+ G_debug(3, "Percentile: %s", param.percentile->answers[p]);
+ sscanf(param.percentile->answers[p], "%lf", &(param.perc[p]));
+ }
+
+ stats = create_univar_stat_struct();
+ for (z = 0; z < n_zones; z++) {
+ stats[z].size = zone_info.len[z];
+ stats[z].zone = z + zone_info.min;
+ stats[z].cat = Rast_get_c_cat(&stats[z].zone, &(zone_info.cats));
+ stats[z].perc = (double *) G_calloc(n_perc, sizeof(double));
+ }
+ 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 unsigned int rows = region->rows;
+ const unsigned int cols = region->cols;
+
+ const RASTER_MAP_TYPE map_type = Rast_get_map_type(fd);
+ const size_t value_sz = Rast_cell_size(map_type);
+ unsigned int row;
+ void *raster_row;
+ CELL *zoneraster_row;
+
+ raster_row = Rast_allocate_buf(map_type);
+ zoneraster_row = Rast_allocate_c_buf();
+
+ for (row = 0; row < rows; row++) {
+ void *ptr;
+ CELL *zptr;
+ unsigned int col;
+
+ Rast_get_row(fd, raster_row, row, map_type);
+ Rast_get_c_row(fdz, zoneraster_row, row);
+ zptr = zoneraster_row;
+
+ ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ double val;
+ int zone = 0;
+
+ if (Rast_is_c_null_value(zptr)) {
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ continue;
+ }
+ zone = *zptr - zone_info.min;
+ /* count all including NULL cells in input map */
+ /*stats[zone].size++;
+ G _debug(3, "Col: %u, zone: %d", col, zone); */
+
+ /* can't do stats with NULL cells in input map */
+ if (Rast_is_null_value(ptr, map_type)) {
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ continue;
+ }
+
+ if (param.extended->answer && stats[zone].array == NULL){
+ stats[zone].n_alloc = zone_info.len[zone] + 1;
+ stats[zone].array = (DCELL *) G_calloc(stats[zone].n_alloc,
+ sizeof(DCELL));
+ }
+
+ val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr)
+ : (map_type == FCELL_TYPE) ? *((FCELL *) ptr)
+ : *((CELL *) ptr));
+
+ compute_stats(&stats[zone], val, zone_info.len[zone]);
+
+ ptr = G_incr_void_ptr(ptr, value_sz);
+ zptr++;
+ }
+ G_percent(row, rows, 2);
+ }
+ G_percent(rows, rows, 2);
+ return;
+}
Added: grass-addons/grass7/raster/r.univar2/sort.c
===================================================================
--- grass-addons/grass7/raster/r.univar2/sort.c (rev 0)
+++ grass-addons/grass7/raster/r.univar2/sort.c 2013-12-10 22:25:27 UTC (rev 58439)
@@ -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/grass7/raster/r.univar2/stats.c
===================================================================
--- grass-addons/grass7/raster/r.univar2/stats.c (rev 0)
+++ grass-addons/grass7/raster/r.univar2/stats.c 2013-12-10 22:25:27 UTC (rev 58439)
@@ -0,0 +1,415 @@
+/*
+ * 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"
+/*
+#include "../../lib/raster/rasterlib.dox"
+*/
+#define I 0
+#define J 1
+
+/* *************************************************************** */
+/* **** univar_stat constructor ********************************** */
+/* *************************************************************** */
+univar_stat *create_univar_stat_struct()
+{
+ univar_stat *stats;
+ int z;
+ int n_zones = zone_info.n_zones;
+
+ if (n_zones == 0)
+ n_zones = 1;
+
+ stats = (univar_stat *) G_calloc(n_zones, sizeof(univar_stat));
+
+ for (z = 0; z < n_zones; z++) {
+ stats[z].zone = 0;
+ stats[z].cat = NULL;
+ stats[z].null_cells = 0;
+ stats[z].sum = 0.0;
+ stats[z].sum2 = 0.0;
+ stats[z].sum3 = 0.0;
+ stats[z].sum4 = 0.0;
+ stats[z].sum_abs = 0.0;
+ stats[z].min = 0.0 / 0.0; /* set to nan as default */
+ stats[z].max = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].range = 0.0 / 0.0; /* set to nan as default */
+ stats[z].mean = 0.0 / 0.0; /* set to nan as default */
+ stats[z].meandev = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].variance = 0.0 / 0.0; /* set to nan as default */
+ stats[z].stddev = 0.0 / 0.0; /* set to nan as default */
+ stats[z].var_coef = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].skewness = 0.0 / 0.0; /* set to nan as default */
+ stats[z].kurtosis = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].variance2 = 0.0 / 0.0; /* set to nan as default */
+ stats[z].stddev2 = 0.0 / 0.0; /* set to nan as default */
+ stats[z].var_coef2 = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].skewness2 = 0.0 / 0.0; /* set to nan as default */
+ stats[z].kurtosis2 = 0.0 / 0.0; /* set to nan as default */
+
+ stats[z].quartile_25 = 0.0 / 0.0; /* set to nan as default */
+ stats[z].median = 0.0 / 0.0; /* set to nan as default */
+ stats[z].quartile_75 = 0.0 / 0.0; /* set to nan as default */
+ stats[z].mode = 0.0 / 0.0; /* set to nan as default */
+ stats[z].occurrences = 0;
+
+ stats[z].n = 0;
+ stats[z].size = 0;
+ stats[z].map_type = 0;
+ stats[z].array = NULL;
+ stats[z].nextp = NULL;
+ stats[z].n_alloc = 0;
+ }
+
+ return stats;
+}
+
+
+/* *************************************************************** */
+/* **** univar_stat destructor *********************************** */
+/* *************************************************************** */
+void free_univar_stat_struct(univar_stat * stats)
+{
+ /*int z, n_zones = zone_info.n_zones;
+
+ for (z = 0; z < n_zones; z++){
+ if (stats[z].perc != NULL)
+ G_free(&stats[z].perc);
+ } */
+
+ return;
+}
+
+int factorial(int n)
+{
+ int res = 1;
+
+ while (n>1)
+ res *= n--;
+ return res;
+}
+
+int num_combinations(int n, int k)
+{
+ return factorial(n) / (factorial(k) * factorial(n - k));
+}
+
+void index_combinations(int n, int **indexes)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++)
+ for (j = i + 1; j < n; j++)
+ indexes[i][I] = i;
+ indexes[i][J] = j;
+ return;
+}
+
+
+void sort_mode_double(double *array, int n, double tol,
+ double *mode, int *occurrences)
+{
+ int previous = array[0];
+ int i = 1, counter = 1;
+ *mode = (double) array[0];
+ *occurrences = 1;
+
+ if (n > 1)
+ {
+ while (i < n)
+ {
+ if ((double) fabs(array[i] - previous) > tol)
+ {
+ if (counter > *occurrences)
+ {
+ *occurrences = counter;
+ *mode = (double) previous;
+ }
+ previous = array[i];
+ counter = 1;
+ } else
+ {
+ counter += 1;
+ }
+ i += 1;
+ }
+ }
+ return;
+}
+
+
+void mode_double(double *array, int n, double tol,
+ double *mode, int *occurrences, int **indexes)
+{
+ int i, i_previous = 0, counter = 0, num_comb = 0;
+ *mode = 0;
+ *occurrences = 0;
+
+ num_comb = num_combinations(n, 2);
+ index_combinations(n, indexes);
+ for (i = 0; i < num_comb; i++)
+ if (array[indexes[i][I]] != i_previous)
+ {
+ if (counter > *occurrences)
+ {
+ *occurrences = counter;
+ *mode = (double) i_previous;
+ }
+ i_previous = array[indexes[i][I]];
+ counter = 0;
+ }
+ if (fabs((double) array[indexes[i][I]] - (double) array[indexes[i][J]]) < tol)
+ counter += 1;
+ return;
+}
+
+
+int stats_extend(univar_stat *stat){
+ int p, qind_25, qind_75;
+ int n = stat->n;
+
+ for (p = 0; p < param.n_perc; p++) {
+ param.index_perc[p] = (int)(n * 1e-2 * param.perc[p] - 0.5);
+ stat->perc[p] = stat->array[param.index_perc[p]];
+ }
+ qind_25 = (int)(n * 0.25 - 0.5);
+ qind_75 = (int)(n * 0.75 - 0.5);
+
+ heapsort_double(stat->array, n);
+
+ stat->quartile_25 = stat->array[qind_25];
+ /* odd ? odd : even */
+ stat->median = ((n % 2)?
+ (stat->array[(int)(n/2)]) :
+ (stat->array[n/2 - 1] + stat->array[n/2]) / 2.0);
+ stat->quartile_75 = stat->array[qind_75];
+
+ sort_mode_double(stat->array, n, param.tol,
+ &stat->mode, &stat->occurrences);
+
+ /* free the memory after compute the extended statistics */
+ G_free(stat->array);
+ stat->array = NULL;
+ stat->n_alloc = 0;
+ return 0;
+}
+
+
+int stats_general(univar_stat *s){
+ double n = s->n;
+
+ s->null_cells = s->size - s->n;
+ s->range = s->max -s->min;
+ s->mean = s->sum / n;
+ s->meandev = s->sum_abs / n;
+
+ s->variance = (s->sum2 - s->sum * s->sum / n) / (n - 1);
+ if (s->variance < GRASS_EPSILON)
+ s->variance = 0.0;
+ s->stddev = sqrt(s->variance);
+ s->var_coef = (s->stddev / s->mean) * 100.;
+
+ s->skewness = (s->sum3 / n
+ - 3 * s->sum * s->sum2 / (n * n)
+ + 2 * s->sum * s->sum * s->sum / (n * n * n)
+ / (pow(s->variance, 1.5)));
+ s->kurtosis = ((s->sum4 / n
+ - 4 * s->sum * s->sum3 / (n * n)
+ + 6 * s->sum * s->sum * s->sum2 / (n * n * n)
+ - 3 * s->sum * s->sum * s->sum* s->sum
+ / (n * n * n * n))
+ / (s->variance * s->variance) - 3);
+
+ s->variance2 = s->sum2 / (n - 1);
+ if (s->variance2 < GRASS_EPSILON)
+ s->variance2 = 0.0;
+ s->stddev2 = sqrt(s->variance2);
+ s->var_coef2 = (s->stddev2 / s->mean) * 100. ;
+
+ s->skewness2 = (s->sum3 / (s->stddev2 * s->stddev2 * s->stddev2) / n);
+ s->kurtosis2 = (s->sum4 / (s->variance2 * s->variance2) / (n - 3));
+ return 0;
+}
+
+
+int compute_stats(univar_stat *stat, double val, int len){
+ stat->sum += val;
+ stat->sum2 += val * val;
+ stat->sum3 += val * val * val;
+ stat->sum4 += val * val * val * val;
+ stat->sum_abs += fabs(val);
+ stat->min = (isnan(stat->min) || val < stat->min)? val: stat->min;
+ stat->max = (isnan(stat->max) || val > stat->max)? val: stat->max;
+
+ if (stat->array != NULL)
+ G_debug(3, "Compute_stats, zone: %d, val=%f, %d/%d, n_alloc:%d", stat->zone, val, stat->n, stat->size, stat->n_alloc);
+ stat->array[stat->n] = val;
+
+ stat->n++;
+
+
+ if (stat->n == stat->size){
+ G_debug(3, " Finish the zone: %d, sum=%f", stat->zone, stat->sum);
+ /* finish this zone */
+ stats_general(stat);
+ if (stat->array != NULL){
+ stats_extend(stat);
+ return 1;
+ }
+ }
+ return 0;
+}
+
+
+/* *************************************************************** */
+/* **** compute and print univar statistics to stdout ************ */
+/* *************************************************************** */
+
+void print_cols_table(){
+ int i;
+
+ fprintf(stdout, "zone");
+ fprintf(stdout, "%slabel", zone_info.sep);
+ fprintf(stdout, "%sall_cells", zone_info.sep);
+ fprintf(stdout, "%snon_null_cells", zone_info.sep);
+ fprintf(stdout, "%snull_cells", zone_info.sep);
+
+ fprintf(stdout, "%ssum", zone_info.sep);
+ fprintf(stdout, "%ssum_abs", zone_info.sep);
+ fprintf(stdout, "%smin", zone_info.sep);
+ fprintf(stdout, "%smax", zone_info.sep);
+
+ fprintf(stdout, "%srange", zone_info.sep);
+ fprintf(stdout, "%smean", zone_info.sep);
+ fprintf(stdout, "%smean_of_abs", zone_info.sep);
+
+ fprintf(stdout, "%svariance", zone_info.sep);
+ fprintf(stdout, "%sstddev", zone_info.sep);
+ fprintf(stdout, "%scoeff_var", zone_info.sep);
+
+ fprintf(stdout, "%sskewness", zone_info.sep);
+ fprintf(stdout, "%skurtosis", zone_info.sep);
+
+ fprintf(stdout, "%svariance2", zone_info.sep);
+ fprintf(stdout, "%sstddev2", zone_info.sep);
+ fprintf(stdout, "%scoeff_var2", zone_info.sep);
+
+ fprintf(stdout, "%sskewness2", zone_info.sep);
+ fprintf(stdout, "%skurtosis2", zone_info.sep);
+
+ if (param.extended->answer) {
+ fprintf(stdout, "%sfirst_quart", zone_info.sep);
+ fprintf(stdout, "%smedian", zone_info.sep);
+ fprintf(stdout, "%sthird_quart", zone_info.sep);
+ for (i = 0; i < param.n_perc; i++) {
+
+ if (param.perc[i] == (int) param.perc[i]) {
+ /* percentile is an exact integer */
+ fprintf(stdout, "%sperc_%d", zone_info.sep, (int)param.perc[i]);
+ }
+ else {
+ /* percentile is not an exact integer */
+ char buf[24];
+ sprintf(buf, "%.15g", param.perc[i]);
+ G_strchg(buf, '.', '_');
+ fprintf(stdout, "%sperc_%s", zone_info.sep, buf);
+ }
+ }
+
+ fprintf(stdout, "%smode", zone_info.sep);
+ fprintf(stdout, "%soccurrences", zone_info.sep);
+ }
+ fprintf(stdout, "\n");
+ return;
+}
+
+int print_row(univar_stat *stat, char *sep){
+ int i;
+ char sum_str[100];
+
+ /* zone number */
+ fprintf(stdout, "%d", stat->zone);
+ /* zone label */
+ fprintf(stdout,"%s%s", sep, stat->cat);
+ /* zone all_cells */
+ fprintf(stdout,"%s%d", sep, stat->size);
+ /* non-null cells */
+ fprintf(stdout, "%s%d", sep, stat->n);
+ /* null cells */
+ fprintf(stdout, "%s%d", sep, stat->null_cells);
+
+
+ sprintf(sum_str, "%.15g", stat->sum);
+ G_trim_decimal(sum_str);
+ fprintf(stdout, "%s%s", sep, sum_str);
+ sprintf(sum_str, "%.15g", stat->sum_abs);
+ G_trim_decimal(sum_str);
+ fprintf(stdout, "%s%s", sep, sum_str);
+ fprintf(stdout, "%s%.15g", sep, stat->min);
+ fprintf(stdout, "%s%.15g", sep, stat->max);
+
+
+ fprintf(stdout, "%s%.15g", sep, stat->range);
+ fprintf(stdout, "%s%.15g", sep, stat->mean);
+ fprintf(stdout, "%s%.15g", sep, stat->meandev);
+
+ fprintf(stdout, "%s%.15g", sep, stat->variance);
+ fprintf(stdout, "%s%.15g", sep, stat->stddev);
+ fprintf(stdout, "%s%.15g", sep, stat->var_coef);
+
+ fprintf(stdout, "%s%g", sep, stat->skewness);
+ fprintf(stdout, "%s%g", sep, stat->kurtosis);
+
+
+ fprintf(stdout, "%s%.15g", sep, stat->variance2);
+ fprintf(stdout, "%s%.15g", sep, stat->stddev2);
+ fprintf(stdout, "%s%.15g", sep, stat->var_coef2);
+
+ fprintf(stdout, "%s%g", sep, stat->skewness2);
+ fprintf(stdout, "%s%g", sep, stat->kurtosis2);
+
+ if (param.extended->answer){
+ fprintf(stdout, "%s%g", sep, stat->quartile_25);
+ fprintf(stdout, "%s%g", sep, stat->median);
+ fprintf(stdout, "%s%g", sep, stat->quartile_75);
+ for (i = 0; i < param.n_perc; i++) {
+ fprintf(stdout, "%s%g", sep, stat->perc[i]);
+ }
+
+ fprintf(stdout, "%s%g", sep, stat->mode);
+ fprintf(stdout, "%s%d", sep, stat->occurrences);
+ }
+ fprintf(stdout, "\n");
+ return 0;
+}
+
+int print_stats_table(univar_stat * stats)
+{
+ int z, n_zones = zone_info.n_zones;
+
+ /* print column headers */
+ print_cols_table();
+
+ /* print stats */
+ for (z = 0; z < n_zones; z++) {
+ /* stats collected for this zone? */
+ if (!stats[z].n == 0)
+ print_row(&stats[z], zone_info.sep);
+ }
+ return 1;
+}
More information about the grass-commit
mailing list