r.out.gdal: improve range checks for float and double

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@72961 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.out.gdal/export_band.c

@@ -229,6 +229,7 @@ int export_band(GDALDatasetH hMEMDS, int band,
     int rows = cellhead->rows;
     int ret = 0;
     char value[200];
+    char *units;
 
     /* Open GRASS raster */
     fd = Rast_open_old(name, mapset);
@@ -244,6 +245,10 @@ int export_band(GDALDatasetH hMEMDS, int band,
     if (!no_metadata)
 	GDALSetDescription(hBand, name);
 
+    units = Rast_read_units(name, mapset);
+    if (units)
+	GDALSetRasterUnitType(hBand, units);
+
     /* Get min/max values. */
     if (Rast_read_fp_range(name, mapset, &sRange) == -1) {
 	bHaveMinMax = FALSE;
@@ -536,9 +541,11 @@ int exact_range_check(double min, double max, GDALDataType datatype,
 
     case GDT_Float32:
     case GDT_CFloat32:
-	if (min < TYPE_FLOAT32_MIN || max > TYPE_FLOAT32_MAX) {
+	/* support export of inf / -inf ? */
+	if ((float)min !=TYPE_FLOAT32_MIN && (float)max != TYPE_FLOAT32_MAX &&
+	    (min < TYPE_FLOAT32_MIN || max > TYPE_FLOAT32_MAX)) {
 	    G_warning(_("Selected GDAL datatype does not cover data range."));
-	    G_warning(_("GDAL datatype: %s, range: %g - %g"),
+	    G_warning(_("GDAL datatype: %s, range: %.7g - %.7g"),
 		      GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
 		      TYPE_FLOAT32_MAX);
 	    G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
@@ -550,9 +557,10 @@ int exact_range_check(double min, double max, GDALDataType datatype,
     case GDT_Float64:
     case GDT_CFloat64:
 	/* not possible because DCELL is FLOAT64, not 128bit floating point, but anyway... */
+	/* support export of inf / -inf ? */
 	if (min < TYPE_FLOAT64_MIN || max > TYPE_FLOAT64_MAX) {
 	    G_warning(_("Selected GDAL datatype does not cover data range."));
-	    G_warning(_("GDAL datatype: %s, range: %g - %g"),
+	    G_warning(_("GDAL datatype: %s, range: %.16g - %.16g"),
 		      GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
 		      TYPE_FLOAT64_MAX);
 	    G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);

raster/r.out.gdal/local_proto.h

@@ -2,6 +2,8 @@
 #define __LOCAL_PROTO_H__
 
 #include <gdal.h>
+#include <math.h>
+/* is it safe to #include <float.h> ? */
 
 /* range limits */
 /*
@@ -36,8 +38,9 @@
 #define TYPE_FLOAT32_MIN	(-MAXFLOAT)
 #define TYPE_FLOAT32_MAX	MAXFLOAT
 #else
-#define TYPE_FLOAT32_MIN	-3.4E38f
-#define TYPE_FLOAT32_MAX	3.4E38f
+/* formula from float.h */
+#define TYPE_FLOAT32_MIN	(float)(-((1 + (1 - pow(2, -23))) * pow(2, 127)))
+#define TYPE_FLOAT32_MAX	(float)((1 + (1 - pow(2, -23))) * pow(2, 127))
 #endif
 
 #ifdef DBL_MAX
@@ -47,8 +50,9 @@
 #define TYPE_FLOAT64_MIN	(-MAXDOUBLE)
 #define TYPE_FLOAT64_MAX	MAXDOUBLE
 #else
-#define TYPE_FLOAT64_MIN	-1.79E308
-#define TYPE_FLOAT64_MAX	1.79E308
+/* formula from float.h */
+#define TYPE_FLOAT64_MIN	(-((1 + (1 - pow(2, -52))) * pow(2, 1023)))
+#define TYPE_FLOAT64_MAX	((1 + (1 - pow(2, -52))) * pow(2, 1023))
 #endif
 
 #define GRASS_MAX_COLORS TYPE_UINT16_MAX    /* ok? */

raster/r.out.gdal/main.c

@@ -775,7 +775,7 @@ int range_check(double min, double max, GDALDataType datatype)
 	if ((!CPLIsInf(max) && max < TYPE_FLOAT32_MIN) ||
 	    (!CPLIsInf(min) && min > TYPE_FLOAT32_MAX)) {
 	    G_warning(_("Selected GDAL datatype does not cover data range."));
-	    G_warning(_("GDAL datatype: %s, range: %g - %g"),
+	    G_warning(_("GDAL datatype: %s, range: %.7g - %.7g"),
 		      GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
 		      TYPE_FLOAT32_MAX);
 	    G_warning(_("Range to be exported: %g - %g"), min, max);