r.external speed-up (#1305)

 * faster linking with r.external
 * super-fast mode for r.external
 * warning for fast mode
 * update r.info if range is missing
 * fix r.support to re-compute range for fp maps

Co-authored-by: Markus Neteler <neteler@gmail.com>

lib/raster/range.c

@@ -120,6 +120,12 @@ int Rast_read_fp_range(const char *name, const char *mapset,
 	Rast_update_fp_range(dcell2, drange);
 	close(fd);
     }
+    else {
+	/* "f_range" file does not exist */
+	G_warning(_("Missing fp range file for <%s> (run r.support -s)"),
+		  G_fully_qualified_name(name, mapset));
+	return -1;
+    }
 
     return 1;
 }
@@ -237,6 +243,12 @@ int Rast_read_range(const char *name, const char *mapset, struct Range *range)
 	}
 	fclose(fd);
     }
+    else {
+	/* "range" file does not exist */
+	G_warning(_("Missing range file for <%s> (run r.support -s)"),
+		  G_fully_qualified_name(name, mapset));
+	return -1;
+    }
 
     return 1;
 }

raster/r.external/link.c

@@ -47,6 +47,9 @@ void query_band(GDALRasterBandH hBand, const char *output,
 	break;
     }
 
+    if (info->have_minmax)
+	GDALComputeRasterMinMax(hBand, 0, info->minmax);
+
     Rast_init_colors(&info->colors);
 
     if (GDALGetRasterColorTable(hBand) != NULL) {
@@ -184,7 +187,6 @@ void create_map(const char *input, int band, const char *output,
     struct History history;
     struct Categories cats;
     char buf[1024];
-    int outfd;
 
     Rast_put_cellhd(output, cellhd);
 
@@ -211,39 +213,24 @@ void create_map(const char *input, int band, const char *output,
     if (title)
 	Rast_put_cell_title(output, title);
 
-    /* get stats for this raster band */
-    G_remove_misc("cell_misc", "stats", output);
-
-    outfd = Rast_open_old(output, G_mapset());
-    if (info->data_type == CELL_TYPE) {
-	int r;
-	struct Range range;
-	CELL *rbuf = Rast_allocate_buf(CELL_TYPE);
+    if (info->have_minmax) {
+	if (info->data_type == CELL_TYPE) {
+	    struct Range range;
 
-	G_remove_misc("cell_misc", "range", output);
-	Rast_init_range(&range);
-	for (r = 0; r < cellhd->rows; r++) {
-	    Rast_get_row(outfd, rbuf, r, CELL_TYPE);
-	    Rast_row_update_range(rbuf, cellhd->cols, &range);
+	    Rast_init_range(&range);
+	    Rast_update_range((CELL)info->minmax[0], &range);
+	    Rast_update_range((CELL)info->minmax[1], &range);
+	    Rast_write_range(output, &range);
 	}
-	Rast_write_range(output, &range);
-	G_free(rbuf);
-    }
-    else {
-	int r;
-	struct FPRange fprange;
-	void *rbuf = Rast_allocate_buf(info->data_type);
-	
-	G_remove_misc("cell_misc", "f_range", output);
-	Rast_init_fp_range(&fprange);
-	for (r = 0; r < cellhd->rows; r++) {
-	    Rast_get_row(outfd, rbuf, r, info->data_type);
-	    Rast_row_update_fp_range(rbuf, cellhd->cols, &fprange, info->data_type);
+	else {
+	    struct FPRange fprange;
+
+	    Rast_init_fp_range(&fprange);
+	    Rast_update_fp_range(info->minmax[0], &fprange);
+	    Rast_update_fp_range(info->minmax[1], &fprange);
+	    Rast_write_fp_range(output, &fprange);
 	}
-	Rast_write_fp_range(output, &fprange);
-	G_free(rbuf);
     }
-    Rast_unopen(outfd);
 
     G_message(_("Link to raster map <%s> created."), output);
 }

raster/r.external/main.c

@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
 	struct Option *input, *source, *output, *band, *title;
     } parm;
     struct {
-	struct Flag *o, *j, *f, *e, *h, *v, *t, *a;
+	struct Flag *o, *j, *f, *e, *h, *v, *t, *a, *r;
     } flag;
     int min_band, max_band, band;
     struct band_info info;
@@ -135,6 +135,11 @@ int main(int argc, char *argv[])
     flag.t->guisection = _("Print");
     flag.t->suppress_required = YES;
 
+    flag.r = G_define_flag();
+    flag.r->key = 'r';
+    flag.r->label = _("Create fast link without data range");
+    flag.r->description = _("WARNING: some modules do not work correctly without known data range");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -219,6 +224,7 @@ int main(int argc, char *argv[])
 	I_init_group_ref(&reference);
     }
 
+    info.have_minmax = !flag.r->answer;
     for (band = min_band; band <= max_band; band++) {
 	char *output2, *title2 = NULL;
 

raster/r.external/proto.h

@@ -11,6 +11,8 @@ struct band_info
     GDALDataType gdal_type;
     int has_null;
     double null_val;
+    DCELL minmax[2];
+    int have_minmax;
     struct Colors colors;
 };
 

raster/r.info/main.c

@@ -47,7 +47,6 @@ int main(int argc, char **argv)
     char *units, *vdatum;
     int i;
     CELL mincat = 0, maxcat = 0, cat;
-    double zmin, zmax;		/* min and max data values */
     FILE *out;
     struct Range crange;
     struct FPRange range;
@@ -130,10 +129,6 @@ int main(int argc, char **argv)
 	    second_time_ok = 1;
     }
 
-    if (Rast_read_fp_range(name, "", &range) < 0)
-	G_fatal_error(_("Unable to read range file"));
-    Rast_get_fp_range_min_max(&range, &zmin, &zmax);
-
     out = stdout;
 
     if (eflag->answer || (!gflag->answer && !rflag->answer && !sflag->answer && !hflag->answer)) {
@@ -225,27 +220,54 @@ int main(int argc, char **argv)
 			 tmp1, tmp2, tmp3);
 
 	    if (data_type == CELL_TYPE) {
-		if (2 == Rast_read_range(name, "", &crange))
+		int ret;
+		CELL min, max;
+
+		/* print range only if available */
+		ret = Rast_read_range(name, "", &crange);
+		if (ret == 2)
 		    compose_line(out,
 				 "  Range of data:    min = NULL  max = NULL");
-		else
-		    compose_line(out,
-				 "  Range of data:    min = %i  max = %i",
-				 (CELL) zmin, (CELL) zmax);
+		else if (ret > 0) {
+		    Rast_get_range_min_max(&crange, &min, &max);
+
+		    if (Rast_is_c_null_value(&min)) {
+			compose_line(out,
+				     "  Range of data:    min = NULL  max = NULL");
+		    }
+		    else {
+			compose_line(out,
+				     "  Range of data:    min = %i  max = %i",
+				     min, max);
+		    }
+		}
 	    }
 	    else {
-		if (Rast_is_d_null_value(&zmin)) {
+		int ret;
+		DCELL min, max;
+
+		/* print range only if available */
+		ret = Rast_read_fp_range(name, "", &range);
+		if (ret == 2) {
 		    compose_line(out,
 				 "  Range of data:    min = NULL  max = NULL");
 		}
-		else {
-		    if (data_type == FCELL_TYPE) {
-			compose_line(out, "  Range of data:    min = %.7g  max = %.7g",
-				     zmin, zmax);
+		else if (ret > 0) {
+		    Rast_get_fp_range_min_max(&range, &min, &max);
+
+		    if (Rast_is_d_null_value(&min)) {
+			compose_line(out,
+				     "  Range of data:    min = NULL  max = NULL");
 		    }
 		    else {
-			compose_line(out, "  Range of data:    min = %.15g  max = %.15g",
-				     zmin, zmax);
+			if (data_type == FCELL_TYPE) {
+			    compose_line(out, "  Range of data:    min = %.7g  max = %.7g",
+					 min, max);
+			}
+			else {
+			    compose_line(out, "  Range of data:    min = %.15g  max = %.15g",
+					 min, max);
+			}
 		    }
 		}
 	    }
@@ -325,7 +347,83 @@ int main(int argc, char **argv)
 
 	fprintf(out, "\n");
     }
-    else {	/* g,r,s, e, or h flags */
+    else {	/* g, r, s, e, or h flags */
+	int need_range, have_range, need_stats, have_stats;
+	
+	need_range = rflag->answer;
+	need_stats = sflag->answer;
+	if (need_stats)
+	    need_range = 1;
+
+	have_range = have_stats = 0;
+	if (need_range) {
+	    if (data_type == CELL_TYPE) {
+		if (Rast_read_range(name, "", &crange) > 0)
+		    have_range = 1;
+	    }
+	    else {
+		if (Rast_read_fp_range(name, "", &range) > 0)
+		    have_range = 1;
+	    }
+	}
+	if (need_stats) {
+	    if (Rast_read_rstats(name, mapset, &rstats) > 0)
+		have_stats = 1;
+	}
+	
+	if ((need_stats && !have_stats) || (need_range && !have_range)) {
+	    DCELL *dbuf, val, min, max;
+	    int fd, r, c;
+	    int first = 1;
+
+	    Rast_set_input_window(&cellhd);
+	    dbuf = Rast_allocate_d_input_buf();
+	    fd = Rast_open_old(name, mapset);
+	    min = max = 0;
+
+	    for (r = 0; r < cellhd.rows; r++) {
+		Rast_get_d_row_nomask(fd, dbuf, r);
+		for (c = 0; c < cellhd.cols; c++) {
+		    val = dbuf[c];
+		    if (Rast_is_d_null_value(&val))
+			continue;
+		    if (first) {
+			rstats.sum = val;
+			rstats.sumsq = (DCELL) val * val;
+			rstats.count = 1;
+			min = max = val;
+
+			first = 0;
+		    }
+		    else {
+			rstats.sum += val;
+			rstats.sumsq += (DCELL) val * val;
+			rstats.count += 1;
+			if (min > val)
+			    min = val;
+			if (max < val)
+			    max = val;
+		    }
+		}
+	    }
+	    Rast_close(fd);
+	    G_free(dbuf);
+
+	    if (data_type == CELL_TYPE) {
+		Rast_init_range(&crange);
+		if (rstats.count > 0) {
+		    Rast_update_range((CELL) min, &crange);
+		    Rast_update_range((CELL) max, &crange);
+		}
+	    }
+	    else {
+		Rast_init_fp_range(&range);
+		if (rstats.count > 0) {
+		    Rast_update_fp_range(min, &range);
+		    Rast_update_fp_range(max, &range);
+		}
+	    }
+	}
 
 	if (gflag->answer) {
 	    G_format_northing(cellhd.north, tmp1, -1);
@@ -360,30 +458,36 @@ int main(int argc, char **argv)
                     cats_ok ? tmp1 : "??");
 	}
 
-	if (rflag->answer) {
+	if (rflag->answer || sflag->answer) {
 	    if (data_type == CELL_TYPE) {
-		if (2 == Rast_read_range(name, "", &crange)) {
+		CELL min, max;
+		
+		Rast_get_range_min_max(&crange, &min, &max);
+		if (Rast_is_c_null_value(&min)) {
 		    fprintf(out, "min=NULL\n");
 		    fprintf(out, "max=NULL\n");
 		}
 		else {
-		    fprintf(out, "min=%i\n", (CELL) zmin);
-		    fprintf(out, "max=%i\n", (CELL) zmax);
+		    fprintf(out, "min=%i\n", min);
+		    fprintf(out, "max=%i\n", max);
 		}
 	    }
 	    else {
-		if (Rast_is_d_null_value(&zmin)) {
+		DCELL min, max;
+
+		Rast_get_fp_range_min_max(&range, &min, &max);
+		if (Rast_is_d_null_value(&min)) {
 		    fprintf(out, "min=NULL\n");
 		    fprintf(out, "max=NULL\n");
 		}
 		else {
 		    if (data_type == FCELL_TYPE) {
-			fprintf(out, "min=%.7g\n", zmin);
-			fprintf(out, "max=%.7g\n", zmax);
+			fprintf(out, "min=%.7g\n", min);
+			fprintf(out, "max=%.7g\n", max);
 		    }
 		    else {
-			fprintf(out, "min=%.15g\n", zmin);
-			fprintf(out, "max=%.15g\n", zmax);
+			fprintf(out, "min=%.15g\n", min);
+			fprintf(out, "max=%.15g\n", max);
 		    }
 		}
 	    }
@@ -391,39 +495,6 @@ int main(int argc, char **argv)
 
 	if (sflag->answer) {
 
-	    if (Rast_read_rstats(name, mapset, &rstats) < 0) {
-		DCELL *dbuf, val;
-		int fd, r, c;
-		int first = 1;
-
-		Rast_set_input_window(&cellhd);
-		dbuf = Rast_allocate_d_input_buf();
-		fd = Rast_open_old(name, mapset);
-
-		for (r = 0; r < cellhd.rows; r++) {
-		    Rast_get_d_row(fd, dbuf, r);
-		    for (c = 0; c < cellhd.cols; c++) {
-			val = dbuf[c];
-			if (Rast_is_d_null_value(&val))
-			    continue;
-			if (first) {
-			    rstats.sum = val;
-			    rstats.sumsq = (DCELL) val * val;
-			    rstats.count = 1;
-
-			    first = 0;
-			}
-			else {
-			    rstats.sum += val;
-			    rstats.sumsq += (DCELL) val * val;
-			    rstats.count += 1;
-			}
-		    }
-		}
-		Rast_close(fd);
-		G_free(dbuf);
-	    }
-
 	    if (!gflag->answer) {
 		/* always report total number of cells */
 		fprintf(out, "cells=%jd\n",
@@ -436,28 +507,33 @@ int main(int argc, char **argv)
 		mean = (double)(rstats.sum / rstats.count);
 		sd = sqrt(rstats.sumsq / rstats.count - (mean * mean));
 
-		fprintf(out, "n=%jd\n", rstats.count);
 
-		if (!rflag->answer) {
-		    fprintf(out, "min=%.15g\n", zmin);
-		    fprintf(out, "max=%.15g\n", zmax);
-		}
-		if (zmin == zmax) {
-		    fprintf(out, "mean=%.15g\n", zmin);
-		    fprintf(out, "stddev=0\n");
+		if (data_type == CELL_TYPE) {
+		    CELL min, max;
+		
+		    Rast_get_range_min_max(&crange, &min, &max);
+		    if (min == max) {
+			mean = min;
+			sd = 0;
+		    }
 		}
 		else {
-		    fprintf(out, "mean=%.15g\n", mean);
-		    fprintf(out, "stddev=%.15g\n", sd);
+		    DCELL min, max;
+
+		    Rast_get_fp_range_min_max(&range, &min, &max);
+		    if (min == max) {
+			mean = min;
+			sd = 0;
+		    }
 		}
+
+		fprintf(out, "n=%jd\n", rstats.count);
+		fprintf(out, "mean=%.15g\n", mean);
+		fprintf(out, "stddev=%.15g\n", sd);
 		fprintf(out, "sum=%.15g\n", rstats.sum);
 	    }
 	    else {
 		fprintf(out, "n=0\n");
-		if (!rflag->answer) {
-		    fprintf(out, "min=NULL\n");
-		    fprintf(out, "max=NULL\n");
-		}
 		fprintf(out, "mean=NULL\n");
 		fprintf(out, "stddev=NULL\n");
 		fprintf(out, "sum=NULL\n");

raster/r.support/check.c

@@ -13,7 +13,6 @@
 int check_stats(const char *name)
 {
     RASTER_MAP_TYPE data_type;
-    struct Histogram histogram;
     struct Categories cats;
     struct Range range;
     struct FPRange fprange;
@@ -27,34 +26,15 @@ int check_stats(const char *name)
     if (do_histogram(name) < 0)
 	return 0;
 
-    if (Rast_read_histogram(name, "", &histogram) <= 0)
-	return 0;
-
-    /* Init histogram range */
-    if (data_type == CELL_TYPE)
-	Rast_init_range(&range);
-    else
-	Rast_init_fp_range(&fprange);
-
-    /* Update histogram range */
-    i = Rast_get_histogram_num(&histogram);
-    while (i >= 0) {
-	if (data_type == CELL_TYPE)
-	    Rast_update_range(Rast_get_histogram_cat(i--, &histogram), &range);
-	else
-	    Rast_update_fp_range((DCELL) Rast_get_histogram_cat(i--, &histogram),
-			      &fprange);
-    }
-
-    /* Write histogram range */
+    /* get range */
     if (data_type == CELL_TYPE)
-	Rast_write_range(name, &range);
+	Rast_read_range(name, "", &range);
     else
-	Rast_write_fp_range(name, &fprange);
+	Rast_read_fp_range(name, "", &fprange);
+    max = (data_type == CELL_TYPE ? range.max : fprange.max);
 
-    /* Get category status and max */
+    /* Get category status */
     cats_ok = (Rast_read_cats(name, "", &cats) >= 0);
-    max = (data_type == CELL_TYPE ? range.max : fprange.max);
 
     /* Further category checks */
     if (!cats_ok)
@@ -69,10 +49,8 @@ int check_stats(const char *name)
 	G_message(_("   Updating the number of categories for "
 		    "[%s]\n\n"), name);
 	Rast_write_cats(name, &cats);
+	Rast_free_cats(&cats);
     }
 
-    Rast_free_histogram(&histogram);
-    Rast_free_cats(&cats);
-
     return 0;
 }

raster/r.support/histo.c

@@ -11,6 +11,7 @@
  */
 int do_histogram(const char *name)
 {
+    RASTER_MAP_TYPE data_type;
     CELL *cell;
     struct Cell_head cellhd;
     struct Cell_stats statf;
@@ -19,6 +20,7 @@ int do_histogram(const char *name)
     int fd;
 
     Rast_get_cellhd(name, "", &cellhd);
+    data_type = Rast_map_type(name, "");
 
     Rast_set_window(&cellhd);
     fd = Rast_open_old(name, "");
@@ -28,20 +30,47 @@ int do_histogram(const char *name)
     cell = Rast_allocate_c_buf();
 
     Rast_init_cell_stats(&statf);
-    for (row = 0; row < nrows; row++) {
-	Rast_get_c_row_nomask(fd, cell, row);
-	Rast_update_cell_stats(cell, ncols, &statf);
+    if (data_type == CELL_TYPE) {
+	struct Range range;
+
+	Rast_init_range(&range);
+	for (row = 0; row < nrows; row++) {
+	    Rast_get_c_row_nomask(fd, cell, row);
+	    Rast_update_cell_stats(cell, ncols, &statf);
+	    Rast_row_update_range(cell, ncols, &range);
+	}
+	Rast_write_range(name, &range);
+    }
+    else {
+	DCELL *dcell;
+	struct FPRange fprange;
+	int col;
+
+	/* histogram for FCELL / DCELL type ? */
+	
+	Rast_init_fp_range(&fprange);
+	dcell = Rast_allocate_d_buf();
+	for (row = 0; row < nrows; row++) {
+	    Rast_get_d_row_nomask(fd, dcell, row);
+	    /* Rast_update_cell_stats wants CELL values */
+	    for (col = 0; col < ncols; col++) {
+		if (Rast_is_d_null_value(&dcell[col]))
+		    Rast_set_c_null_value(&cell[col], 1);
+		else
+		    cell[col] = (CELL) dcell[col];
+	    }
+	    Rast_update_cell_stats(cell, ncols, &statf);
+	    Rast_row_update_fp_range(dcell, ncols, &fprange, data_type);
+	}
+	Rast_write_fp_range(name, &fprange);
+	G_free(dcell);
     }
 
-    if (row == nrows)
-	Rast_write_histogram_cs(name, &statf);
+    Rast_write_histogram_cs(name, &statf);
 
     Rast_free_cell_stats(&statf);
     Rast_close(fd);
     G_free(cell);
 
-    if (row < nrows)
-	return -1;
-
     return 0;
 }