[GRASS-SVN] r33926 - in grass/trunk/raster: . r.statistics3
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Oct 17 12:37:26 EDT 2008
Author: glynn
Date: 2008-10-17 12:37:26 -0400 (Fri, 17 Oct 2008)
New Revision: 33926
Added:
grass/trunk/raster/r.statistics3/
grass/trunk/raster/r.statistics3/Makefile
grass/trunk/raster/r.statistics3/main.c
grass/trunk/raster/r.statistics3/r.statistics3.html
Log:
Add r.statistics3 (quantiles)
Property changes on: grass/trunk/raster/r.statistics3
___________________________________________________________________
Name: svn:ignore
+ *OBJ*
Added: grass/trunk/raster/r.statistics3/Makefile
===================================================================
--- grass/trunk/raster/r.statistics3/Makefile (rev 0)
+++ grass/trunk/raster/r.statistics3/Makefile 2008-10-17 16:37:26 UTC (rev 33926)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.statistics3
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass/trunk/raster/r.statistics3/main.c
===================================================================
--- grass/trunk/raster/r.statistics3/main.c (rev 0)
+++ grass/trunk/raster/r.statistics3/main.c 2008-10-17 16:37:26 UTC (rev 33926)
@@ -0,0 +1,525 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.statistics3
+ * AUTHOR(S): Glynn Clements <glynn gclements.plus.com> (original contributor),
+ * PURPOSE: Compute category or object oriented quantiles using two passes
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/spawn.h>
+
+#define MAX_CATS 1000
+
+struct bin
+{
+ unsigned long origin;
+ DCELL min, max;
+ int base, count;
+};
+
+struct basecat
+{
+ unsigned int *slots;
+ unsigned long total;
+ int num_values;
+ unsigned char *slot_bins;
+ int num_bins;
+ struct bin *bins;
+ DCELL *values;
+ DCELL *quants;
+};
+
+static int num_quants;
+static DCELL *quants;
+static DCELL fmin, fmax;
+static int num_slots;
+static DCELL slot_size;
+
+static int rows, cols;
+
+static int min, max;
+static int num_cats;
+static struct basecat *basecats;
+
+static inline int get_slot(DCELL c)
+{
+ int i = (int)floor((c - min) / slot_size);
+
+ if (i < 0)
+ i = 0;
+ if (i > num_slots - 1)
+ i = num_slots - 1;
+ return i;
+}
+
+static inline double get_quantile(struct basecat *bc, int n)
+{
+ return (double)bc->total * quants[n];
+}
+
+static void get_slot_counts(int basefile, int coverfile)
+{
+ CELL *basebuf = G_allocate_c_raster_buf();
+ DCELL *coverbuf = G_allocate_d_raster_buf();
+ int row, col;
+
+ G_message(_("Computing histograms"));
+
+ for (row = 0; row < rows; row++) {
+ G_get_c_raster_row(basefile, basebuf, row);
+ G_get_d_raster_row(coverfile, coverbuf, row);
+
+ for (col = 0; col < cols; col++) {
+ struct basecat *bc;
+ int i;
+
+ if (G_is_c_null_value(&basebuf[col]))
+ continue;
+
+ if (G_is_d_null_value(&coverbuf[col]))
+ continue;
+
+ i = get_slot(coverbuf[col]);
+ bc = &basecats[basebuf[col] - min];
+
+ bc->slots[i]++;
+ bc->total++;
+ }
+
+ G_percent(row, rows, 2);
+ }
+
+ G_percent(rows, rows, 2);
+ G_free(basebuf);
+ G_free(coverbuf);
+}
+
+static void initialize_bins(void)
+{
+ int cat;
+
+ G_message(_("Computing bins"));
+
+ for (cat = 0; cat < num_cats; cat++) {
+ struct basecat *bc = &basecats[cat];
+ int slot;
+ double next;
+ int num_values = 0;
+ int bin = 0;
+ unsigned long accum = 0;
+ int quant = 0;
+
+ bc->bins = G_calloc(num_quants, sizeof(struct bin));
+ bc->slot_bins = G_calloc(num_slots, sizeof(unsigned char));
+
+ next = get_quantile(bc, quant);
+
+ for (slot = 0; slot < num_slots; slot++) {
+ unsigned int count = bc->slots[slot];
+ unsigned long accum2 = accum + count;
+
+ if (accum2 > next) {
+ struct bin *b = &bc->bins[bin];
+
+ bc->slot_bins[slot] = ++bin;
+
+ b->origin = accum;
+ b->base = num_values;
+ b->count = 0;
+ b->min = min + slot_size * slot;
+ b->max = min + slot_size * (slot + 1);
+
+ while (accum2 > next)
+ next = get_quantile(bc, ++quant);
+
+ num_values += count;
+ }
+
+ accum = accum2;
+ }
+
+ bc->num_values = num_values;
+ bc->num_bins = bin;
+
+ G_free(bc->slots);
+
+ bc->values = G_calloc(num_values, sizeof(DCELL));
+ }
+}
+
+static void fill_bins(int basefile, int coverfile)
+{
+ CELL *basebuf = G_allocate_c_raster_buf();
+ DCELL *coverbuf = G_allocate_d_raster_buf();
+ int row, col;
+
+ G_message(_("Binning data"));
+
+ for (row = 0; row < rows; row++) {
+ G_get_c_raster_row(basefile, basebuf, row);
+ G_get_d_raster_row(coverfile, coverbuf, row);
+
+ for (col = 0; col < cols; col++) {
+ struct basecat *bc;
+ int i, bin;
+ struct bin *b;
+
+ if (G_is_c_null_value(&basebuf[col]))
+ continue;
+
+ if (G_is_d_null_value(&coverbuf[col]))
+ continue;
+
+ i = get_slot(coverbuf[col]);
+ bc = &basecats[basebuf[col] - min];
+ if (!bc->slot_bins[i])
+ continue;
+
+ bin = bc->slot_bins[i] - 1;
+ b = &bc->bins[bin];
+
+ bc->values[b->base + b->count++] = coverbuf[col];
+ }
+
+ G_percent(row, rows, 2);
+ }
+
+ G_percent(rows, rows, 2);
+ G_free(basebuf);
+ G_free(coverbuf);
+}
+
+static int compare_dcell(const void *aa, const void *bb)
+{
+ DCELL a = *(const DCELL *)aa;
+ DCELL b = *(const DCELL *)bb;
+
+ if (a < b)
+ return -1;
+ if (a > b)
+ return 1;
+ return 0;
+}
+
+static void sort_bins(void)
+{
+ int cat;
+
+ G_message(_("Sorting bins"));
+
+ for (cat = 0; cat < num_cats; cat++) {
+ struct basecat *bc = &basecats[cat];
+ int bin;
+
+ G_free(bc->slot_bins);
+
+ for (bin = 0; bin < bc->num_bins; bin++) {
+ struct bin *b = &bc->bins[bin];
+
+ qsort(&bc->values[b->base], b->count, sizeof(DCELL), compare_dcell);
+ }
+
+ G_percent(cat, num_cats, 2);
+ }
+
+ G_percent(cat, num_cats, 2);
+}
+
+static void print_quantiles(void)
+{
+ int cat;
+
+ G_message(_("Printing quantiles"));
+
+ for (cat = 0; cat < num_cats; cat++) {
+ struct basecat *bc = &basecats[cat];
+ int quant;
+
+ for (quant = 0; quant < num_quants; quant++)
+ printf("%d:%d:%f:%f\n", min + cat, quant, 100 * quants[quant], bc->quants[quant]);
+ }
+}
+
+static void compute_quantiles(void)
+{
+ int cat;
+
+ G_message(_("Computing quantiles"));
+
+ for (cat = 0; cat < num_cats; cat++) {
+ struct basecat *bc = &basecats[cat];
+ struct bin *b = &bc->bins[0];
+ int quant;
+
+ bc->quants = G_malloc(num_quants * sizeof(DCELL));
+
+ for (quant = 0; quant < num_quants; quant++) {
+ double next = get_quantile(bc, quant);
+ double k, v;
+ int i0, i1;
+
+ while (b->origin + b->count < next)
+ b++;
+
+ k = next - b->origin;
+ i0 = (int)floor(k);
+ i1 = (int)ceil(k);
+
+ v = (i0 == i1)
+ ? bc->values[b->base + i0]
+ : bc->values[b->base + i0] * (i1 - k) + bc->values[b->base + i1] * (k - i0);
+
+ bc->quants[quant] = v;
+ }
+ }
+}
+
+static void do_reclass(const char *basemap, char **outputs)
+{
+ const char *tempfile = G_tempfile();
+ char *input_arg = G_malloc(strlen(basemap) + 7);
+ char *rules_arg = G_malloc(strlen(tempfile) + 7);
+ int quant;
+
+ G_message(_("Generating reclass maps"));
+
+ sprintf(input_arg, "input=%s", basemap);
+ sprintf(rules_arg, "rules=%s", tempfile);
+
+ for (quant = 0; quant < num_quants; quant++) {
+ const char *output = outputs[quant];
+ char *output_arg = G_malloc(strlen(output) + 8);
+ FILE *fp;
+ int cat;
+
+ sprintf(output_arg, "output=%s", output);
+
+ fp = fopen(tempfile, "w");
+ if (!fp)
+ G_fatal_error(_("Unable to open temporary file"));
+
+ for (cat = 0; cat < num_cats; cat++)
+ fprintf(fp, "%d = %d %f\n", min + cat, min + cat, basecats[cat].quants[quant]);
+
+ fclose(fp);
+
+ G_spawn("r.reclass", "r.reclass", input_arg, output_arg, rules_arg, NULL);
+ }
+
+ remove(tempfile);
+}
+
+static void do_output(int base_fd, char **outputs, const char *covermap)
+{
+ int *out_fd = G_malloc(num_quants * sizeof(int));
+ CELL *base_buf = G_allocate_c_raster_buf();
+ DCELL *out_buf = G_allocate_d_raster_buf();
+ const char *mapset = G_mapset();
+ struct Colors colors;
+ int have_colors;
+ int quant;
+ int row, col;
+
+ G_message(_("Writing output maps"));
+
+ for (quant = 0; quant < num_quants; quant++) {
+ const char *output = outputs[quant];
+
+ out_fd[quant] = G_open_fp_cell_new(output);
+ if (out_fd[quant] < 0)
+ G_fatal_error(_("Unable to open output map <%s>"), output);
+ }
+
+ have_colors = G_read_colors(covermap, "", &colors) > 0;
+
+ for (row = 0; row < rows; row++) {
+ G_get_c_raster_row(base_fd, base_buf, row);
+
+ for (quant = 0; quant < num_quants; quant++) {
+ for (col = 0; col < cols; col++)
+ if (G_is_c_null_value(&base_buf[col]))
+ G_set_d_null_value(&out_buf[col], 1);
+ else
+ out_buf[col] = basecats[base_buf[col] - min].quants[quant];
+
+ G_put_d_raster_row(out_fd[quant], out_buf);
+ }
+
+ G_percent(row, rows, 2);
+ }
+
+ G_percent(row, rows, 2);
+
+ for (quant = 0; quant < num_quants; quant++) {
+ G_close_cell(out_fd[quant]);
+ if (have_colors)
+ G_write_colors(outputs[quant], mapset, &colors);
+ }
+}
+
+int main(int argc, char *argv[])
+{
+ struct GModule *module;
+ struct
+ {
+ struct Option *quant, *perc, *slots, *basemap, *covermap, *output;
+ } opt;
+ struct {
+ struct Flag *r, *p;
+ } flag;
+ const char *basemap, *covermap;
+ char **outputs;
+ int reclass, print;
+ int cover_fd, base_fd;
+ struct Range range;
+ struct FPRange fprange;
+ int i;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster, statistics");
+ module->description = _("Compute category quantiles using two passes.");
+
+ opt.basemap = G_define_standard_option(G_OPT_R_BASE);
+
+ opt.covermap = G_define_standard_option(G_OPT_R_COVER);
+
+ opt.quant = G_define_option();
+ opt.quant->key = "quantiles";
+ opt.quant->type = TYPE_INTEGER;
+ opt.quant->required = NO;
+ opt.quant->description = _("Number of quantiles");
+
+ opt.perc = G_define_option();
+ opt.perc->key = "percentiles";
+ opt.perc->type = TYPE_DOUBLE;
+ opt.perc->multiple = YES;
+ opt.perc->description = _("List of percentiles");
+ opt.perc->answer = "50";
+
+ opt.slots = G_define_option();
+ opt.slots->key = "bins";
+ opt.slots->type = TYPE_INTEGER;
+ opt.slots->required = NO;
+ opt.slots->description = _("Number of bins to use");
+ opt.slots->answer = "1000";
+
+ opt.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.output->description = _("Resultant raster map(s)");
+ opt.output->required = NO;
+ opt.output->multiple = YES;
+
+ flag.r = G_define_flag();
+ flag.r->key = 'r';
+ flag.r->description =
+ _("Create reclass map with statistics as category labels");
+
+ flag.p = G_define_flag();
+ flag.p->key = 'p';
+ flag.p->description =
+ _("Don't create output maps; just print statistics");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ basemap = opt.basemap->answer;
+ covermap = opt.covermap->answer;
+ outputs = opt.output->answers;
+ reclass = flag.r->answer;
+ print = flag.p->answer;
+
+ if (!print && !opt.output->answers)
+ G_fatal_error(_("Either -p or output= must be given"));
+
+ if (print && opt.output->answers)
+ G_fatal_error(_("-p and output= are mutually exclusive"));
+
+ num_slots = atoi(opt.slots->answer);
+
+ if (opt.quant->answer) {
+ num_quants = atoi(opt.quant->answer) - 1;
+ quants = G_calloc(num_quants, sizeof(DCELL));
+ for (i = 0; i < num_quants; i++)
+ quants[i] = 1.0 * (i + 1) / (num_quants + 1);
+ }
+ else {
+ for (i = 0; opt.perc->answers[i]; i++)
+ ;
+ num_quants = i;
+ quants = G_calloc(num_quants, sizeof(DCELL));
+ for (i = 0; i < num_quants; i++)
+ quants[i] = atof(opt.perc->answers[i]) / 100;
+ qsort(quants, num_quants, sizeof(DCELL), compare_dcell);
+ }
+
+ if (opt.output->answer) {
+ for (i = 0; opt.output->answers[i]; i++)
+ ;
+ if (i != num_quants)
+ G_fatal_error(_("Number of quantiles (%d) does not match number of output maps (%d)"),
+ num_quants, i);
+ }
+
+ base_fd = G_open_cell_old(basemap, "");
+ if (base_fd < 0)
+ G_fatal_error(_("Unable to open base map <%s>"), basemap);
+
+ cover_fd = G_open_cell_old(covermap, "");
+ if (cover_fd < 0)
+ G_fatal_error(_("Unable to open cover map <%s>"), covermap);
+
+ if (G_raster_map_is_fp(basemap, "") != 0)
+ G_fatal_error(_("The base map must be an integer (CELL) map"));
+
+ if (G_read_range(basemap, "", &range) < 0)
+ G_fatal_error(_("Unable to read range of base map <%s>"), basemap);
+
+ G_get_range_min_max(&range, &min, &max);
+ num_cats = max - min + 1;
+ if (num_cats > MAX_CATS)
+ G_fatal_error(_("Base map <%s> has too many categories (max: %d)"),
+ basemap, MAX_CATS);
+
+ G_read_fp_range(covermap, "", &fprange);
+ G_get_fp_range_min_max(&fprange, &fmin, &fmax);
+ slot_size = (fmax - fmin) / num_slots;
+
+ basecats = G_calloc(num_cats, sizeof(struct basecat));
+
+ for (i = 0; i < num_cats; i++)
+ basecats[i].slots = G_calloc(num_slots, sizeof(unsigned int));
+
+ rows = G_window_rows();
+ cols = G_window_cols();
+
+ get_slot_counts(base_fd, cover_fd);
+
+ initialize_bins();
+ fill_bins(base_fd, cover_fd);
+
+
+ sort_bins();
+ compute_quantiles();
+
+ if (print)
+ print_quantiles();
+ else if (reclass)
+ do_reclass(basemap, outputs);
+ else
+ do_output(base_fd, outputs, covermap);
+
+ G_close_cell(cover_fd);
+ G_close_cell(base_fd);
+
+ return (EXIT_SUCCESS);
+}
Added: grass/trunk/raster/r.statistics3/r.statistics3.html
===================================================================
--- grass/trunk/raster/r.statistics3/r.statistics3.html (rev 0)
+++ grass/trunk/raster/r.statistics3/r.statistics3.html 2008-10-17 16:37:26 UTC (rev 33926)
@@ -0,0 +1,19 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.statistics3</em> is intended to be a partial replacement for
+r.statistics, with support for floating-point cover maps. It provides
+quantile calculations, which are absent from <em>r.statistics2</em>.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.statistics.html">r.quant</a>
+<a href="r.statistics2.html">r.quant</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Glynn Clements
+<p>
+<i>Last changed: $Date: 2007-04-06 21:49:26 +0200 (Fri, 06 Apr 2007) $</i>
More information about the grass-commit
mailing list