[GRASS-SVN] r42991 - in grass-addons/raster/r.pi: . r.pi.export
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Aug 4 06:53:56 EDT 2010
Author: wegmann
Date: 2010-08-04 10:53:56 +0000 (Wed, 04 Aug 2010)
New Revision: 42991
Added:
grass-addons/raster/r.pi/r.pi.export/
grass-addons/raster/r.pi/r.pi.export/Makefile
grass-addons/raster/r.pi/r.pi.export/description.html
grass-addons/raster/r.pi/r.pi.export/frag.c
grass-addons/raster/r.pi/r.pi.export/helpers.c
grass-addons/raster/r.pi/r.pi.export/local_proto.h
grass-addons/raster/r.pi/r.pi.export/main.c
grass-addons/raster/r.pi/r.pi.export/stat_method.c
Log:
new module for exporting raster values based on patches
Added: grass-addons/raster/r.pi/r.pi.export/Makefile
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/Makefile (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/Makefile 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.pi.export
+
+LIBES = $(STATSLIB) $(GISLIB)
+DEPENDENCIES = $(STATSDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/raster/r.pi/r.pi.export/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/raster/r.pi/r.pi.export/description.html
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/description.html (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/description.html 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,27 @@
+<H2>DESCRIPTION</H2>
+
+The <EM>r.pi.export</EM> module computes statistics of values within a patch and exports them to an ascii file. Additionally the raw patch values split up by patch ID are exported. The <EM>patch ID</EM>
+raster can be generated in order to visually link patch statistics and individual patches.
+
+<P>
+
+The program will be run non-interactively if the user specifies program arguments (see OPTIONS) on the command
+line. Alternately, the user can simply type <B>r.pi.export</B> on the command line, without program
+arguments. In this case, the user will be prompted for flag settings and parameter values.
+
+
+<H2>NOTES</H2>
+
+For export and import of patch values after e.g. R calculation <EM>r.pi.export</EM> and <EM>r.pi.import</EM> can be applied using the same raster (same extent, resolution is mandatory).
+
+
+
+<H2>SEE ALSO</H2>
+
+<EM><A HREF="r.pi.import/description.html">r.pi.import</A></EM><br>
+
+<H2>AUTHOR</H2>
+Programming: Elshad Shirinov<br>
+Scientific concept: Dr. Martin Wegmann <br>
+Department of Remote Sensing <br>Remote Sensing and Biodiversity Unit<br> University of Wuerzburg, Germany
+
Property changes on: grass-addons/raster/r.pi/r.pi.export/description.html
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/raster/r.pi/r.pi.export/frag.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/frag.c (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/frag.c 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,196 @@
+#include "local_proto.h"
+
+typedef struct
+{
+ int x, y;
+} Position;
+
+Coords *writeFrag(DCELL * flagbuf, Coords * actpos, int row, int col,
+ int nrows, int ncols, int nbr_cnt);
+int getNeighbors(Position * res, DCELL * flagbuf, int x, int y, int nx,
+ int ny, int nbr_cnt);
+
+int getNeighbors(Position * res, DCELL * flagbuf, int x, int y, int nx,
+ int ny, int nbr_cnt)
+{
+ int left, right, top, bottom;
+
+ int i, j;
+
+ int cnt = 0;
+
+ switch (nbr_cnt) {
+ case 4: /* von Neumann neighborhood */
+ if (x > 0 && !G_is_d_null_value(&flagbuf[y * nx + x - 1])) {
+ res[cnt].x = x - 1;
+ res[cnt].y = y;
+ cnt++;
+ }
+ if (y > 0 && !G_is_d_null_value(&flagbuf[(y - 1) * nx + x])) {
+ res[cnt].x = x;
+ res[cnt].y = y - 1;
+ cnt++;
+ }
+ if (x < nx - 1 && !G_is_d_null_value(&flagbuf[y * nx + x + 1])) {
+ res[cnt].x = x + 1;
+ res[cnt].y = y;
+ cnt++;
+ }
+ if (y < ny - 1 && !G_is_d_null_value(&flagbuf[(y + 1) * nx + x])) {
+ res[cnt].x = x;
+ res[cnt].y = y + 1;
+ cnt++;
+ }
+ break;
+ case 8: /* Moore neighborhood */
+ left = x > 0 ? x - 1 : 0;
+ top = y > 0 ? y - 1 : 0;
+ right = x < nx - 1 ? x + 1 : nx - 1;
+ bottom = y < ny - 1 ? y + 1 : ny - 1;
+ for (i = left; i <= right; i++) {
+ for (j = top; j <= bottom; j++) {
+ if (!(i == x && j == y) &&
+ !G_is_d_null_value(&flagbuf[j * nx + i])) {
+ res[cnt].x = i;
+ res[cnt].y = j;
+ cnt++;
+ }
+ }
+ }
+ break;
+ }
+
+ return cnt;
+}
+
+Coords *writeFrag(DCELL * flagbuf, Coords * actpos, int row, int col,
+ int nrows, int ncols, int nbr_cnt)
+{
+ int x, y, i;
+
+ Position *list = (Position *) G_malloc(nrows * ncols * sizeof(Position));
+
+ Position *first = list;
+
+ Position *last = list;
+
+ Position *nbr_list = (Position *) G_malloc(8 * sizeof(Position));
+
+ /* count neighbors */
+ int neighbors = 0;
+
+ if (col <= 0 || !G_is_d_null_value(&flagbuf[row * ncols + col - 1]))
+ neighbors++;
+ if (row <= 0 || !G_is_d_null_value(&flagbuf[(row - 1) * ncols + col]))
+ neighbors++;
+ if (col >= ncols - 1 ||
+ !G_is_d_null_value(&flagbuf[row * ncols + col + 1]))
+ neighbors++;
+ if (row >= nrows - 1 ||
+ !G_is_d_null_value(&flagbuf[(row + 1) * ncols + col]))
+ neighbors++;
+
+ /* write first cell */
+ actpos->x = col;
+ actpos->y = row;
+ actpos->neighbors = neighbors;
+ actpos->value = flagbuf[row * ncols + col];
+ actpos++;
+
+ /* write position to id map */
+ id_map[row * ncols + col] = fragcount - 1;
+
+ /* delete position from buffer */
+ G_set_d_null_value(&flagbuf[row * ncols + col], 1);
+
+ /* push position on fifo-list */
+ last->x = col;
+ last->y = row;
+ last++;
+
+ while (first < last) {
+ /* get position from fifo-list */
+ int r = first->y;
+
+ int c = first->x;
+
+ first++;
+
+ int left = c > 0 ? c - 1 : 0;
+
+ int top = r > 0 ? r - 1 : 0;
+
+ int right = c < ncols - 1 ? c + 1 : ncols - 1;
+
+ int bottom = r < nrows - 1 ? r + 1 : nrows - 1;
+
+ /* add neighbors to fifo-list */
+ int cnt =
+ getNeighbors(nbr_list, flagbuf, c, r, ncols, nrows, nbr_cnt);
+
+ for (i = 0; i < cnt; i++) {
+ x = nbr_list[i].x;
+ y = nbr_list[i].y;
+
+ /* add position to fifo-list */
+ last->x = x;
+ last->y = y;
+ last++;
+
+ /* count neighbors */
+ neighbors = 0;
+ if (x <= 0 || !G_is_d_null_value(&flagbuf[y * ncols + x - 1]))
+ neighbors++;
+ if (y <= 0 || !G_is_d_null_value(&flagbuf[(y - 1) * ncols + x]))
+ neighbors++;
+ if (x >= ncols - 1 ||
+ !G_is_d_null_value(&flagbuf[y * ncols + x + 1]))
+ neighbors++;
+ if (y >= nrows - 1 ||
+ !G_is_d_null_value(&flagbuf[(y + 1) * ncols + x]))
+ neighbors++;
+
+ /* set values */
+ actpos->x = x;
+ actpos->y = y;
+ actpos->neighbors = neighbors;
+ actpos->value = flagbuf[y * ncols + x];
+ actpos++;
+
+ /* write position to id map */
+ id_map[y * ncols + x] = fragcount - 1;
+
+ /* delete position from buffer */
+ G_set_d_null_value(&flagbuf[y * ncols + x], 1);
+ }
+ }
+
+ G_free(list);
+ G_free(nbr_list);
+ return actpos;
+}
+
+void writeFragments(DCELL * flagbuf, int nrows, int ncols, int nbr_cnt)
+{
+ int row, col, i;
+
+ Coords *p;
+
+ fragcount = 0;
+ Coords *actpos = fragments[0];
+
+ /* find fragments */
+ for (row = 0; row < nrows; row++) {
+ for (col = 0; col < ncols; col++) {
+ if (!G_is_d_null_value(&flagbuf[row * ncols + col])) {
+ fragcount++;
+
+ fragments[fragcount] =
+ writeFrag(flagbuf, fragments[fragcount - 1], row, col,
+ nrows, ncols, nbr_cnt);
+ }
+ }
+ }
+
+ return;
+}
Added: grass-addons/raster/r.pi/r.pi.export/helpers.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/helpers.c (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/helpers.c 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,95 @@
+#include "local_proto.h"
+
+int Round(double d)
+{
+ return d < 0 ? d - 0.5 : d + 0.5;
+}
+
+int Random(int max)
+{
+ return max <=
+ RAND_MAX ? rand() % max : floor((double)rand() /
+ (double)(RAND_MAX + 1) * max);
+}
+
+double Randomf()
+{
+ return ((double)rand()) / ((double)RAND_MAX);
+}
+
+void print_buffer(int *buffer, int sx, int sy)
+{
+ int x, y;
+
+ fprintf(stderr, "buffer:\n");
+ for (y = 0; y < sy; y++) {
+ for (x = 0; x < sx; x++) {
+ switch (buffer[x + y * sx]) {
+ case NOTHING:
+ fprintf(stderr, ".");
+ break;
+ default:
+ if (buffer[x + y * sx] < 0) {
+ fprintf(stderr, "%d", buffer[x + y * sx]);
+ }
+ else {
+ fprintf(stderr, "%d", buffer[x + y * sx]);
+ }
+ break;
+ }
+ }
+ fprintf(stderr, "\n");
+ }
+}
+
+void print_d_buffer(DCELL * buffer, int sx, int sy)
+{
+ int x, y;
+
+ fprintf(stderr, "buffer:\n");
+ for (y = 0; y < sy; y++) {
+ for (x = 0; x < sx; x++) {
+ fprintf(stderr, "%0.2f ", buffer[y * sx + x]);
+ }
+ fprintf(stderr, "\n");
+ }
+}
+
+void print_array(DCELL * buffer, int size)
+{
+ int i;
+
+ for (i = 0; i < size; i++) {
+ fprintf(stderr, "%0.2f ", buffer[i]);
+ }
+ fprintf(stderr, "\n");
+}
+
+void print_fragments()
+{
+ int f;
+
+ Coords *p;
+
+ for (f = 0; f < fragcount; f++) {
+ fprintf(stderr, "frag%d: ", f);
+ for (p = fragments[f]; p < fragments[f + 1]; p++) {
+ fprintf(stderr, "(%d,%d), n=%d, val=%0.2f;", p->x, p->y,
+ p->neighbors, p->value);
+ }
+ fprintf(stderr, "\n");
+ }
+}
+
+void print_map(double *map, int size)
+{
+ int x, y;
+
+ fprintf(stderr, "map:\n");
+ for (y = 0; y < size; y++) {
+ for (x = 0; x < size; x++) {
+ fprintf(stderr, "%0.0f ", map[x + y * size]);
+ }
+ fprintf(stderr, "\n");
+ }
+}
Added: grass-addons/raster/r.pi/r.pi.export/local_proto.h
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/local_proto.h (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/local_proto.h 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,88 @@
+#ifndef LOCAL_PROTO_H
+#define LOCAL_PROTO_H
+
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/stats.h>
+#include <math.h>
+#include <time.h>
+
+#define NOTHING 0
+
+#define RESOLUTION 10000
+
+#define MIN_DOUBLE -1000000
+#define MAX_DOUBLE 1000000
+
+typedef struct
+{
+ int x, y;
+ int neighbors;
+ DCELL value;
+} Coords;
+
+typedef struct
+{
+ int x, y;
+} Point;
+
+typedef struct
+{
+ int x, y;
+ int patch;
+} PatchPoint;
+
+typedef DCELL(*f_statmethod) (DCELL *, int);
+
+/* helpers.c */
+int Round(double d);
+
+int Random(int max);
+
+double Randomf();
+
+void print_buffer(int *buffer, int sx, int sy);
+
+void print_d_buffer(DCELL * buffer, int sx, int sy);
+
+void print_map(double *map, int size);
+
+void print_array(DCELL * buffer, int size);
+
+void print_fragments();
+
+/* frag.c */
+void writeFragments(DCELL * flagbuf, int nrows, int ncols, int nbr_cnt);
+
+/* stat_method.c */
+DCELL average(DCELL * vals, int count);
+
+DCELL variance(DCELL * vals, int count);
+
+DCELL std_deviat(DCELL * vals, int count);
+
+DCELL median(DCELL * vals, int count);
+
+DCELL min(DCELL * vals, int count);
+
+DCELL max(DCELL * vals, int count);
+
+/* global variables */
+int verbose;
+
+Coords **fragments;
+
+Coords *cells;
+
+int fragcount;
+
+int sx, sy;
+
+int *id_map;
+
+int *adj_matrix;
+
+#endif /* LOCAL_PROTO_H */
Added: grass-addons/raster/r.pi/r.pi.export/main.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/main.c (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/main.c 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,436 @@
+#include "local_proto.h"
+
+/*
+ Export of patch based values
+
+ by Elshad Shirinov.
+ */
+
+struct statmethod
+{
+ f_statmethod *method; /* routine to compute new value */
+ char *name; /* method name */
+ char *text; /* menu display - full description */
+ char *suffix; /* output suffix */
+};
+
+static struct statmethod statmethods[] = {
+ {average, "average", "average of patch values", "avg"},
+ {variance, "variance", "variance of patch values", "var"},
+ {std_deviat, "standard deviation", "standard deviation of patch values",
+ "dev"},
+ {median, "median", "median of patch values", "med"},
+ {min, "min", "minimum of patch values", "min"},
+ {max, "max", "maximum of patch values", "max"},
+ {0, 0, 0, 0}
+};
+
+int main(int argc, char *argv[])
+{
+ /* input */
+ char *oldname, *oldmapset;
+
+ /* id raster */
+ char *idname, *idmapset;
+
+ /* patch raster */
+ char *patchname, *patchmapset;
+
+ /* in and out file pointers */
+ int in_fd;
+
+ FILE *out_fp; /* ASCII - output */
+
+ int out_fd;
+
+ /* parameters */
+ int keyval;
+
+ int stats[256];
+
+ int stat_count;
+
+ int ratio_flag;
+
+ int neighb_count;
+
+ /* maps */
+ DCELL *map;
+
+ /* other parameters */
+ char *title;
+
+ /* helper variables */
+ int row, col;
+
+ int *result;
+
+ DCELL *d_res;
+
+ int i, j, n;
+
+ int x, y;
+
+ Coords *p;
+
+ char output_name[256];
+
+ char *str;
+
+ int method;
+
+ f_statmethod perform_method;
+
+ DCELL val;
+
+ DCELL *values;
+
+ int count;
+
+ int area;
+
+ RASTER_MAP_TYPE map_type;
+
+ struct Cell_head ch, window;
+
+ struct GModule *module;
+
+ struct
+ {
+ struct Option *input, *output;
+ struct Option *values, *id_rast;
+ struct Option *patch_rast, *stats;
+ struct Option *title;
+ } parm;
+
+ struct
+ {
+ struct Flag *adjacent, *quiet;
+ } flag;
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster");
+ module->description =
+ _("Performs statistical analysis and export of values within a patch");
+
+ parm.input = G_define_option();
+ parm.input->key = "input";
+ parm.input->type = TYPE_STRING;
+ parm.input->required = YES;
+ parm.input->gisprompt = "old,cell,raster";
+ parm.input->description = _("Name of existing raster file");
+
+ parm.output = G_define_option();
+ parm.output->key = "output";
+ parm.output->type = TYPE_STRING;
+ parm.output->required = YES;
+ parm.output->gisprompt = "new_file,file,output";
+ parm.output->description =
+ _("Name for the output ASCII-file with statistical values");
+
+ parm.values = G_define_option();
+ parm.values->key = "values";
+ parm.values->type = TYPE_STRING;
+ parm.values->required = NO;
+ parm.values->gisprompt = "new_file,file,output";
+ parm.values->description =
+ _("Name for the output ASCII-file with patch values");
+
+ parm.id_rast = G_define_option();
+ parm.id_rast->key = "id_raster";
+ parm.id_rast->type = TYPE_STRING;
+ parm.id_rast->required = NO;
+ parm.id_rast->gisprompt = "new,cell,raster";
+ parm.id_rast->description = _("Name for the id raster");
+
+ parm.patch_rast = G_define_option();
+ parm.patch_rast->key = "patch_raster";
+ parm.patch_rast->type = TYPE_STRING;
+ parm.patch_rast->required = NO;
+ parm.patch_rast->gisprompt = "new,cell,raster";
+ parm.patch_rast->description = _("Name for the binary patch raster");
+
+ parm.stats = G_define_option();
+ parm.stats->key = "stats";
+ parm.stats->type = TYPE_STRING;
+ parm.stats->required = YES;
+ parm.stats->multiple = YES;
+ str = parm.stats->options = G_malloc(1024);
+ for (n = 0; statmethods[n].name; n++) {
+ if (n)
+ strcat(str, ",");
+ else
+ *str = 0;
+ strcat(str, statmethods[n].name);
+ }
+ parm.stats->description =
+ _("Statistical method to perform on the values");
+
+ parm.title = G_define_option();
+ parm.title->key = "title";
+ parm.title->key_desc = "\"phrase\"";
+ parm.title->type = TYPE_STRING;
+ parm.title->required = NO;
+ parm.title->description = _("Title of the output raster file");
+
+ flag.adjacent = G_define_flag();
+ flag.adjacent->key = 'a';
+ flag.adjacent->description =
+ _("Set for 8 cell-neighbors. 4 cell-neighbors are default.");
+
+ flag.quiet = G_define_flag();
+ flag.quiet->key = 'q';
+ flag.quiet->description = _("Run quietly");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get name of input file */
+ oldname = parm.input->answer;
+
+ /* test input files existance */
+ if ((oldmapset = G_find_cell2(oldname, "")) == NULL) {
+ G_fatal_error(_("%s: <%s> raster file not found\n"), G_program_name(),
+ oldname);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+
+ /* get verbose */
+ verbose = !flag.quiet->answer;
+
+ /* get number of cell-neighbors */
+ neighb_count = flag.adjacent->answer ? 8 : 4;
+
+ /* check if the id file name is correct */
+ idname = parm.id_rast->answer;
+ if (idname) {
+ if (G_legal_filename(idname) < 0) {
+ G_fatal_error(_("%s: <%s> illegal file name\n"), G_program_name(),
+ idname);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ idmapset = G_mapset();
+ }
+
+ /* check if the patch file name is correct */
+ patchname = parm.patch_rast->answer;
+ if (patchname) {
+ if (G_legal_filename(patchname) < 0) {
+ G_fatal_error(_("%s: <%s> illegal file name\n"), G_program_name(),
+ patchname);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ patchmapset = G_mapset();
+ }
+
+ /* get size */
+ sx = G_window_cols();
+ sy = G_window_rows();
+
+ /* scan all statmethod answers */
+ stat_count = 0;
+ while (parm.stats->answers[stat_count] != NULL) {
+ /* get actual method */
+ for (method = 0; (str = statmethods[method].name); method++)
+ if (strcmp(str, parm.stats->answers[stat_count]) == 0)
+ break;
+ if (!str) {
+ G_warning(_("<%s=%s> unknown %s"),
+ parm.stats->key, parm.stats->answers[stat_count],
+ parm.stats->key);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+
+ stats[stat_count] = method;
+
+ stat_count++;
+ }
+
+ /* allocate map buffers */
+ cells = (Coords *) G_malloc(sx * sy * sizeof(Coords));
+ fragments = (Coords **) G_malloc(sx * sy * sizeof(Coords *));
+ fragments[0] = cells;
+ map = (DCELL *) G_malloc(sx * sy * sizeof(DCELL));
+ id_map = (int *)G_malloc(sx * sy * sizeof(int));
+ d_res = G_allocate_d_raster_buf();
+ result = G_allocate_c_raster_buf();
+
+ G_set_c_null_value(id_map, sx * sy);
+
+ /* open map */
+ if ((in_fd = G_open_cell_old(oldname, oldmapset)) < 0) {
+ G_fatal_error(_("can't open cell file <%s> in mapset %s\n"), oldname,
+ oldmapset);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+
+ /* read map */
+ if (verbose)
+ G_message("Reading map file:\n");
+ for (row = 0; row < sy; row++) {
+ G_get_d_raster_row(in_fd, map + row * sx, row);
+
+ if (verbose)
+ G_percent(row, sy, 2);
+ }
+ if (verbose)
+ G_percent(1, 1, 2);
+
+ /* close map */
+ G_close_cell(in_fd);
+
+ /* find fragments */
+ writeFragments(map, sy, sx, neighb_count);
+
+ if (verbose)
+ G_message("Writing output ... ");
+
+ /* open ASCII-file or use stdout */
+ if (parm.values->answer) {
+ if (strcmp(parm.values->answer, "-") != 0) {
+ if (!(out_fp = fopen(parm.values->answer, "w"))) {
+ G_fatal_error(_("Error creating file <%s>."),
+ parm.values->answer);
+ }
+ }
+ else {
+ out_fp = stdout;
+ }
+
+ /* write values */
+ for (i = 0; i < fragcount; i++) {
+ double area = fragments[i + 1] - fragments[i];
+
+ values = (DCELL *) G_malloc(area * sizeof(DCELL));
+ for (p = fragments[i], j = 0; p < fragments[i + 1]; p++, j++) {
+ values[j] = p->value;
+ }
+
+ /* write patch index */
+ fprintf(out_fp, "%d %lf", i, area / (sx * sy));
+
+ /* write patch value */
+ for (method = 0; method < stat_count; method++) {
+ DCELL outval;
+
+ perform_method = statmethods[stats[method]].method;
+ outval = perform_method(values, area);
+
+ fprintf(out_fp, " %lf", outval);
+ }
+ fprintf(out_fp, "\n");
+
+ G_free(values);
+ }
+
+ /* close file */
+ if (out_fp != stdout) {
+ fclose(out_fp);
+ }
+ else {
+ fprintf(out_fp, "\n");
+ }
+ }
+
+ /* open ASCII-file or use stdout */
+ if (parm.output->answer && strcmp(parm.output->answer, "-") != 0) {
+ if (!(out_fp = fopen(parm.output->answer, "w"))) {
+ G_fatal_error(_("Error creating file <%s>."),
+ parm.output->answer);
+ }
+ }
+ else {
+ out_fp = stdout;
+ }
+
+ /* write method names */
+ for (method = 0; method < stat_count; method++) {
+ fprintf(out_fp, "%s ", statmethods[stats[method]].suffix);
+ }
+ fprintf(out_fp, "landcover number\n");
+
+ /* allocate an array for cell values */
+ count = fragments[fragcount] - fragments[0];
+ values = (DCELL *) G_malloc(count * sizeof(DCELL));
+ for (p = fragments[0], j = 0; p < fragments[fragcount]; p++, j++) {
+ values[j] = p->value;
+ }
+
+ /* write values */
+ for (method = 0; method < stat_count; method++) {
+ perform_method = statmethods[stats[method]].method;
+ val = perform_method(values, count);
+ fprintf(out_fp, "%f ", val);
+ }
+
+ /* write landcover and number of patches */
+ area = fragments[fragcount] - fragments[0];
+ fprintf(out_fp, "%lf %d\n", (DCELL) area / (DCELL) (sx * sy), fragcount);
+
+ G_free(values);
+
+ /* close file */
+ if (parm.output->answer && strcmp(parm.output->answer, "-") != 0) {
+ fclose(out_fp);
+ }
+
+ /* write id raster */
+ if (idname) {
+ /* open new cellfile */
+ out_fd = G_open_raster_new(idname, CELL_TYPE);
+ if (out_fd < 0) {
+ G_fatal_error(_("can't create new cell file <%s> in mapset %s\n"),
+ idname, idmapset);
+ exit(EXIT_FAILURE);
+ }
+
+ /* write the output file */
+ for (row = 0; row < sy; row++) {
+ G_put_c_raster_row(out_fd, id_map + row * sx);
+
+ G_percent(row + 1, 2 * sy, 1);
+ }
+
+ /* close output */
+ G_close_cell(out_fd);
+ }
+
+ /* write patch raster */
+ if (patchname) {
+ /* open new cellfile */
+ out_fd = G_open_raster_new(patchname, CELL_TYPE);
+ if (out_fd < 0) {
+ G_fatal_error(_("can't create new cell file <%s> in mapset %s\n"),
+ patchname, patchmapset);
+ exit(EXIT_FAILURE);
+ }
+
+ /* write the output file */
+ for (row = 0; row < sy; row++) {
+ G_set_c_null_value(result, sx);
+ for (col = 0; col < sx; col++) {
+ if (!G_is_c_null_value(id_map + row * sx + col)) {
+ result[col] = 1;
+ }
+ }
+ G_put_c_raster_row(out_fd, result);
+
+ G_percent(sy + row + 1, 2 * sy, 1);
+ }
+
+ /* close output */
+ G_close_cell(out_fd);
+ }
+
+ /* free allocated resources */
+ G_free(map);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/raster/r.pi/r.pi.export/stat_method.c
===================================================================
--- grass-addons/raster/r.pi/r.pi.export/stat_method.c (rev 0)
+++ grass-addons/raster/r.pi/r.pi.export/stat_method.c 2010-08-04 10:53:56 UTC (rev 42991)
@@ -0,0 +1,125 @@
+#include "local_proto.h"
+
+DCELL average(DCELL * vals, int count)
+{
+ int i;
+
+ DCELL res = 0;
+
+ if (count <= 0)
+ return 0;
+
+ for (i = 0; i < count; i++)
+ res += vals[i];
+
+ return res / count;
+}
+
+DCELL variance(DCELL * vals, int count)
+{
+ int i;
+
+ DCELL mean;
+
+ DCELL s = 0;
+
+ DCELL ss = 0;
+
+ if (count <= 0)
+ return 0;
+
+ for (i = 0; i < count; i++) {
+ DCELL val = vals[i];
+
+ s += val;
+ ss += val * val;
+ }
+
+ mean = s / (DCELL) count;
+ return ss / count - mean * mean;
+}
+
+DCELL std_deviat(DCELL * vals, int count)
+{
+ if (count <= 0)
+ return 0;
+
+ return sqrt(variance(vals, count));
+}
+
+DCELL median(DCELL * vals, int count)
+{
+ int k = (count - 1) / 2;
+
+ int l = 0;
+
+ int h = count - 1;
+
+ DCELL pivot, tmp;
+
+ int i, j, z;
+
+ if (count <= 0)
+ return 0;
+
+ while (l < h) {
+ pivot = vals[k];
+ i = l;
+ j = h;
+
+ do {
+ while (vals[i] < pivot)
+ i++;
+ while (vals[j] > pivot)
+ j--;
+ if (i <= j) {
+ tmp = vals[i];
+ vals[i] = vals[j];
+ vals[j] = tmp;
+ i++;
+ j--;
+ }
+ } while (i <= j);
+
+ if (j < k)
+ l = i;
+ if (i > k)
+ h = j;
+ }
+
+ return vals[k];
+}
+
+DCELL min(DCELL * vals, int count)
+{
+ int i;
+
+ DCELL res = 0;
+
+ if (count <= 0)
+ return 0;
+
+ res = vals[0];
+ for (i = 0; i < count; i++)
+ if (vals[i] < res)
+ res = vals[i];
+
+ return res;
+}
+
+DCELL max(DCELL * vals, int count)
+{
+ int i;
+
+ DCELL res = 0;
+
+ if (count <= 0)
+ return 0;
+
+ res = vals[0];
+ for (i = 0; i < count; i++)
+ if (vals[i] > res)
+ res = vals[i];
+
+ return res;
+}
More information about the grass-commit
mailing list