[GRASS-SVN] r42117 -
grass/branches/releasebranch_6_4/raster/r.out.gdal
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 5 09:08:24 EDT 2010
Author: mmetz
Date: 2010-05-05 09:08:23 -0400 (Wed, 05 May 2010)
New Revision: 42117
Modified:
grass/branches/releasebranch_6_4/raster/r.out.gdal/description.html
grass/branches/releasebranch_6_4/raster/r.out.gdal/export_band.c
grass/branches/releasebranch_6_4/raster/r.out.gdal/local_proto.h
grass/branches/releasebranch_6_4/raster/r.out.gdal/main.c
Log:
r.out.gdal fix
Modified: grass/branches/releasebranch_6_4/raster/r.out.gdal/description.html
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.out.gdal/description.html 2010-05-04 19:17:13 UTC (rev 42116)
+++ grass/branches/releasebranch_6_4/raster/r.out.gdal/description.html 2010-05-05 13:08:23 UTC (rev 42117)
@@ -1,11 +1,18 @@
<h2>DESCRIPTION</h2>
<em>r.out.gdal</em> allows a user to export a GRASS raster map layer
-into any GDAL supported raster map format.
-
-For possible <em>metaopt</em> parameters see the individual
-<a href="http://www.gdal.org/formats_list.html">supported formats</a> pages
-on the GDAL website.
+into any GDAL supported raster map format. If a GRASS raster map is
+exported for a particular application, the application's native format
+would be preferrable. GeoTIFF is supported by a wide range of
+applications (see also <b>NOTES</b> on GeoTIFF below).
+<p>
+To specify multiple creation options use a comma separated list
+(<em>createopt="TFW=YES,COMPRESS=DEFLATE"</em>).
+<p>
+For possible <em>createopt</em> and <em>metaopt</em> parameters please
+consult the individual
+<a href="http://www.gdal.org/formats_list.html">supported formats</a>
+pages on the GDAL website.
The <em>createopt</em> parameter may be used to create TFW or World files
("TFW=YES","WORLDFILE=ON").
<p>
@@ -56,19 +63,7 @@
<h2>NOTES</h2>
-When writing out multi-band GeoTIFF images for users of ESRI software or
-ImageMagick, the interleaving mode should be set to "pixel" using
-<em>createopt="INTERLEAVE=PIXEL"</em>. BAND interleaving is slightly more
-efficient, but not supported by some applications.
-<!-- GDAL switched default from BAND to PIXEL interleave on 08/01/07 (r11823) -->
-This issue only arises when writing out multi-band imagery groups.
-Some software may not recognize all of the compression methods available.
-<!-- e.g. data destined for ESRI software should use COMPRESS=LZW/PACKBITS/DEFLATE ??? -->
-
<p>
-To specify multiple options use a comma separated list
-(<em>createopt="TFW=YES,COMPRESS=DEFLATE"</em>).
-<p>
Out of the GDAL data types, the closest match for GRASS CELL, FCELL and
DCELL rasters are respectively Int32, Float32 and Float64. These are not
exact equivalents, but they will preserve the maximum possible data range
@@ -80,8 +75,8 @@
check the data type and range for your GRASS raster, refer to specific
format documentation (on the <a href="http://www.gdal.org/">GDAL website</a>),
format vendor's documentation, and e.g. the Wikipedia article
-<em><a href="http://en.wikipedia.org/wiki/C_syntax#Typical_boundaries_of_primitive_integral_types">Typical
- boundaries of primitive integral types</a></em>
+<em><a href="http://en.wikipedia.org/wiki/C_syntax#Typical_boundaries_of_primitive_integral_types">
+Typical boundaries of primitive integral types</a></em>
for details.
@@ -104,16 +99,52 @@
use Byte rather than UInt16; use Int16 rather than Int32; or use Float32
rather than Float64. In addition, the COMPRESS <b>createopt</b> used can
have a very large impact on the size of the output file.
+<p>
+Some software may not recognize all of the compression methods
+available for a given file format, and certain compression methods may
+only be supported for certain data types (depends on vendor and version).
+<!-- e.g. data destined for ESRI software should use COMPRESS=LZW/PACKBITS/DEFLATE ??? -->
+<p>
+If the export settings are set such that data loss would occur in the output
+file (i.e, due to the particular choice of data type and/or file type), the
+normal behaviour of <em>r.out.gdal</em> in this case would be to issue an error
+message describing the problem and exit without exporting. The <b>-f</b> flag
+allows raster export even if some of the data loss tests are not passed, and
+warnings are issued instead of errors.
+<p>
+<em>r.out.gdal</em> exports may appear all black or gray on initial
+display in other GIS software. This is not a bug of <em>r.out.gdal</em>,
+but often caused by the default color table assigned by that software.
+The default color table may be grayscale covering the whole range of
+possible values which is very large for e.g. Int32 or Float32. E.g.
+stretching the color table to actual min/max would help (sometimes under
+symbology).
+<h4>GeoTIFF caveats</h4>
+GeoTIFF exports can only be displayed by standard image viewers
+if the GDAL data type was set to Byte and the GeoTIFF contains either
+one or three bands. All other data types and numbers of bands can be
+properly read with GIS software only. Although GeoTIFF files usually
+have a .tif extension, these files are not necessarily images but
+first of all spatial raster datasets, e.g. SRTM DEM version 4.
+<p>
+When writing out multi-band GeoTIFF images for users of ESRI software or
+ImageMagick, the interleaving mode should be set to "pixel" using
+<em>createopt="INTERLEAVE=PIXEL"</em>. BAND interleaving is slightly more
+efficient, but not supported by some applications.
+<!-- GDAL switched default from BAND to PIXEL interleave on 08/01/07 (r11823) -->
+This issue only arises when writing out multi-band imagery groups.
+<p>
+
<h4>Improving GeoTIFF compatibility</h4>
-To create a highly compatibility GeoTIFF with various other GIS software
-packages, it is recommended to keep the GeoTIFF file as simple as possible.
-You will have to experiment with which options your software is compatible
-with, as this varies widely between vendors and versions. Long term, the
-less metadata you have to remove the more self-documenting (and useful)
-the dataset will be.
+To create a GeoTIFF that is highly compatible with various other GIS
+software packages, it is recommended to keep the GeoTIFF file as simple
+as possible. You will have to experiment with which options your
+software is compatible with, as this varies widely between vendors and
+versions. Long term, the less metadata you have to remove the more
+self-documenting (and useful) the dataset will be.
<p>
Here are some things to try:
@@ -127,8 +158,10 @@
compressed with software like <tt>zip</tt>, <tt>gnuzip</tt>, or <tt>bzip2</tt>.
<li>Skip exporting the color table. Color tables are not always properly
-rendered and the GeoTIFF file can appear completely black. If you are lucky
-the problematic software package has a method to reset the color table.
+rendered, particularly for type UInt16, and the GeoTIFF file can appear
+completely black. If you are lucky the problematic software package has
+a method to reset the color table and assign a new color table
+(sometimes called symbology).
<li>Keep metadata simple with <tt>createopt="PROFILE=GeoTIFF"</tt> or
<tt>createopt="PROFILE=BASELINE"</tt>. With BASELINE no GDAL or GeoTIFF
@@ -150,14 +183,14 @@
<h4>Export a DCELL raster map in GeoTIFF format suitable for ESRI software:</h4>
<div class="code"><pre>
-r.out.gdal in=elevation.10m out=ned_elev10m.tif type=Float64 createopt="TFW=YES"
+r.out.gdal in=elevation.10m out=ned_elev10m.tif type=Float64 createopt="PROFILE=GeoTIFF,TFW=YES"
</pre></div>
<p>
<h4>Export R,G,B imagery bands in GeoTIFF format suitable for ESRI software:</h4>
<div class="code"><pre>
i.group group=nc_landsat_rgb input=lsat7_2002_30,lsat7_2002_20,lsat7_2002_10
-r.out.gdal in=nc_landsat_rgb out=nc_landsat_rgb.tif type=Byte createopt="INTERLEAVE=PIXEL,TFW=YES"
+r.out.gdal in=nc_landsat_rgb out=nc_landsat_rgb.tif type=Byte createopt="PROFILE=GeoTIFF,INTERLEAVE=PIXEL,TFW=YES"
</pre></div>
<p>
@@ -216,9 +249,10 @@
GDAL Pages: <a href="http://www.gdal.org">http://www.gdal.org</a>
-<h2>AUTHOR</h2>
+<h2>AUTHORS</h2>
-Vytautas Vebra (oliver4grass at gmail.com)
+Vytautas Vebra (oliver4grass at gmail.com)<br>
+Markus Metz (improved nodata logic)
<p>
<i>Last changed: $Date$</i>
Modified: grass/branches/releasebranch_6_4/raster/r.out.gdal/export_band.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.out.gdal/export_band.c 2010-05-04 19:17:13 UTC (rev 42116)
+++ grass/branches/releasebranch_6_4/raster/r.out.gdal/export_band.c 2010-05-05 13:08:23 UTC (rev 42117)
@@ -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
@@ -16,18 +16,223 @@
#include <grass/gis.h>
#include <grass/glocale.h>
-
#include "cpl_string.h"
#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)
+/* exact check for each band
+ * returns 0 on success
+ * -1 if given nodata value was present in data
+ * -2 if selected GDAL datatype could not hold all values
+ * */
+int exact_checks(GDALDataType export_datatype,
+ const char *name, const char *mapset,
+ struct Cell_head *cellhead, RASTER_MAP_TYPE maptype,
+ double nodataval, const char *nodatakey,
+ int default_nodataval)
{
+ int bHaveMinMax;
+ double dfCellMin;
+ double dfCellMax;
+ struct FPRange sRange;
+ int fd;
+ int cols = cellhead->cols;
+ int rows = cellhead->rows;
+ int ret = 0;
+ /* Open GRASS raster */
+ fd = G_open_cell_old(name, mapset);
+ if (fd < 0) {
+ G_warning(_("Unable to open raster map <%s>"), name);
+ return -1;
+ }
+
+ /* Get min/max values. */
+ if (G_read_fp_range(name, mapset, &sRange) == -1) {
+ bHaveMinMax = FALSE;
+ }
+ else {
+ bHaveMinMax = TRUE;
+ G_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
+ }
+
+ /* Create GRASS raster buffer */
+ void *bufer = G_allocate_raster_buf(maptype);
+
+ if (bufer == NULL) {
+ G_warning(_("Unable to allocate buffer for reading raster map"));
+ return -1;
+ }
+ char *nulls = (char *)G_malloc(cols);
+
+ if (nulls == NULL) {
+ G_warning(_("Unable to allocate buffer for reading raster map"));
+ return -1;
+ }
+
+ /* 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) {
+
+ 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) {
+ G_warning(_("Unable to read raster map <%s> row %d"),
+ name, row);
+ return -1;
+ }
+ G_get_null_value_row(fd, nulls, row);
+ for (col = 0; col < cols; col++) {
+ if (nulls[col]) {
+ ((FCELL *) bufer)[col] = fnullval;
+ n_nulls++;
+ }
+ 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];
+ }
+ }
+ G_percent(row + 1, rows, 2);
+ }
+ }
+ else if (maptype == DCELL_TYPE) {
+
+ 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) {
+ G_warning(_("Unable to read raster map <%s> row %d"),
+ name, row);
+ return -1;
+ }
+ G_get_null_value_row(fd, nulls, row);
+ for (col = 0; col < cols; col++) {
+ if (nulls[col]) {
+ ((DCELL *) bufer)[col] = dnullval;
+ n_nulls++;
+ }
+ 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];
+ }
+ }
+ G_percent(row + 1, rows, 2);
+ }
+ }
+ else {
+
+ 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) {
+ G_warning(_("Unable to read raster map <%s> row %d"),
+ name, row);
+ return -1;
+ }
+ G_get_null_value_row(fd, nulls, row);
+ for (col = 0; col < cols; col++) {
+ if (nulls[col]) {
+ ((CELL *) bufer)[col] = inullval;
+ n_nulls++;
+ }
+ 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];
+ }
+ }
+ G_percent(row + 1, rows, 2);
+ }
+ }
+ G_debug(1, "min %g max %g", dfCellMin, dfCellMax);
+
+ /* can the GDAL datatype hold the data range to be exported ? */
+ /* f-flag does not override */
+ if (exact_range_check(dfCellMin, dfCellMax, export_datatype, name)) {
+ G_warning("Raster export results in data loss.");
+ ret = -2;
+ }
+
+ /* 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 will be 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 %f will be 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) {
+ /* 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 "
+ "custom nodata value with the %s parameter."),
+ name, nodatakey);
+ }
+ /* 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);
+ }
+ ret = -1;
+ }
+
+ return ret;
+}
+
+/* actual raster band export
+ * returns 0 on success
+ * -1 on raster data read/write error
+ * */
+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;
GDALColorTableH hCT;
int iColor;
@@ -38,6 +243,7 @@
int fd;
int cols = cellhead->cols;
int rows = cellhead->rows;
+ int ret = 0;
/* Open GRASS raster */
fd = G_open_cell_old(name, mapset);
@@ -63,14 +269,12 @@
G_get_fp_range_min_max(&sRange, &dfCellMin, &dfCellMax);
}
- /* TODO: if data range exceeds export TYPE range exit with an error.
- Otherwise the module doesn't know what to write for those values */
-
/* suppress useless warnings */
CPLPushErrorHandler(CPLQuietErrorHandler);
GDALSetRasterColorInterpretation(hBand, GPI_RGB);
CPLPopErrorHandler();
+ /* use default color rules if no color rules are given */
if (G_read_colors(name, mapset, &sGrassColors) >= 0) {
int maxcolor, i;
CELL min, max;
@@ -114,8 +318,8 @@
int nRed, nGreen, nBlue;
GDALColorEntry sColor;
- if (G_get_color
- (iColor, &nRed, &nGreen, &nBlue, &sGrassColors)) {
+ if (G_get_color(iColor, &nRed, &nGreen, &nBlue,
+ &sGrassColors)) {
sColor.c1 = nRed;
sColor.c2 = nGreen;
sColor.c3 = nBlue;
@@ -138,6 +342,8 @@
GDALSetColorEntry(hCT, iColor, &sColor);
}
}
+
+ GDALSetRasterColorTable(hBand, hCT);
}
if (rcount > 0) {
@@ -148,6 +354,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;
@@ -161,17 +369,7 @@
r2, g2, b2);
GDALSetMetadataItem(hBand, key, value, NULL);
}
-
- if (!suppress_main_colortable)
- GDALSetRasterColorTable(hBand, hCT);
}
- else {
- if (!suppress_main_colortable) {
- hCT = GDALCreateColorTable(GPI_RGB);
- GDALSetMetadataItem(hBand, "COLOR_TABLE_RULES_COUNT", "0", NULL);
- GDALSetRasterColorTable(hBand, hCT);
- }
- }
/* Create GRASS raster buffer */
void *bufer = G_allocate_raster_buf(maptype);
@@ -187,16 +385,21 @@
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;
+ /* Better use selected GDAL datatype instead of
+ * the best match with GRASS raster map types ? */
+
if (maptype == FCELL_TYPE) {
- /* Source datatype understandible by GDAL */
+ /* 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,7 +408,7 @@
return -1;
}
G_get_null_value_row(fd, nulls, row);
- for (col = 0; col < cols; col++)
+ for (col = 0; col < cols; col++) {
if (nulls[col]) {
((FCELL *) bufer)[col] = fnullval;
if (n_nulls == 0) {
@@ -213,6 +416,7 @@
}
n_nulls++;
}
+ }
if (GDALRasterIO
(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
@@ -228,6 +432,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) {
@@ -236,7 +442,7 @@
return -1;
}
G_get_null_value_row(fd, nulls, row);
- for (col = 0; col < cols; col++)
+ for (col = 0; col < cols; col++) {
if (nulls[col]) {
((DCELL *) bufer)[col] = dnullval;
if (n_nulls == 0) {
@@ -244,6 +450,7 @@
}
n_nulls++;
}
+ }
if (GDALRasterIO
(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
@@ -259,6 +466,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) {
@@ -267,7 +476,7 @@
return -1;
}
G_get_null_value_row(fd, nulls, row);
- for (col = 0; col < cols; col++)
+ for (col = 0; col < cols; col++) {
if (nulls[col]) {
((CELL *) bufer)[col] = inullval;
if (n_nulls == 0) {
@@ -275,6 +484,7 @@
}
n_nulls++;
}
+ }
if (GDALRasterIO
(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
@@ -286,18 +496,103 @@
}
}
- if (n_nulls > 0) { /* TODO: && nodata_param NOT specified */
- 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 nodata value by %s parameter."),
- (int)nodataval, nodatakey);
+ 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: %g - %g"), name, min, max);
+ return 1;
+ }
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 nodata value by %s parameter."),
- nodataval, nodatakey);
+ 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: %g - %g"), 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: %g - %g"), 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: %g - %g"), 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: %g - %g"), 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: %g - %g"),
+ GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
+ TYPE_FLOAT32_MAX);
+ G_warning(_("Raster map <%s> range: %g - %g"), 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: %g - %g"),
+ GDALGetDataTypeName(datatype), TYPE_FLOAT64_MIN,
+ TYPE_FLOAT64_MAX);
+ G_warning(_("Raster map <%s> range: %g - %g"), name, min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ default:
+ return 0;
}
-
- return 0;
}
Modified: grass/branches/releasebranch_6_4/raster/r.out.gdal/local_proto.h
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.out.gdal/local_proto.h 2010-05-04 19:17:13 UTC (rev 42116)
+++ grass/branches/releasebranch_6_4/raster/r.out.gdal/local_proto.h 2010-05-05 13:08:23 UTC (rev 42117)
@@ -1,10 +1,8 @@
#ifndef __LOCAL_PROTO_H__
#define __LOCAL_PROTO_H__
+#include "gdal.h"
-#define GRASS_MAX_COLORS 100000 /* what is the right value? UInt16 -> 65535 ? */
-/* TODO: better- set from lesser of above and TYPE limits listed below */
-
/* range limits */
/*
GDAL data type minimum maximum
@@ -18,38 +16,49 @@
Float64, CFloat64 -1.79E308 1.79E308
*/
-/* TODO: for data export range checks and nodata value check
+/* copied from limits.h, checked against gdal-1.6.0/gcore/gdalrasterband.cpp */
#define TYPE_BYTE_MIN 0
#define TYPE_BYTE_MAX 255
-#define TYPE_INT16_MIN -32768
+#define TYPE_INT16_MIN (-32768)
#define TYPE_INT16_MAX 32767
#define TYPE_UINT16_MIN 0
#define TYPE_UINT16_MAX 65535
#define TYPE_UINT32_MIN 0
-#define TYPE_UINT32_MAX 4294967295 // better to use (double)(unsigned int)0xFFFFFFFFu ?
-#define TYPE_INT32_MIN -2147483648
+#define TYPE_UINT32_MAX 4294967295U
+#define TYPE_INT32_MIN (-TYPE_INT32_MAX - 1)
#define TYPE_INT32_MAX 2147483647
+
+/* new systems: FLT_MAX, DBL_MAX, old systems: MAXFLOAT, MAXDOUBLE, fallback: 3.4E38 and 1.79E308f */
+#ifdef FLT_MAX
+#define TYPE_FLOAT32_MIN (-FLT_MAX)
+#define TYPE_FLOAT32_MAX FLT_MAX
+#elif defined(MAX_FLOAT)
+#define TYPE_FLOAT32_MIN (-MAXFLOAT)
+#define TYPE_FLOAT32_MAX MAXFLOAT
+#else
#define TYPE_FLOAT32_MIN -3.4E38f
#define TYPE_FLOAT32_MAX 3.4E38f
-#define TYPE_FLOAT64_MIN -1.79E308f
-#define TYPE_FLOAT64_MAX 1.79E308f
+#endif
- better to not define Ctypes here, rather in the code use
- if( type == F16 || type== CF16 )
-#define TYPE_CINT16_MIN TYPE_INT16_MIN
-#define TYPE_CINT16_MAX TYPE_INT16_MAX
-#define TYPE_CINT32_MIN TYPE_INT32_MIN
-#define TYPE_CINT32_MAX TYPE_INT32_MAX
-#define TYPE_CFLOAT32_MIN TYPE_FLOAT32_MIN
-#define TYPE_CFLOAT32_MAX TYPE_FLOAT32_MAX
-#define TYPE_CFLOAT64_MIN TYPE_FLOAT64_MIN
-#define TYPE_CFLOAT64_MAX TYPE_FLOAT64_MAX
-*/
+#ifdef DBL_MAX
+#define TYPE_FLOAT64_MIN (-DBL_MAX)
+#define TYPE_FLOAT64_MAX DBL_MAX
+#elif defined(MAXDOUBLE)
+#define TYPE_FLOAT64_MIN (-MAXDOUBLE)
+#define TYPE_FLOAT64_MAX MAXDOUBLE
+#else
+#define TYPE_FLOAT64_MIN -1.79E308
+#define TYPE_FLOAT64_MAX 1.79E308
+#endif
+#define GRASS_MAX_COLORS TYPE_UINT16_MAX /* ok? */
/* export_band.c */
-int export_band(GDALDatasetH, int, const char *, const char *,
- struct Cell_head *, RASTER_MAP_TYPE, double,
- const char *, int);
+int export_band(GDALDatasetH, GDALDataType, int, const char *,
+ const char *, struct Cell_head *, RASTER_MAP_TYPE,
+ double, const char *, int, int);
+int exact_checks(GDALDataType, const char *, const char *,
+ struct Cell_head *, RASTER_MAP_TYPE, double,
+ const char *, int);
#endif /* __LOCAL_PROTO_H__ */
Modified: grass/branches/releasebranch_6_4/raster/r.out.gdal/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.out.gdal/main.c 2010-05-04 19:17:13 UTC (rev 42116)
+++ grass/branches/releasebranch_6_4/raster/r.out.gdal/main.c 2010-05-05 13:08:23 UTC (rev 42117)
@@ -2,12 +2,12 @@
/****************************************************************************
*
* MODULE: r.out.gdal
-* AUTHOR(S): Vytautas Vebra <olivership at gmail.com>
+* AUTHOR(S): Vytautas Vebra <olivership at gmail.com>, Markus Metz
* PURPOSE: Exports GRASS raster to GDAL suported formats;
* 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
@@ -30,16 +30,17 @@
#include <grass/dbmi.h>
#include "cpl_string.h"
-#include "gdal.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(char **formats)
+void supported_formats(const char **formats)
{
/* Code taken from r.in.gdal */
int iDr;
-
dbString gdal_formats;
db_init_string(&gdal_formats);
@@ -105,18 +106,19 @@
{
struct GModule *module;
- struct Flag *flag_l, *flag_c;
+ struct Flag *flag_l, *flag_c, *flag_f;
struct Option *input, *format, *type, *output, *createopt, *metaopt,
*nodataopt;
struct Cell_head cellhead;
struct Ref ref;
- char *mapset, *gdal_formats = NULL;
- RASTER_MAP_TYPE maptype;
+ const char *mapset, *gdal_formats = NULL;
+ RASTER_MAP_TYPE maptype, testmaptype;
int bHaveMinMax;
- double dfCellMin;
- double dfCellMax;
+ double dfCellMin, export_min;
+ double dfCellMax, export_max;
struct FPRange sRange;
+ int retval = 0, check_range;
G_gisinit(argv[0]);
@@ -132,9 +134,14 @@
flag_c = G_define_flag();
flag_c->key = 'c';
- flag_c->label = _("Do not export long colortable");
+ flag_c->label = _("Do not write GDAL standard colortable");
flag_c->description = _("Only applicable to Byte or UInt16 data types.");
+ flag_f = G_define_flag();
+ flag_f->key = 'f';
+ flag_f->label = _("Force raster export despite any warnings of data loss");
+ flag_f->description = _("Overrides nodata safety check.");
+
input = G_define_standard_option(G_OPT_R_INPUT);
input->required = NO;
input->description = _("Name of raster map (or group) to export");
@@ -270,139 +277,249 @@
format->answer);
}
- /* Determine raster data type */
+ /* Determine GDAL data type */
GDALDataType datatype = GDT_Unknown;
- double nodataval;
+ maptype = CELL_TYPE;
+
if (type->answer) {
/* reduce number of strcmps ... */
if (type->answer[0] == 'B') {
datatype = GDT_Byte;
maptype = CELL_TYPE;
- nodataval = (double)(unsigned char)0xFFu;
}
else if (type->answer[0] == 'I') {
if (strcmp(type->answer, "Int16") == 0) {
datatype = GDT_Int16;
maptype = CELL_TYPE;
- nodataval = (double)(short)0x8000;
}
else if (strcmp(type->answer, "Int32") == 0) {
datatype = GDT_Int32;
maptype = CELL_TYPE;
- nodataval = (double)(int)0x80000000;
}
}
else if (type->answer[0] == 'U') {
if (strcmp(type->answer, "UInt16") == 0) {
datatype = GDT_UInt16;
maptype = CELL_TYPE;
- nodataval = (double)(unsigned short)0xFFFFu;
}
else if (strcmp(type->answer, "UInt32") == 0) {
datatype = GDT_UInt32;
maptype = DCELL_TYPE;
- nodataval = (double)(unsigned int)0xFFFFFFFFu;
}
}
else if (type->answer[0] == 'F') {
if (strcmp(type->answer, "Float32") == 0) {
datatype = GDT_Float32;
maptype = FCELL_TYPE;
- nodataval = -1E37f;
}
else if (strcmp(type->answer, "Float64") == 0) {
datatype = GDT_Float64;
maptype = DCELL_TYPE;
- nodataval = -1E37f;
}
}
else if (type->answer[0] == 'C') {
if (strcmp(type->answer, "CInt16") == 0) {
datatype = GDT_CInt16;
maptype = CELL_TYPE;
- nodataval = (double)(short)0x8000;
}
else if (strcmp(type->answer, "CInt32") == 0) {
datatype = GDT_CInt32;
maptype = CELL_TYPE;
- nodataval = (double)(int)0x80000000;
}
else if (strcmp(type->answer, "CFloat32") == 0) {
datatype = GDT_CFloat32;
maptype = FCELL_TYPE;
- nodataval = -1E37f;
}
else if (strcmp(type->answer, "CFloat64") == 0) {
datatype = GDT_CFloat64;
maptype = DCELL_TYPE;
- nodataval = -1E37f;
}
}
}
- /* 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;
+ export_min = TYPE_FLOAT64_MIN;
+ export_max = TYPE_FLOAT64_MAX;
+ 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 = -1E37f;
}
else if (maptype == DCELL_TYPE) {
datatype = GDT_Float64;
- nodataval = -1E37f;
}
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 = (double)(unsigned char)0xFFu;
}
else {
- if (dfCellMin >= 0 && dfCellMax < 65536) {
+ if (export_min >= TYPE_UINT16_MIN &&
+ export_max <= TYPE_UINT16_MAX) {
datatype = GDT_UInt16;
- nodataval = (double)(short)0x8000;
}
+ 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 = (double)(int)0x80000000;
}
}
}
}
- /* force nodata-value if needed */
- if (nodataopt->answer) {
- nodataval = atof(nodataopt->answer);
- }
+ /* got a GDAL datatype, report to user */
+ G_message(_("Exporting to GDAL data type: %s"),
+ GDALGetDataTypeName(datatype));
- /* TODO:
- if( nodata_used && nodataval is out_of_type_range )
- G_fatal_error(_("No-data value out of range for type %s"), type->answer);
- */
-
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));
+
+ /* 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;
+ }
+ }
+ /* disabled to not break R-GRASS interface in 6.4 */
+ /*
+ 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 if given nodata value can be represented by selected GDAL datatype */
+ if (nodataval_check(nodataval, datatype)) {
+ G_fatal_error(_("Raster export aborted."));
+ }
+ }
+ /* Set reasonable default nodata value */
+ else {
+ nodataval =
+ set_default_nodata_value(datatype, export_min, export_max);
+ }
+
+ /* exact range and nodata checks for each band */
+ G_message(_("Checking GDAl data type and nodata value"));
+ for (band = 0; band < ref.nfiles; band++) {
+ if (ref.nfiles > 1) {
+ G_verbose_message(_("Checking options for raster map <%s> (band %d)..."),
+ G_fully_qualified_name(ref.file[band].name,
+ ref.file[band].mapset),
+ band + 1);
+ }
+
+ retval = exact_checks
+ (datatype, ref.file[band].name, ref.file[band].mapset,
+ &cellhead, maptype, nodataval, nodataopt->key,
+ default_nodataval);
+
+ /* nodata value is present in the data to be exported */
+ if (retval == -1) {
+ 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 == -2) {
+ G_fatal_error(_("Raster export aborted."));
+ }
+ }
+
/* 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;
@@ -460,8 +577,7 @@
AttachMetadata(hCurrDS, metaopt->answers);
/* Export to GDAL raster */
- int band;
-
+ G_message(_("Exporting to GDAL raster"));
for (band = 0; band < ref.nfiles; band++) {
if (ref.nfiles > 1) {
G_verbose_message(_("Exporting raster map <%s> (band %d)..."),
@@ -469,15 +585,21 @@
ref.file[band].mapset),
band + 1);
}
- if (export_band
- (hCurrDS, band + 1, ref.file[band].name, ref.file[band].mapset,
- &cellhead, maptype, nodataval, nodataopt->key,
- flag_c->answer) < 0)
+
+ retval = export_band
+ (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);
+ }
}
- /* 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,
@@ -497,3 +619,250 @@
G_done_msg(" ");
exit(EXIT_SUCCESS);
}
+
+
+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 (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(_("Range to be exported: %f - %f"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt16:
+ 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(_("Range to be exported: %f - %f"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int16:
+ case GDT_CInt16:
+ 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(_("Range to be exported: %f - %f"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int32:
+ case GDT_CInt32:
+ 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(_("Range to be exported: %f - %f"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt32:
+ 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(_("Range to be exported: %f - %f"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float32:
+ case GDT_CFloat32:
+ 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(_("Range to be exported: %g - %g"), min, max);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float64:
+ case GDT_CFloat64:
+ /* not needed because FLOAT64 should always cover the data range */
+ return 0;
+
+ default:
+ return 0;
+ }
+}
+
+int nodataval_check(double nodataval, GDALDataType datatype)
+{
+
+ switch (datatype) {
+ case GDT_Byte:
+ /* the additional cast to CELL is what happens in export_band()
+ * accordingly below for the other GDAL types */
+ if (nodataval != (double)(GByte)(CELL) 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)(CELL) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_BYTE_MIN,
+ TYPE_BYTE_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt16:
+ if (nodataval != (double)(GUInt16)(CELL) 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)(CELL) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %u - %u"),
+ GDALGetDataTypeName(datatype), TYPE_UINT16_MIN,
+ TYPE_UINT16_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int16:
+ case GDT_CInt16:
+ if (nodataval != (double)(GInt16)(CELL) 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)(CELL) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_INT16_MIN,
+ TYPE_INT16_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_UInt32:
+ if (nodataval != (double)(GUInt32)(DCELL) 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)(DCELL) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %u - %u"),
+ GDALGetDataTypeName(datatype), TYPE_UINT32_MIN,
+ TYPE_UINT32_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Int32:
+ case GDT_CInt32:
+ /* GInt32 is equal to CELL, but that may change in the future */
+ if (nodataval != (double)(GInt32)(CELL) 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)(CELL) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %d - %d"),
+ GDALGetDataTypeName(datatype), TYPE_INT32_MIN,
+ TYPE_INT32_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float32:
+ case GDT_CFloat32:
+ if (nodataval != (double)(float) nodataval) {
+ G_warning(_("Mismatch between metadata nodata value and actual nodata value in exported raster: "
+ "specified nodata value %g gets converted to %g by selected GDAL datatype."),
+ nodataval, (float) nodataval);
+ G_warning(_("GDAL datatype: %s, valid range: %g - %g"),
+ GDALGetDataTypeName(datatype), TYPE_FLOAT32_MIN,
+ TYPE_FLOAT32_MAX);
+ return 1;
+ }
+ else
+ return 0;
+
+ case GDT_Float64:
+ case GDT_CFloat64:
+ /* not needed because FLOAT64 is equal to double */
+ 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 (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;
+
+ default:
+ return 0;
+ }
+}
More information about the grass-commit
mailing list