[GRASS-SVN] r36002 - grass-addons/raster/r.univar2.zonal

svn_grass at osgeo.org svn_grass at osgeo.org
Sat Feb 21 02:57:19 EST 2009


Author: mmetz
Date: 2009-02-21 02:57:19 -0500 (Sat, 21 Feb 2009)
New Revision: 36002

Modified:
   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/stats.c
Log:
option for table format ouput instead of standard r.univar output

Modified: grass-addons/raster/r.univar2.zonal/globals.h
===================================================================
--- grass-addons/raster/r.univar2.zonal/globals.h	2009-02-20 22:34:02 UTC (rev 36001)
+++ grass-addons/raster/r.univar2.zonal/globals.h	2009-02-21 07:57:19 UTC (rev 36002)
@@ -45,13 +45,14 @@
 {
     CELL min, max, n_zones;
     struct Categories cats;
+    char *sep;
 } zone_type;
 
 /* 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;
+    struct Flag *shell_style, *extended, *table;
 } param_type;
 
 #ifdef MAIN
@@ -67,6 +68,7 @@
 void heapsort_float(float *data, int n);
 void heapsort_int(int *data, int n);
 int print_stats(univar_stat * stats);
+int print_stats2(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);
 

Modified: grass-addons/raster/r.univar2.zonal/r.univar.html
===================================================================
--- grass-addons/raster/r.univar2.zonal/r.univar.html	2009-02-20 22:34:02 UTC (rev 36001)
+++ grass-addons/raster/r.univar2.zonal/r.univar.html	2009-02-21 07:57:19 UTC (rev 36002)
@@ -1,13 +1,18 @@
 <h2>DESCRIPTION</h2>
 
-<em>r.univar</em> calculates the univariate statistics of a raster map.
+<em>r.univar.zonal</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.
+coefficient of variation. Statistics are calculated separately for every
+category/zone found in the <b>zones</b> input map.
 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 ";" as field separator. The table can immediately be converted to a
+proper 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.
 
 <h2>NOTES</h2>
 

Modified: grass-addons/raster/r.univar2.zonal/r.univar_main.c
===================================================================
--- grass-addons/raster/r.univar2.zonal/r.univar_main.c	2009-02-20 22:34:02 UTC (rev 36001)
+++ grass-addons/raster/r.univar2.zonal/r.univar_main.c	2009-02-21 07:57:19 UTC (rev 36002)
@@ -36,7 +36,7 @@
 	_("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->required = YES;
     param.output_file->description =
 	_("Name for output file (if omitted or \"-\" output to stdout)");
 
@@ -59,6 +59,10 @@
     param.extended->key = 'e';
     param.extended->description = _("Calculate extended statistics");
 
+    param.table = G_define_flag();
+    param.table->key = 't';
+    param.table->description = _("Table output format instead of r.univar like output format");
+
     return;
 }
 
@@ -104,6 +108,9 @@
 	}
     }
 
+    /* TODO: make it an option */
+    zone_info.sep = ";";
+
     G_get_window(&region);
     rows = region.rows;
     cols = region.cols;
@@ -161,12 +168,16 @@
 
     process_raster(stats, fd, fdz, &region);
 
-    if (!(param.shell_style->answer))
-	G_percent(rows, rows, 2);	/* finish it off */
+    /* closing raster maps */
+    G_close_cell(fd);
+    G_close_cell(fdz);
 
     /* create the output */
-    print_stats(stats);
-
+    if (param.table->answer)
+	print_stats2(stats);
+    else
+	print_stats(stats);
+	
     /* release memory */
     free_univar_stat_struct(stats);
 
@@ -279,19 +290,22 @@
 
 	for (col = 0; col < cols; col++) {
 
-	    if (G_is_null_value(ptr, map_type)) {
+	    if (G_is_c_null_value(zptr)) {
 		ptr = G_incr_void_ptr(ptr, value_sz);
 		zptr++;
 		continue;
 	    }
-	    if (G_is_c_null_value(zptr)) {
+
+	    /* also count NULL cell in input map */
+	    zone = *zptr - zone_info.min;
+	    stats[zone].size++;
+	    
+	    if (G_is_null_value(ptr, map_type)) {
 		ptr = G_incr_void_ptr(ptr, value_sz);
 		zptr++;
 		continue;
 	    }
 
-	    zone = *zptr - zone_info.min;
-
 	    if (stats[zone].nextp) {
 		/* put the value into stats->XXXcell_array */
 		memcpy(stats[zone].nextp, ptr, value_sz);
@@ -324,7 +338,8 @@
 	    zptr++;
 	    stats[zone].n++;
 	}
-	if (!(param.shell_style->answer))
-	    G_percent(row, rows, 2);
+	G_percent(row, rows, 2);
     }
+    G_percent(rows, rows, 2);	/* finish it off */
+
 }

Modified: grass-addons/raster/r.univar2.zonal/stats.c
===================================================================
--- grass-addons/raster/r.univar2.zonal/stats.c	2009-02-20 22:34:02 UTC (rev 36001)
+++ grass-addons/raster/r.univar2.zonal/stats.c	2009-02-21 07:57:19 UTC (rev 36002)
@@ -36,7 +36,7 @@
 	    stats[i].perc = NULL;
 	stats[i].sum_abs = 0.0;
 	stats[i].n = 0;
-	stats[i].size = size;
+	stats[i].size = 0;
 	stats[i].dcell_array = NULL;
 	stats[i].fcell_array = NULL;
 	stats[i].cell_array = NULL;
@@ -46,12 +46,12 @@
 	if (param.extended->answer) {
 	    if (map_type == DCELL_TYPE)
 		stats[i].dcell_array =
-		    (DCELL *) G_calloc(stats[i].size, sizeof(DCELL));
+		    (DCELL *) G_calloc(size, sizeof(DCELL));
 	    if (map_type == FCELL_TYPE)
 		stats[i].fcell_array =
-		    (FCELL *) G_calloc(stats[i].size, sizeof(FCELL));
+		    (FCELL *) G_calloc(size, sizeof(FCELL));
 	    if (map_type == CELL_TYPE)
-		stats[i].cell_array = (CELL *) G_calloc(stats[i].size, sizeof(CELL));
+		stats[i].cell_array = (CELL *) G_calloc(size, sizeof(CELL));
 	}
     }
 
@@ -125,7 +125,7 @@
 	}
 
 	if (param.shell_style->answer) {
-	    fprintf(stdout, "n=%d\n", stats[z].n);
+	    fprintf(stdout, "non_null_cells=%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);
@@ -260,3 +260,182 @@
 
     return 1;
 }
+
+int print_stats2(univar_stat * stats)
+{
+    int z;
+    unsigned int i;
+
+    /* print column headers */
+
+    fprintf(stdout, "zone%s", zone_info.sep);
+    fprintf(stdout, "label%s", zone_info.sep);
+    fprintf(stdout, "non_null_cells%s", zone_info.sep);
+    fprintf(stdout, "null_cells%s", zone_info.sep);
+    fprintf(stdout, "min%s", zone_info.sep);
+    fprintf(stdout, "max%s", zone_info.sep);
+    fprintf(stdout, "range%s", zone_info.sep);
+    fprintf(stdout, "mean%s", zone_info.sep);
+    fprintf(stdout, "mean_of_abs%s", zone_info.sep);
+    fprintf(stdout, "stddev%s", zone_info.sep);
+    fprintf(stdout, "variance%s", zone_info.sep);
+    fprintf(stdout, "coeff_var%s", zone_info.sep);
+    fprintf(stdout, "sum%s", zone_info.sep);
+    fprintf(stdout, "sum_abs%s", zone_info.sep);
+
+    fprintf(stdout, "first_quart%s", zone_info.sep);
+    fprintf(stdout, "median%s", zone_info.sep);
+    fprintf(stdout, "third_quart%s", zone_info.sep);
+    for (i = 0; i < stats[0].n_perc; i++) {
+	fprintf(stdout, "perc_%d%s", stats[0].perc[i],
+		zone_info.sep);
+    }
+    fprintf(stdout, "\n");
+
+    /* print stats */
+
+    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;
+	int qpos_25, qpos_75, *qpos_perc;
+
+	/* stats collected for this zone? */
+	if (stats[z].n == 0)
+	    continue;
+
+	i = 0;
+
+
+	/* 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) ? */
+
+	/* zone number */
+	fprintf(stdout, "%d%s", z + zone_info.min, zone_info.sep);
+	/* zone label */
+	fprintf(stdout,"%s%s", G_get_c_raster_cat(&z + zone_info.min, &(zone_info.cats)), zone_info.sep);
+
+	/* total cells */
+	fprintf(stdout, "%d%s", stats[z].n, zone_info.sep);
+	/* null cells */
+	fprintf(stdout, "%d%s", stats[z].size - stats[z].n, zone_info.sep);
+	/* min */
+	fprintf(stdout, "%.15g%s", stats[z].min, zone_info.sep);
+	/* max */
+	fprintf(stdout, "%.15g%s", stats[z].max, zone_info.sep);
+	/* range */
+	fprintf(stdout, "%.15g%s", stats[z].max - stats[z].min, zone_info.sep);
+	/* mean */
+	fprintf(stdout, "%.15g%s", mean, zone_info.sep);
+	/* mean of abs */
+	fprintf(stdout, "%.15g%s", stats[z].sum_abs / stats[z].n, zone_info.sep);
+	/* stddev */
+	fprintf(stdout, "%.15g%s", stdev, zone_info.sep);
+	/* variance */
+	fprintf(stdout, "%.15g%s", variance, zone_info.sep);
+	/* coefficient of variance */
+	fprintf(stdout, "%.15g%s", var_coef, zone_info.sep);
+	/* sum */
+	sprintf(sum_str, "%.10f", stats[z].sum);
+	G_trim_decimal(sum_str);
+	fprintf(stdout, "%s%s", sum_str, zone_info.sep);
+	/* absolute sum */
+	sprintf(sum_str, "%.10f", stats[z].sum_abs);
+	G_trim_decimal(sum_str);
+	fprintf(stdout, "%s%s", sum_str, zone_info.sep);
+
+	/* 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;
+	    }
+
+	    /* first quartile */
+	    fprintf(stdout, "%g%s", quartile_25, zone_info.sep);
+	    /* median */
+	    fprintf(stdout, "%g%s", median, zone_info.sep);
+	    /* third quartile */
+	    fprintf(stdout, "%g%s", quartile_75, zone_info.sep);
+	    /* percentiles */
+	    for (i = 0; i < stats[z].n_perc; i++) {
+		fprintf(stdout, "%g%s", 
+			quartile_perc[i], zone_info.sep);
+	    }
+
+	    G_free((void *)quartile_perc);
+	    G_free((void *)qpos_perc);
+	}
+
+	fprintf(stdout, "\n");
+
+	/* zone z finished */
+
+    }
+
+    return 1;
+}



More information about the grass-commit mailing list