[GRASS-SVN] r37763 - grass/trunk/raster/r.out.gdal
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jun 6 12:36:20 EDT 2009
Author: mmetz
Date: 2009-06-06 12:36:19 -0400 (Sat, 06 Jun 2009)
New Revision: 37763
Modified:
grass/trunk/raster/r.out.gdal/export_band.c
grass/trunk/raster/r.out.gdal/local_proto.h
grass/trunk/raster/r.out.gdal/main.c
Log:
attempt to solve ticket #73
Modified: grass/trunk/raster/r.out.gdal/export_band.c
===================================================================
--- grass/trunk/raster/r.out.gdal/export_band.c 2009-06-06 12:31:34 UTC (rev 37762)
+++ grass/trunk/raster/r.out.gdal/export_band.c 2009-06-06 16:36:19 UTC (rev 37763)
@@ -6,7 +6,7 @@
* PURPOSE: Exports GRASS raster to GDAL suported formats;
* based on GDAL library.
*
-* COPYRIGHT: (C) 2006-2008 by the GRASS Development Team
+* COPYRIGHT: (C) 2006-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -21,12 +21,19 @@
#include "gdal.h"
#include "local_proto.h"
+int exact_range_check(double, double, GDALDataType, const char *);
-int export_band(GDALDatasetH hMEMDS, int band, const char *name,
- const char *mapset, struct Cell_head *cellhead,
- RASTER_MAP_TYPE maptype, double nodataval,
- const char *nodatakey, int suppress_main_colortable,
- int default_nodataval)
+/* actual raster band export
+ * returns 0 on success
+ * -1 on raster data read/write error
+ * -2 if given nodata value was present in data
+ * -3 if selected GDAL datatype could not hold all values
+ * */
+int export_band(GDALDatasetH hMEMDS, GDALDataType export_datatype, int band,
+ const char *name, const char *mapset,
+ struct Cell_head *cellhead, RASTER_MAP_TYPE maptype,
+ double nodataval, const char *nodatakey,
+ int suppress_main_colortable, int default_nodataval)
{
struct Colors sGrassColors;
@@ -39,6 +46,7 @@
int fd;
int cols = cellhead->cols;
int rows = cellhead->rows;
+ int ret = 0;
/* Open GRASS raster */
fd = G_open_cell_old(name, mapset);
@@ -137,6 +145,8 @@
GDALSetColorEntry(hCT, iColor, &sColor);
}
}
+
+ GDALSetRasterColorTable(hBand, hCT);
}
if (rcount > 0) {
@@ -147,6 +157,8 @@
}
/* Add the rules in reverse order */
+ /* This can cause a GDAL warning with many rules, something like
+ * Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag. */
for (i = rcount - 1; i >= 0; i--) {
DCELL val1, val2;
unsigned char r1, g1, b1, r2, g2, b2;
@@ -160,9 +172,6 @@
r2, g2, b2);
GDALSetMetadataItem(hBand, key, value, NULL);
}
-
- if (!suppress_main_colortable)
- GDALSetRasterColorTable(hBand, hCT);
}
/* Create GRASS raster buffer */
@@ -179,16 +188,24 @@
return -1;
}
- /* Copy data form GRASS raster to memory raster */
+ /* Copy data form GRASS raster to GDAL raster */
int row, col;
int n_nulls = 0, nodatavalmatch = 0;
+ dfCellMin = TYPE_FLOAT64_MAX;
+ dfCellMax = TYPE_FLOAT64_MIN;
+
+ /* Better use selected GDAL datatype instead of
+ * the best match with GRASS raster map types ? */
+
if (maptype == FCELL_TYPE) {
/* Source datatype understandable by GDAL */
GDALDataType datatype = GDT_Float32;
FCELL fnullval = (FCELL) nodataval;
+ G_debug(1, "FCELL nodata val: %f", fnullval);
+
for (row = 0; row < rows; row++) {
if (G_get_raster_row(fd, bufer, row, maptype) < 0) {
@@ -205,8 +222,14 @@
}
n_nulls++;
}
- else if (((FCELL *) bufer)[col] == fnullval) {
- nodatavalmatch = 1;
+ else {
+ if (((FCELL *) bufer)[col] == fnullval) {
+ nodatavalmatch = 1;
+ }
+ if (dfCellMin > ((FCELL *) bufer)[col])
+ dfCellMin = ((FCELL *) bufer)[col];
+ if (dfCellMax < ((FCELL *) bufer)[col])
+ dfCellMax = ((FCELL *) bufer)[col];
}
if (GDALRasterIO
@@ -223,6 +246,8 @@
GDALDataType datatype = GDT_Float64;
DCELL dnullval = (DCELL) nodataval;
+ G_debug(1, "DCELL nodata val: %f", dnullval);
+
for (row = 0; row < rows; row++) {
if (G_get_raster_row(fd, bufer, row, maptype) < 0) {
@@ -239,8 +264,14 @@
}
n_nulls++;
}
- else if (((DCELL *) bufer)[col] == dnullval) {
- nodatavalmatch = 1;
+ else {
+ if (((DCELL *) bufer)[col] == dnullval) {
+ nodatavalmatch = 1;
+ }
+ if (dfCellMin > ((DCELL *) bufer)[col])
+ dfCellMin = ((DCELL *) bufer)[col];
+ if (dfCellMax < ((DCELL *) bufer)[col])
+ dfCellMax = ((DCELL *) bufer)[col];
}
if (GDALRasterIO
@@ -257,6 +288,8 @@
GDALDataType datatype = GDT_Int32;
CELL inullval = (CELL) nodataval;
+ G_debug(1, "CELL nodata val: %d", inullval);
+
for (row = 0; row < rows; row++) {
if (G_get_raster_row(fd, bufer, row, maptype) < 0) {
@@ -273,8 +306,14 @@
}
n_nulls++;
}
- else if (((CELL *) bufer)[col] == inullval) {
- nodatavalmatch = 1;
+ else {
+ if (((CELL *) bufer)[col] == inullval) {
+ nodatavalmatch = 1;
+ }
+ if (dfCellMin > ((CELL *) bufer)[col])
+ dfCellMin = ((CELL *) bufer)[col];
+ if (dfCellMax < ((CELL *) bufer)[col])
+ dfCellMax = ((CELL *) bufer)[col];
}
if (GDALRasterIO
@@ -287,35 +326,144 @@
}
}
- if (n_nulls > 0 && default_nodataval) {
+ /* can the GDAL datatype hold the data range to be exported ? */
+ /* f-flag does not override */
+ if (exact_range_check(export_datatype, dfCellMin, dfCellMax, name)) {
+ G_warning("Raster export results in data loss.");
+ ret = -3;
+ }
+
+ /* a default nodata value was used and NULL cells were present */
+ if (n_nulls && default_nodataval) {
if (maptype == CELL_TYPE)
G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
- "The value %d was used to represent no-data values in the input map. "
- "You can specify a nodata value with the %s option."),
- (int)nodataval, nodatakey);
+ "The value %d was used to represent no-data values in the input map. "
+ "You can specify a nodata value with the %s option."),
+ (int)nodataval, nodatakey);
else
G_important_message(_("Input raster map contains cells with NULL-value (no-data). "
- "The value %g was used to represent no-data values in the input map. "
- "You can specify a nodata value with the %s option."),
- nodataval, nodatakey);
+ "The value %f was used to represent no-data values in the input map. "
+ "You can specify a nodata value with the %s option."),
+ nodataval, nodatakey);
}
+ /* the nodata value was present in the exported data */
if (nodatavalmatch && n_nulls) {
- if (default_nodataval) { /* default nodataval didn't work */
+ /* default nodataval didn't work */
+ if (default_nodataval) {
G_warning(_("The default nodata value is present in raster"
- "band <%s> and would lead to data loss. Please specify a "
- "different nodata value with the %s parameter."),
- name, nodatakey);
+ "band <%s> and would lead to data loss. Please specify a "
+ "custom nodata value with the %s parameter."),
+ name, nodatakey);
}
- else { /* user-specified nodataval didn't work */
+ /* user-specified nodataval didn't work */
+ else {
G_warning(_("The given nodata value is present in raster"
- "band <%s> and would lead to data loss. Please specify a "
- "different nodata value with the %s parameter."),
- name, nodatakey);
+ "band <%s> and would lead to data loss. Please specify a "
+ "different nodata value with the %s parameter."),
+ name, nodatakey);
}
-
- return -2;
+ ret = -2;
}
- return 0;
+ return ret;
}
+
+int exact_range_check(double min, double max, GDALDataType datatype,
+ const char *name)
+{
+
+ switch (datatype) {
+ case GDT_Byte:
+ if (min < TYPE_BYTE_MIN || max > TYPE_BYTE_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_BYTE_MIN,
+ TYPE_BYTE_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt16:
+ if (min < TYPE_UINT16_MIN || max > TYPE_UINT16_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_UINT16_MIN,
+ TYPE_UINT16_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int16:
+ case GDT_CInt16:
+ if (min < TYPE_INT16_MIN || max > TYPE_INT16_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_INT16_MIN,
+ TYPE_INT16_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int32:
+ case GDT_CInt32:
+ if (min < TYPE_INT32_MIN || max > TYPE_INT32_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_INT32_MIN,
+ TYPE_INT32_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt32:
+ if (min < TYPE_UINT32_MIN || max > TYPE_UINT32_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %u - %u"),
+ GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
+ TYPE_UINT32_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float32:
+ case GDT_CFloat32:
+ if (min < TYPE_FLOAT32_MIN || max > TYPE_FLOAT32_MAX) {
+ G_warning(_("Selected GDAL datatype does not cover data range."));
+ G_warning(_("GDAL datatype: %s, range: %f - %f"),
+ GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
+ TYPE_FLOAT32_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float64:
+ case GDT_CFloat64:
+ /* not possible because DCELL is FLOAT64, not 128bit floating point, but anyway... */
+ 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: %f - %f"),
+ GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
+ TYPE_FLOAT64_MAX);
+ G_warning(_("Raster map <%s> range: %f - %f"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ default:
+ return 0;
+ }
+}
Modified: grass/trunk/raster/r.out.gdal/local_proto.h
===================================================================
--- grass/trunk/raster/r.out.gdal/local_proto.h 2009-06-06 12:31:34 UTC (rev 37762)
+++ grass/trunk/raster/r.out.gdal/local_proto.h 2009-06-06 16:36:19 UTC (rev 37763)
@@ -51,15 +51,11 @@
#define TYPE_FLOAT64_MAX 1.79E308f
#endif
-#define GRASS_MAX_COLORS TYPE_UINT16_MAX
+#define GRASS_MAX_COLORS TYPE_UINT16_MAX /* ok? */
-/* main.c */
-int range_check(double, double, GDALDataType, char *);
-int nullvalue_check(double, GDALDataType);
-
/* export_band.c */
-int export_band(GDALDatasetH, int, const char *, const char *,
- struct Cell_head *, RASTER_MAP_TYPE, double,
- const char *, int, int);
+int export_band(GDALDatasetH, GDALDataType, int, const char *,
+ const char *, struct Cell_head *, RASTER_MAP_TYPE,
+ double, const char *, int, int);
#endif /* __LOCAL_PROTO_H__ */
Modified: grass/trunk/raster/r.out.gdal/main.c
===================================================================
--- grass/trunk/raster/r.out.gdal/main.c 2009-06-06 12:31:34 UTC (rev 37762)
+++ grass/trunk/raster/r.out.gdal/main.c 2009-06-06 16:36:19 UTC (rev 37763)
@@ -7,7 +7,7 @@
* based on GDAL library.
* Replaces r.out.gdal.sh script which used the gdal_translate
* executable and GDAL grass-format plugin.
-* COPYRIGHT: (C) 2006-2008 by the GRASS Development Team
+* COPYRIGHT: (C) 2006-2009 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -32,13 +32,15 @@
#include "cpl_string.h"
#include "local_proto.h"
+int range_check(double, double, GDALDataType);
+int nodataval_check(double, GDALDataType);
+double set_default_nodata_value(GDALDataType, double, double);
void supported_formats(const char **formats)
{
/* Code taken from r.in.gdal */
int iDr;
-
dbString gdal_formats;
db_init_string(&gdal_formats);
@@ -111,12 +113,12 @@
struct Cell_head cellhead;
struct Ref ref;
const char *mapset, *gdal_formats = NULL;
- RASTER_MAP_TYPE maptype;
+ RASTER_MAP_TYPE maptype, testmaptype;
int bHaveMinMax;
- double dfCellMin;
- double dfCellMax;
+ double dfCellMin, export_min;
+ double dfCellMax, export_max;
struct FPRange sRange;
- int retval;
+ int retval = 0, check_range;
G_gisinit(argv[0]);
@@ -138,7 +140,7 @@
flag_f = G_define_flag();
flag_f->key = 'f';
flag_f->label = _("Force raster export also if data loss may occur");
- flag_f->description = _("Overrides saftey checks.");
+ flag_f->description = _("Overrides nodata saftey check.");
input = G_define_standard_option(G_OPT_R_INPUT);
input->required = NO;
@@ -275,9 +277,8 @@
format->answer);
}
- /* Determine raster data type */
+ /* Determine GDAL data type */
GDALDataType datatype = GDT_Unknown;
- double nodataval;
maptype = CELL_TYPE;
@@ -286,141 +287,211 @@
if (type->answer[0] == 'B') {
datatype = GDT_Byte;
maptype = CELL_TYPE;
- nodataval = TYPE_BYTE_MIN;
}
else if (type->answer[0] == 'I') {
if (strcmp(type->answer, "Int16") == 0) {
datatype = GDT_Int16;
maptype = CELL_TYPE;
- nodataval = TYPE_INT16_MIN;
}
else if (strcmp(type->answer, "Int32") == 0) {
datatype = GDT_Int32;
maptype = CELL_TYPE;
- nodataval = TYPE_INT32_MIN;
}
}
else if (type->answer[0] == 'U') {
if (strcmp(type->answer, "UInt16") == 0) {
datatype = GDT_UInt16;
maptype = CELL_TYPE;
- nodataval = TYPE_UINT16_MIN;
}
else if (strcmp(type->answer, "UInt32") == 0) {
datatype = GDT_UInt32;
maptype = DCELL_TYPE;
- nodataval = TYPE_UINT32_MIN;
}
}
else if (type->answer[0] == 'F') {
if (strcmp(type->answer, "Float32") == 0) {
datatype = GDT_Float32;
maptype = FCELL_TYPE;
- /* nodataval = TYPE_FLOAT32_MIN; */
- nodataval = 0.0/0.0;
}
else if (strcmp(type->answer, "Float64") == 0) {
datatype = GDT_Float64;
maptype = DCELL_TYPE;
- /* nodataval = TYPE_FLOAT64_MIN; */
- nodataval = 0.0/0.0;
}
}
else if (type->answer[0] == 'C') {
if (strcmp(type->answer, "CInt16") == 0) {
datatype = GDT_CInt16;
maptype = CELL_TYPE;
- nodataval = TYPE_INT16_MIN;
}
else if (strcmp(type->answer, "CInt32") == 0) {
datatype = GDT_CInt32;
maptype = CELL_TYPE;
- nodataval = TYPE_INT32_MIN;
}
else if (strcmp(type->answer, "CFloat32") == 0) {
datatype = GDT_CFloat32;
maptype = FCELL_TYPE;
- /* nodataval = TYPE_FLOAT32_MIN; */
- nodataval = 0.0/0.0;
}
else if (strcmp(type->answer, "CFloat64") == 0) {
datatype = GDT_CFloat64;
maptype = DCELL_TYPE;
- /* nodataval = TYPE_FLOAT64_MIN; */
- nodataval = 0.0/0.0;
}
}
}
- /* If file type not set by user ... */
- /* Get min/max values. */
- if (G_read_fp_range(ref.file[0].name, ref.file[0].mapset, &sRange) == -1) {
- bHaveMinMax = FALSE;
+ /* get min/max values */
+ int band;
+
+ check_range = 0;
+ bHaveMinMax = TRUE;
+ for (band = 0; band < ref.nfiles; band++) {
+ if (G_read_fp_range
+ (ref.file[band].name, ref.file[band].mapset, &sRange) == -1) {
+ bHaveMinMax = FALSE;
+ G_warning(_("Could not read data range of raster <%s>"),
+ ref.file[band].name);
+ }
+ else {
+ G_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
+ if (band == 0) {
+ export_min = dfCellMin;
+ export_max = dfCellMax;
+ }
+ else {
+ if (export_min > dfCellMin)
+ export_min = dfCellMin;
+ if (export_max < dfCellMax)
+ export_max = dfCellMax;
+ }
+
+ }
+ G_debug(3, "Range of <%s>: min: %f, max: %f", ref.file[band].name,
+ dfCellMin, dfCellMax);
+
}
- else {
- bHaveMinMax = TRUE;
- G_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
+ if (bHaveMinMax == FALSE) {
+ export_min = TYPE_FLOAT64_MIN;
+ export_max = TYPE_FLOAT64_MAX;
}
- G_debug(3, "Range: min: %f, max: %f", dfCellMin, dfCellMax);
+ G_debug(3, "Total range: min: %f, max: %f", export_min, export_max);
+ /* GDAL datatype not set by user, determine suitable datatype */
if (datatype == GDT_Unknown) {
- /* ... determine raster data type from first GRASS raster in a group */
+ /* Use raster data type from first GRASS raster in a group */
maptype = G_raster_map_type(ref.file[0].name, ref.file[0].mapset);
if (maptype == FCELL_TYPE) {
datatype = GDT_Float32;
- /* nodataval = TYPE_FLOAT32_MIN; */
- nodataval = 0.0/0.0;
}
else if (maptype == DCELL_TYPE) {
datatype = GDT_Float64;
- /* nodataval = TYPE_FLOAT64_MIN; */
- nodataval = 0.0/0.0;
}
else {
/* Special tricks for GeoTIFF color table support and such */
- if (dfCellMin >= 0 && dfCellMax < 256) {
+ if (export_min >= TYPE_BYTE_MIN && export_max <= TYPE_BYTE_MAX) {
datatype = GDT_Byte;
- nodataval = TYPE_BYTE_MIN;
}
else {
- if (dfCellMin >= 0 && dfCellMax < 65536) {
+ if (export_min >= TYPE_UINT16_MIN &&
+ export_max <= TYPE_UINT16_MAX) {
datatype = GDT_UInt16;
- nodataval = TYPE_UINT16_MIN;
}
+ else if (export_min >= TYPE_INT16_MIN &&
+ export_max <= TYPE_INT16_MAX) {
+ datatype = GDT_Int16;
+ }
else {
datatype = GDT_Int32; /* need to fine tune this more? */
- nodataval = TYPE_INT32_MIN;
}
}
}
}
- /* default or user-specified nodata-value ? */
+ /* got a GDAL datatype, report to user */
+ G_message(_("Exporting to GDAL data type: %s"),
+ GDALGetDataTypeName(datatype));
+
+ G_debug(3, "Input map datatype=%s\n",
+ (maptype == CELL_TYPE ? "CELL" :
+ (maptype == DCELL_TYPE ? "DCELL" :
+ (maptype == FCELL_TYPE ? "FCELL" : "??"))));
+
+
+ /* if GDAL datatype set by user, do checks */
+ if (type->answer) {
+
+ /* Check if raster data range is outside of the range of
+ * given GDAL datatype, not even overlapping */
+ if (range_check(export_min, export_max, datatype))
+ G_fatal_error(_("Raster export would result in complete data loss, aborting."));
+
+ /* Precision tests */
+ for (band = 0; band < ref.nfiles; band++) {
+ testmaptype =
+ G_raster_map_type(ref.file[band].name, ref.file[band].mapset);
+ /* Exporting floating point rasters to some integer type ? */
+ if ((testmaptype == FCELL_TYPE || testmaptype == DCELL_TYPE) &&
+ (datatype == GDT_Byte || datatype == GDT_Int16 ||
+ datatype == GDT_UInt16 || datatype == GDT_Int32 ||
+ datatype == GDT_UInt32)) {
+ G_warning(_("Precision loss: Raster map <%s> of type %s to be exported as %s. "
+ "This can be avoided by using %s."),
+ ref.file[band].name,
+ (maptype == FCELL_TYPE ? "FCELL" : "DCELL"),
+ GDALGetDataTypeName(datatype),
+ (maptype == FCELL_TYPE ? "Float32" : "Float64"));
+ retval = -1;
+ }
+ /* Exporting CELL with large values to GDT_Float32 ? Cap is 2^24 */
+ if (testmaptype == CELL_TYPE && datatype == GDT_Float32 &&
+ (dfCellMin < -16777216 || dfCellMax > 16777216)) {
+ G_warning(_("Precision loss: The range of <%s> can not be "
+ "accurately preserved with GDAL datatype Float32. "
+ "This can be avoided by exporting to Int32 or Float64."),
+ ref.file[band].name);
+ retval = -1;
+ }
+ /* Exporting DCELL to GDT_Float32 ? */
+ if (testmaptype == DCELL_TYPE && datatype == GDT_Float32) {
+ G_warning(_("Precision loss: Float32 can not preserve the "
+ "DCELL precision of raster <%s>. "
+ "This can be avoided by using Float64"),
+ ref.file[band].name);
+ retval = -1;
+ }
+ }
+ if (retval == -1) {
+ if (flag_f->answer)
+ G_warning(_("Forcing raster export."));
+ else
+ G_fatal_error(_("Raster export aborted."));
+ }
+ }
+
+ /* Nodata value */
+ double nodataval;
int default_nodataval = 1;
+ /* User-specified nodata-value ? */
if (nodataopt->answer) {
nodataval = atof(nodataopt->answer);
default_nodataval = 0;
- /* check nodataval here, not in export_band(), because it is specified only once */
- if (nullvalue_check(nodataval, datatype)) {
+ /* Check if given nodata value can be represented by selected GDAL datatype */
+ if (nodataval_check(nodataval, datatype)) {
if (flag_f->answer)
G_warning(_("Forcing raster export."));
else
G_fatal_error(_("Raster export aborted."));
}
}
+ /* Set reasonable default nodata value */
+ else {
+ nodataval =
+ set_default_nodata_value(datatype, export_min, export_max);
+ }
- G_debug(3, "Input map datatype=%s\n",
- (maptype == CELL_TYPE ? "CELL" :
- (maptype == DCELL_TYPE ? "DCELL" :
- (maptype == FCELL_TYPE ? "FCELL" : "??"))));
- G_message(_("Exporting to GDAL data type: %s"),
- GDALGetDataTypeName(datatype));
-
/* Create dataset for output with target driver or, if needed, with in-memory driver */
char **papszOptions = NULL;
- /* parse dataset creation options */
+ /* Parse dataset creation options */
if (createopt->answer) {
int i;
@@ -478,8 +549,6 @@
AttachMetadata(hCurrDS, metaopt->answers);
/* Export to GDAL raster */
- int band;
-
for (band = 0; band < ref.nfiles; band++) {
if (ref.nfiles > 1) {
G_verbose_message(_("Exporting raster map <%s> (band %d)..."),
@@ -487,44 +556,32 @@
ref.file[band].mapset),
band + 1);
}
- /* do range check here because we don't have GDALDataType in export_band() */
- if (G_read_fp_range
- (ref.file[band].name, ref.file[band].mapset, &sRange) == -1) {
- bHaveMinMax = FALSE;
- }
- else {
- bHaveMinMax = TRUE;
- G_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
- }
- G_debug(3, "Range: min: %f, max: %f", dfCellMin, dfCellMax);
- if (bHaveMinMax == TRUE) {
- if (range_check
- (dfCellMin, dfCellMax, datatype, ref.file[band].name)) {
- if (flag_f->answer)
- G_warning(_("Forcing raster export."));
- else
- G_fatal_error(_("Raster export aborted."));
- }
- }
- /* ready to export */
retval = export_band
- (hCurrDS, band + 1, ref.file[band].name, ref.file[band].mapset,
- &cellhead, maptype, nodataval, nodataopt->key, flag_c->answer,
- default_nodataval);
+ (hCurrDS, datatype, band + 1, ref.file[band].name,
+ ref.file[band].mapset, &cellhead, maptype, nodataval,
+ nodataopt->key, flag_c->answer, default_nodataval);
+
+ /* read/write error */
if (retval == -1) {
G_warning(_("Unable to export raster map <%s>"),
ref.file[band].name);
}
+ /* nodata problem */
else if (retval == -2) {
if (flag_f->answer)
G_warning(_("Forcing raster export."));
else
G_fatal_error(_("Raster export aborted."));
}
+ /* data don't fit into range of GDAL datatype */
+ else if (retval == -3) {
+ G_fatal_error(_("Raster export aborted."));
+ }
}
- /* Finaly create user required raster format from memory raster if in-memory driver was used */
+ /* Finaly create user requested raster format from memory raster
+ * if in-memory driver was used */
if (hMEMDS) {
hDstDS =
GDALCreateCopy(hDriver, output->answer, hMEMDS, FALSE,
@@ -546,30 +603,30 @@
}
-int range_check(double min, double max, GDALDataType datatype, char *name)
+int range_check(double min, double max, GDALDataType datatype)
{
/* what accuracy to use to print min max for FLOAT32 and FLOAT64? %g enough? */
switch (datatype) {
case GDT_Byte:
- if (min < TYPE_BYTE_MIN || max > TYPE_BYTE_MAX) {
+ if (max < TYPE_BYTE_MIN || min > TYPE_BYTE_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_BYTE_MIN,
TYPE_BYTE_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
return 0;
case GDT_UInt16:
- if (min < TYPE_UINT16_MIN || max > TYPE_UINT16_MAX) {
+ if (max < TYPE_UINT16_MIN || min > TYPE_UINT16_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_UINT16_MIN,
TYPE_UINT16_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
@@ -577,12 +634,12 @@
case GDT_Int16:
case GDT_CInt16:
- if (min < TYPE_INT16_MIN || max > TYPE_INT16_MAX) {
+ if (max < TYPE_INT16_MIN || min > TYPE_INT16_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_INT16_MIN,
TYPE_INT16_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
@@ -590,24 +647,24 @@
case GDT_Int32:
case GDT_CInt32:
- if (min < TYPE_INT32_MIN || max > TYPE_INT32_MAX) {
+ if (max < TYPE_INT32_MIN || min > TYPE_INT32_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_INT32_MIN,
TYPE_INT32_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
return 0;
case GDT_UInt32:
- if (min < TYPE_UINT32_MIN || max > TYPE_UINT32_MAX) {
+ if (max < TYPE_UINT32_MIN || min > TYPE_UINT32_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %u - %u"),
GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
TYPE_UINT32_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
@@ -615,12 +672,12 @@
case GDT_Float32:
case GDT_CFloat32:
- if (min < TYPE_FLOAT32_MIN || max > TYPE_FLOAT32_MAX) {
+ if (max < TYPE_FLOAT32_MIN || min > TYPE_FLOAT32_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %g - %g"),
GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
TYPE_FLOAT32_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
@@ -629,12 +686,12 @@
case GDT_Float64:
case GDT_CFloat64:
/* not possible because DCELL is FLOAT64, not 128bit floating point, but anyway... */
- if (min < TYPE_FLOAT64_MIN || max > TYPE_FLOAT64_MAX) {
+ if (max < TYPE_FLOAT64_MIN || min > TYPE_FLOAT64_MAX) {
G_warning(_("Selected GDAL datatype does not cover data range."));
G_warning(_("GDAL datatype: %s, range: %g - %g"),
GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
TYPE_FLOAT64_MAX);
- G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ G_warning(_("Range to be exported: %f - %f"), min, max);
return 1;
}
else
@@ -645,15 +702,15 @@
}
}
-int nullvalue_check(double nodataval, GDALDataType datatype)
+int nodataval_check(double nodataval, GDALDataType datatype)
{
- /* what accuracy to use to print nodataval for FLOAT32 and FLOAT64? %g enough? */
switch (datatype) {
case GDT_Byte:
- if (nodataval < TYPE_BYTE_MIN || nodataval > TYPE_BYTE_MAX) {
- G_warning(_("Specified nodata value %d is not covered by range of selected GDAL datatype."),
- (int) nodataval);
+ if (nodataval != (double)(GByte) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %d by selected GDAL datatype."),
+ nodataval, (GByte) nodataval);
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_BYTE_MIN,
TYPE_BYTE_MAX);
@@ -663,9 +720,10 @@
return 0;
case GDT_UInt16:
- if (nodataval < TYPE_UINT16_MIN || nodataval > TYPE_UINT16_MAX) {
- G_warning(_("Specified nodata value %d is not covered by range of selected GDAL datatype."),
- (int) nodataval);
+ if (nodataval != (double)(GUInt16) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %d by selected GDAL datatype."),
+ nodataval, (GUInt16) nodataval);
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_UINT16_MIN,
TYPE_UINT16_MAX);
@@ -676,9 +734,10 @@
case GDT_Int16:
case GDT_CInt16:
- if (nodataval < TYPE_INT16_MIN || nodataval > TYPE_INT16_MAX) {
- G_warning(_("Specified nodata value %d is not covered by range of selected GDAL datatype."),
- (int) nodataval);
+ if (nodataval != (double)(GInt16) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %d by selected GDAL datatype."),
+ nodataval, (GInt16) nodataval);
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_INT16_MIN,
TYPE_INT16_MAX);
@@ -687,11 +746,25 @@
else
return 0;
+ case GDT_UInt32:
+ if (nodataval != (double)(GUInt32) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %d by selected GDAL datatype."),
+ nodataval, (GUInt32) nodataval);
+ G_warning(_("GDAL datatype: %s, range: %u - %u"),
+ GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
+ TYPE_UINT32_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
case GDT_Int32:
case GDT_CInt32:
- if (nodataval < TYPE_INT32_MIN || nodataval > TYPE_INT32_MAX) {
- G_warning(_("Specified nodata value %lld is not covered by range of selected GDAL datatype."),
- (long long) nodataval);
+ if (nodataval != (double)(GInt32) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %d by selected GDAL datatype."),
+ nodataval, (GInt32) nodataval);
G_warning(_("GDAL datatype: %s, range: %d - %d"),
GDALGetDataTypeName(datatype), TYPE_INT32_MIN,
TYPE_INT32_MAX);
@@ -700,23 +773,12 @@
else
return 0;
- case GDT_UInt32:
- if (nodataval < TYPE_UINT32_MIN || nodataval > TYPE_UINT32_MAX) {
- G_warning(_("Specified nodata value %lld is not covered by range of selected GDAL datatype."),
- (long long) nodataval);
- G_warning(_("GDAL datatype: %s, range: %u - %u"),
- GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
- TYPE_UINT32_MAX);
- return 1;
- }
- else
- return 0;
-
case GDT_Float32:
case GDT_CFloat32:
- if (nodataval < TYPE_FLOAT32_MIN || nodataval > TYPE_FLOAT32_MAX) {
- G_warning(_("Specified nodata value %g is not covered by range of selected GDAL datatype."),
- nodataval);
+ if (nodataval != (double)(float)nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %f gets converted to %f by selected GDAL datatype."),
+ nodataval, (float)nodataval);
G_warning(_("GDAL datatype: %s, range: %g - %g"),
GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
TYPE_FLOAT32_MAX);
@@ -727,18 +789,68 @@
case GDT_Float64:
case GDT_CFloat64:
- /* not possible because double is FLOAT64, not 128bit floating point */
- if (nodataval < TYPE_FLOAT64_MIN || nodataval > TYPE_FLOAT64_MAX) {
- G_warning(_("Specified nodata value %g is not covered by range of selected GDAL datatype."),
- nodataval);
- G_warning(_("GDAL datatype: %s, range: %g - %g"),
- GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
- TYPE_FLOAT64_MAX);
- return 1;
- }
+ /* not needed because double is FLOAT64, not 128bit floating point */
+ return 0;
+
+ default:
+ return 0;
+ }
+}
+
+double set_default_nodata_value(GDALDataType datatype, double min, double max)
+{
+ switch (datatype) {
+ case GDT_Byte:
+ if (max < TYPE_BYTE_MAX)
+ return (double)TYPE_BYTE_MAX;
+ else if (min > TYPE_BYTE_MIN)
+ return (double)TYPE_BYTE_MIN;
else
- return 0;
+ return (double)TYPE_BYTE_MAX;
+ case GDT_UInt16:
+ if (max < TYPE_UINT16_MAX)
+ return (double)TYPE_UINT16_MAX;
+ else if (min > TYPE_UINT16_MIN)
+ return (double)TYPE_UINT16_MIN;
+ else
+ return (double)TYPE_UINT16_MAX;
+
+ case GDT_Int16:
+ case GDT_CInt16:
+ if (min > TYPE_INT16_MIN)
+ return (double)TYPE_INT16_MIN;
+ else if (max < TYPE_INT16_MAX)
+ return (double)TYPE_INT16_MAX;
+ else
+ return (double)TYPE_INT16_MIN;
+
+ case GDT_UInt32:
+ if (max < TYPE_UINT32_MAX)
+ return (double)TYPE_UINT32_MAX;
+ else if (min > TYPE_UINT32_MIN)
+ return (double)TYPE_UINT32_MIN;
+ else
+ return (double)TYPE_UINT32_MAX;
+
+ case GDT_Int32:
+ case GDT_CInt32:
+ if (min > TYPE_INT32_MIN)
+ return (double)TYPE_INT32_MIN;
+ else if (max < TYPE_INT32_MAX)
+ return (double)TYPE_INT32_MAX;
+ else
+ return (double)TYPE_INT32_MIN;
+
+ case GDT_Float32:
+ case GDT_CFloat32:
+ return 0.0 / 0.0;
+
+ case GDT_Float64:
+ case GDT_CFloat64:
+ return 0.0 / 0.0;
+
+ /* should not happen: */
default:
return 0;
}
More information about the grass-commit
mailing list