[GRASS-SVN] r34015 - grass/branches/develbranch_6/raster/r.out.gdal

svn_grass at osgeo.org svn_grass at osgeo.org
Sun Oct 26 23:56:42 EDT 2008


Author: hamish
Date: 2008-10-26 23:56:41 -0400 (Sun, 26 Oct 2008)
New Revision: 34015

Added:
   grass/branches/develbranch_6/raster/r.out.gdal/export_band.c
   grass/branches/develbranch_6/raster/r.out.gdal/local_proto.h
Modified:
   grass/branches/develbranch_6/raster/r.out.gdal/description.html
   grass/branches/develbranch_6/raster/r.out.gdal/main.c
Log:
cleanups in pursuit of bug #73 in collaboration with Markus Metz
* rename & move export band to its own file
* flag to disable writing long (UInt16) color tables
* add type length defines (needs competion & review)
* add some TODO comments where future work is required
* help page cleanups


Modified: grass/branches/develbranch_6/raster/r.out.gdal/description.html
===================================================================
--- grass/branches/develbranch_6/raster/r.out.gdal/description.html	2008-10-26 23:52:01 UTC (rev 34014)
+++ grass/branches/develbranch_6/raster/r.out.gdal/description.html	2008-10-27 03:56:41 UTC (rev 34015)
@@ -3,15 +3,24 @@
 <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 'supported formats' pages
-of GDAL.
-The <em>createopt</em> may be used to create TFW or World files ("TFW=YES",
-"WORLDFILE=ON").
+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.
+The <em>createopt</em> parameter may be used to create TFW or World files
+("TFW=YES","WORLDFILE=ON").
 <p>
-<em>r.out.gdal</em> also supports the export of multiband rasters through
-a group (created e.g. with <em>i.group</em>), when the group's name is entered
-as input.
+<em>r.out.gdal</em> also supports the export of multiband rasters as
+a group, when the imagery group's name is entered as input.
+(created imagery groups with the <em><a href="i.group.html">i.group</a></em>
+module)
+<p>
+As with most GRASS raster modules, the current region extents and region
+resolution are used, and a MASK is respected if present.
+Use <em><a href="g.region.html">g.region</a></em>'s "align=", or "rast="
+options if you need to realign the region settings the original map's
+before export.
 
+
 <h2>SUPPORTED RASTER FORMATS</h2>
 
 The set of <a href="http://www.gdal.org/formats_list.html">supported
@@ -44,29 +53,102 @@
   XPM: X11 PixMap Format
 </pre>
 
+
 <h2>NOTES</h2>
 
-When writing out GeoTIFF format for users of ESRI software or ImageMagick,
-the band interleaving should be switched to pixel interleaving using
-<em>createopt="INTERLEAVE=PIXEL"</em>.
+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 max possible data range and number of
-decimal places for each respective GRASS raster data type. Please keep in mind, that
-not all CELL rasters will require Int32 - e.g., 0-255 CELL raster are covered
-by the Byte <em>type</em> as well. Moreover, some GDAL-supported formats do not
-support all the data types possible in GDAL and GRASS. Use
-<em><a href="r.info">r.info</a></em> to check the data type and range for your
-GRASS raster, refer to specific format documentation
-(<a href="http://www.gdal.org/">GDAL website</a>, format vendor's docs) 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>
+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
+and number of decimal places for each respective GRASS raster data type.
+Please keep in mind that not all CELL rasters will require Int32 - e.g.,
+0-255 CELL raster are covered by the Byte <em>type</em> as well.
+Moreover, some GDAL-supported formats do not support all the data types
+possible in GDAL and GRASS. Use <em><a href="r.info">r.info</a></em> to
+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>
 for details.
 
+
+<h4>Ranges of GDAL data types</h4>
+<pre>
+  GDAL data type	       minimum  	maximum
+
+  Byte  			     0  	    255
+  UInt16			     0  	 65,535
+  Int16, CInt16 	       -32,768  	 32,767
+  UInt32			     0    4,294,967,295
+  Int32, CInt32 	-2,147,483,648    2,147,483,647
+  Float32, CFloat32	       -3.4E38  	 3.4E38
+  Float64, CFloat64	     -1.79E308         1.79E308
+</pre>
+
+<p>
+Float32 precision is 7 decimal places, e.g. 1.0000001E20, so 1.000000001E20
+is truncated to 1.0000000E20.
+<p>
+Float64 precision is 16 decimal places, e.g. 1.0000000000000001E20,
+equivalent to DCELL.
+<p>
+Note that e.g. 0.0000000012345678E20 is stored as 1.2345678000000000E11,
+so Float32 would be ok.
+<p>
+If there is a need to keep file sizes small, use the simplest data type
+covering the data range of the raster(s) to be exported, e.g., if suitable
+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.
+
+
+<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.
+<p>
+Here are some things to try:
+
+<ul>
+<li>Create a World file with <tt>createopt="TFW=YES"</tt>.
+
+<li>Do not use GeoTIFF internal compression. Other GIS software often 
+supports only a subset of the available compression methods with the 
+supported methods differing between GIS software packages. Unfortunately
+this means the output image can be rather huge, but the file can be
+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.
+
+<li>Keep metadata simple with <tt>createopt="PROFILE=GeoTIFF"</tt> or 
+<tt>createopt="PROFILE=BASELINE"</tt>. With BASELINE no GDAL or GeoTIFF
+tags will be written and a World file is required (<em>createopt="TFW=YES"</em>).
+
+<li>Adding overviews with <tt>gdaladdo</tt> after exporting can speed up display.
+Note that other software might create their own overviews, ignoring existing
+overviews.
+</ul>
+
+
 <h2>EXAMPLES</h2>
 
 <h4>Export the integer raster roads map to GeoTIFF format:</h4>
@@ -77,10 +159,17 @@
 
 <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="INTERLEAVE=PIXEL,TFW=YES"
+r.out.gdal in=elevation.10m out=ned_elev10m.tif type=Float64 createopt="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"
+</pre></div>
+<p>
+
 <h4>Export the floating point raster elevation map to ERDAS/IMG format:</h4>
 <div class="code"><pre>
 r.out.gdal input=elevation.10m output=elev_dem10.img format=HFA type=Float32
@@ -99,13 +188,17 @@
 <h2>GDAL RELATED ERROR MESSAGES</h2>
 
 <ul>
-<li> "ERROR 6: SetColorInterpretation() not supported for this dataset.": This <i>may</i>
- indicate that the color table was not written properly. But usually it will be
- correct and the message can be ignored.</li>
-<li> "ERROR 6: SetNoDataValue() not supported for this dataset.": The selected output format
- does not support "no data". It is recommended to use a different output format.</li>
-<li> "Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in tag.": The color table metadata
-  are too large. It is recommended to use a different output format.</li>
+<li> "ERROR 6: SetColorInterpretation() not supported for this dataset.":
+ This <i>may</i> indicate that the color table was not written properly.
+ But usually it will be correct and the message can be ignored.</li>
+
+<li> "ERROR 6: SetNoDataValue() not supported for this dataset.":
+ The selected output format does not support "no data". It is recommended
+ to use a different output format if your data contains NULLs.</li>
+
+<li> "Warning 1: Lost metadata writing to GeoTIFF ... too large to fit in
+ tag.": The color table metadata may be too large. It is recommended to
+ simplify or not write the color table, or use a different output format.</li>
 </ul>
 
 

Added: grass/branches/develbranch_6/raster/r.out.gdal/export_band.c
===================================================================
--- grass/branches/develbranch_6/raster/r.out.gdal/export_band.c	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.out.gdal/export_band.c	2008-10-27 03:56:41 UTC (rev 34015)
@@ -0,0 +1,301 @@
+/****************************************************************************
+*
+* MODULE:       r.out.gdal
+* AUTHOR(S):    Vytautas Vebra <olivership at gmail.com>
+* PURPOSE:      Exports GRASS raster to GDAL suported formats;
+*               based on GDAL library.
+*
+* COPYRIGHT:    (C) 2006-2008 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
+*   	    	for details.
+*
+*****************************************************************************/
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "cpl_string.h"
+#include "gdal.h"
+#include "local_proto.h"
+
+
+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)
+{
+
+    struct Colors sGrassColors;
+    GDALColorTableH hCT;
+    int iColor;
+    int bHaveMinMax;
+    double dfCellMin;
+    double dfCellMax;
+    struct FPRange sRange;
+    int fd;
+    int cols = cellhead->cols;
+    int rows = cellhead->rows;
+
+    /* 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 raster band  */
+    GDALRasterBandH hBand = GDALGetRasterBand(hMEMDS, band);
+
+    if (hBand == NULL) {
+	G_warning(_("Unable to get raster band"));
+	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);
+    }
+
+/* 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();
+
+    if (G_read_colors(name, mapset, &sGrassColors) >= 0) {
+	int maxcolor, i;
+	CELL min, max;
+	char key[200], value[200];
+	int rcount;
+
+	G_get_color_range(&min, &max, &sGrassColors);
+	if (bHaveMinMax) {
+	    if (max < dfCellMax) {
+		maxcolor = max;
+	    }
+	    else {
+		maxcolor = (int)ceil(dfCellMax);
+	    }
+	    if (maxcolor > GRASS_MAX_COLORS) {
+		maxcolor = GRASS_MAX_COLORS;
+		G_warning("Too many values, color table cut to %d entries",
+			  maxcolor);
+	    }
+	}
+	else {
+	    if (max < GRASS_MAX_COLORS) {
+		maxcolor = max;
+	    }
+	    else {
+		maxcolor = GRASS_MAX_COLORS;
+		G_warning("Too many values, color table set to %d entries",
+			  maxcolor);
+	    }
+	}
+
+	rcount = G_colors_count(&sGrassColors);
+
+	G_debug(3, "dfCellMin: %f, dfCellMax: %f, maxcolor: %d", dfCellMin,
+		dfCellMax, maxcolor);
+
+	if(!suppress_main_colortable) {
+	    hCT = GDALCreateColorTable(GPI_RGB);
+
+	    for (iColor = 0; iColor <= maxcolor; iColor++) {
+		int nRed, nGreen, nBlue;
+		GDALColorEntry sColor;
+
+		if (G_get_color(iColor, &nRed, &nGreen, &nBlue, &sGrassColors)) {
+		    sColor.c1 = nRed;
+		    sColor.c2 = nGreen;
+		    sColor.c3 = nBlue;
+		    sColor.c4 = 255;
+
+		    G_debug(3,
+			"G_get_color: Y, rcount %d, nRed %d, nGreen %d, nBlue %d",
+			rcount, nRed, nGreen, nBlue);
+		    GDALSetColorEntry(hCT, iColor, &sColor);
+		}
+		else {
+		    sColor.c1 = 0;
+		    sColor.c2 = 0;
+		    sColor.c3 = 0;
+		    sColor.c4 = 0;
+
+		    G_debug(3,
+			"G_get_color: N, rcount %d, nRed %d, nGreen %d, nBlue %d",
+			rcount, nRed, nGreen, nBlue);
+		    GDALSetColorEntry(hCT, iColor, &sColor);
+		}
+	    }
+	}
+
+	if (rcount > 0) {
+	    /* Create metadata entries for color table rules */
+	    sprintf(value, "%d", rcount);
+	    GDALSetMetadataItem(hBand, "COLOR_TABLE_RULES_COUNT", value,
+				NULL);
+	}
+
+	/* Add the rules in reverse order */
+	for (i = rcount - 1; i >= 0; i--) {
+	    DCELL val1, val2;
+	    unsigned char r1, g1, b1, r2, g2, b2;
+
+	    G_get_f_color_rule(&val1, &r1, &g1, &b1, &val2, &r2, &g2, &b2,
+			       &sGrassColors, i);
+
+
+	    sprintf(key, "COLOR_TABLE_RULE_RGB_%d", rcount - i - 1);
+	    sprintf(value, "%e %e %d %d %d %d %d %d", val1, val2, r1, g1, b1,
+		    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);
+
+    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 memory raster */
+    int row, col;
+    int n_nulls = 0;
+
+    if (maptype == FCELL_TYPE) {
+
+	/* Source datatype understandible by GDAL */
+	GDALDataType datatype = GDT_Float32;
+	FCELL fnullval = (FCELL) nodataval;
+
+	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;
+		    if (n_nulls == 0) {
+			GDALSetRasterNoDataValue(hBand, nodataval);
+		    }
+		    n_nulls++;
+		}
+
+	    if (GDALRasterIO
+		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
+		 0, 0) >= CE_Failure) {
+		G_warning(_("Unable to write GDAL raster file"));
+		return -1;
+	    }
+	    G_percent(row + 1, rows, 2);
+	}
+    }
+    else if (maptype == DCELL_TYPE) {
+
+	GDALDataType datatype = GDT_Float64;
+	DCELL dnullval = (DCELL) nodataval;
+
+	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;
+		    if (n_nulls == 0) {
+			GDALSetRasterNoDataValue(hBand, nodataval);
+		    }
+		    n_nulls++;
+		}
+
+	    if (GDALRasterIO
+		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
+		 0, 0) >= CE_Failure) {
+		G_warning(_("Unable to write GDAL raster file"));
+		return -1;
+	    }
+	    G_percent(row + 1, rows, 2);
+	}
+    }
+    else {
+
+	GDALDataType datatype = GDT_Int32;
+	CELL inullval = (CELL) nodataval;
+
+	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;
+		    if (n_nulls == 0) {
+			GDALSetRasterNoDataValue(hBand, nodataval);
+		    }
+		    n_nulls++;
+		}
+
+	    if (GDALRasterIO
+		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
+		 0, 0) >= CE_Failure) {
+		G_warning(_("Unable to write GDAL raster file"));
+		return -1;
+	    }
+	    G_percent(row + 1, rows, 2);
+	}
+    }
+
+    if (n_nulls > 0) {  /* TODO: && nodata_param NOT specified */
+	if (maptype == CELL_TYPE)
+	    G_warning(_("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);
+	else
+	    G_warning(_("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;
+}

Added: grass/branches/develbranch_6/raster/r.out.gdal/local_proto.h
===================================================================
--- grass/branches/develbranch_6/raster/r.out.gdal/local_proto.h	                        (rev 0)
+++ grass/branches/develbranch_6/raster/r.out.gdal/local_proto.h	2008-10-27 03:56:41 UTC (rev 34015)
@@ -0,0 +1,55 @@
+#ifndef __LOCAL_PROTO_H__
+#define __LOCAL_PROTO_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
+
+  Byte                               0              255
+  UInt16                             0           65,535
+  Int16, CInt16                -32,768           32,767
+  UInt32                             0    4,294,967,295
+  Int32, CInt32         -2,147,483,648    2,147,483,647
+  Float32, CFloat32            -3.4E38           3.4E38
+  Float64, CFloat64          -1.79E308         1.79E308
+*/
+
+/* TODO: for data export range checks and nodata value check
+#define TYPE_BYTE_MIN		0
+#define TYPE_BYTE_MAX		255
+#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_INT32_MAX		2147483647
+#define TYPE_FLOAT32_MIN	-3.4E38f
+#define TYPE_FLOAT32_MAX	3.4E38f
+#define TYPE_FLOAT64_MIN	-1.79E308f
+#define TYPE_FLOAT64_MAX	1.79E308f
+
+ 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
+*/
+
+
+/* export_band.c */
+int export_band(GDALDatasetH, int, const char *, const char *,
+		struct Cell_head *, RASTER_MAP_TYPE, double,
+		const char *, int);
+
+#endif /* __LOCAL_PROTO_H__ */

Modified: grass/branches/develbranch_6/raster/r.out.gdal/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.out.gdal/main.c	2008-10-26 23:52:01 UTC (rev 34014)
+++ grass/branches/develbranch_6/raster/r.out.gdal/main.c	2008-10-27 03:56:41 UTC (rev 34015)
@@ -3,11 +3,11 @@
 *
 * MODULE:       r.out.gdal
 * AUTHOR(S):    Vytautas Vebra <olivership at gmail.com>
-* PURPOSE:      exports GRASS raster to many GDAL suported formats;
+* PURPOSE:      Exports GRASS raster to GDAL suported formats;
 *               based on GDAL library.
-*               Replace r.out.gdal.sh script based on gdal_translate
-*               executable.
-* COPYRIGHT:    (C) 2006-2007 by the GRASS Development Team
+*               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
 *
 *               This program is free software under the GNU General Public
 *   	    	License (>=v2). Read the file COPYING that comes with GRASS
@@ -31,8 +31,8 @@
 
 #include "cpl_string.h"
 #include "gdal.h"
+#include "local_proto.h"
 
-#define GRASS_MAX_COLORS 100000	/* what is the right value ? */
 
 void supported_formats(char **formats)
 {
@@ -81,275 +81,7 @@
     return;
 }
 
-int import_band(GDALDatasetH hMEMDS, int band, const char *name,
-		const char *mapset, struct Cell_head *cellhead,
-		RASTER_MAP_TYPE maptype, double nodataval,
-		const char *nodatakey)
-{
 
-    struct Colors sGrassColors;
-    GDALColorTableH hCT;
-    int iColor;
-    int bHaveMinMax;
-    double dfCellMin;
-    double dfCellMax;
-    struct FPRange sRange;
-    int fd;
-    int cols = cellhead->cols;
-    int rows = cellhead->rows;
-
-    /* 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 raster band  */
-    GDALRasterBandH hBand = GDALGetRasterBand(hMEMDS, band);
-
-    if (hBand == NULL) {
-	G_warning(_("Unable to get raster band"));
-	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);
-    }
-
-    /* suppress useless warnings */
-    CPLPushErrorHandler(CPLQuietErrorHandler);
-    GDALSetRasterColorInterpretation(hBand, GPI_RGB);
-    CPLPopErrorHandler();
-
-    if (G_read_colors(name, mapset, &sGrassColors) >= 0) {
-	int maxcolor, i;
-	CELL min, max;
-	char key[200], value[200];
-	int rcount;
-
-	G_get_color_range(&min, &max, &sGrassColors);
-	if (bHaveMinMax) {
-	    if (max < dfCellMax) {
-		maxcolor = max;
-	    }
-	    else {
-		maxcolor = (int)ceil(dfCellMax);
-	    }
-	    if (maxcolor > GRASS_MAX_COLORS) {
-		maxcolor = GRASS_MAX_COLORS;
-		G_warning("Too many values, color table cut to %d entries",
-			  maxcolor);
-	    }
-	}
-	else {
-	    if (max < GRASS_MAX_COLORS) {
-		maxcolor = max;
-	    }
-	    else {
-		maxcolor = GRASS_MAX_COLORS;
-		G_warning("Too many values, color table set to %d entries",
-			  maxcolor);
-	    }
-	}
-
-	rcount = G_colors_count(&sGrassColors);
-
-	G_debug(3, "dfCellMin: %f, dfCellMax: %f, maxcolor: %d", dfCellMin,
-		dfCellMax, maxcolor);
-
-	hCT = GDALCreateColorTable(GPI_RGB);
-	for (iColor = 0; iColor <= maxcolor; iColor++) {
-	    int nRed, nGreen, nBlue;
-	    GDALColorEntry sColor;
-
-	    if (G_get_color(iColor, &nRed, &nGreen, &nBlue, &sGrassColors)) {
-		sColor.c1 = nRed;
-		sColor.c2 = nGreen;
-		sColor.c3 = nBlue;
-		sColor.c4 = 255;
-
-		G_debug(3,
-			"G_get_color: Y, rcount %d, nRed %d, nGreen %d, nBlue %d",
-			rcount, nRed, nGreen, nBlue);
-		GDALSetColorEntry(hCT, iColor, &sColor);
-	    }
-	    else {
-		sColor.c1 = 0;
-		sColor.c2 = 0;
-		sColor.c3 = 0;
-		sColor.c4 = 0;
-
-		G_debug(3,
-			"G_get_color: N, rcount %d, nRed %d, nGreen %d, nBlue %d",
-			rcount, nRed, nGreen, nBlue);
-		GDALSetColorEntry(hCT, iColor, &sColor);
-	    }
-	}
-
-	if (rcount > 0) {
-	    /* Create metadata entries for color table rules */
-	    sprintf(value, "%d", rcount);
-	    GDALSetMetadataItem(hBand, "COLOR_TABLE_RULES_COUNT", value,
-				NULL);
-	}
-
-	/* Add the rules in reverse order */
-	for (i = rcount - 1; i >= 0; i--) {
-	    DCELL val1, val2;
-	    unsigned char r1, g1, b1, r2, g2, b2;
-
-	    G_get_f_color_rule(&val1, &r1, &g1, &b1, &val2, &r2, &g2, &b2,
-			       &sGrassColors, i);
-
-
-	    sprintf(key, "COLOR_TABLE_RULE_RGB_%d", rcount - i - 1);
-	    sprintf(value, "%e %e %d %d %d %d %d %d", val1, val2, r1, g1, b1,
-		    r2, g2, b2);
-	    GDALSetMetadataItem(hBand, key, value, NULL);
-	}
-	GDALSetRasterColorTable(hBand, hCT);
-    }
-    else {
-	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);
-
-    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 memory raster */
-    int row, col;
-    int n_nulls = 0;
-
-    if (maptype == FCELL_TYPE) {
-
-	/* Source datatype understandible by GDAL */
-	GDALDataType datatype = GDT_Float32;
-	FCELL fnullval = (FCELL) nodataval;
-
-	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;
-		    if (n_nulls == 0) {
-			GDALSetRasterNoDataValue(hBand, nodataval);
-		    }
-		    n_nulls++;
-		}
-
-	    if (GDALRasterIO
-		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
-		 0, 0) >= CE_Failure) {
-		G_warning(_("Unable to write GDAL raster file"));
-		return -1;
-	    }
-	    G_percent(row + 1, rows, 2);
-	}
-    }
-    else if (maptype == DCELL_TYPE) {
-
-	GDALDataType datatype = GDT_Float64;
-	DCELL dnullval = (DCELL) nodataval;
-
-	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;
-		    if (n_nulls == 0) {
-			GDALSetRasterNoDataValue(hBand, nodataval);
-		    }
-		    n_nulls++;
-		}
-
-	    if (GDALRasterIO
-		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
-		 0, 0) >= CE_Failure) {
-		G_warning(_("Unable to write GDAL raster file"));
-		return -1;
-	    }
-	    G_percent(row + 1, rows, 2);
-	}
-    }
-    else {
-
-	GDALDataType datatype = GDT_Int32;
-	CELL inullval = (CELL) nodataval;
-
-	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;
-		    if (n_nulls == 0) {
-			GDALSetRasterNoDataValue(hBand, nodataval);
-		    }
-		    n_nulls++;
-		}
-
-	    if (GDALRasterIO
-		(hBand, GF_Write, 0, row, cols, 1, bufer, cols, 1, datatype,
-		 0, 0) >= CE_Failure) {
-		G_warning(_("Unable to write GDAL raster file"));
-		return -1;
-	    }
-	    G_percent(row + 1, rows, 2);
-	}
-    }
-
-    if (n_nulls > 0) {  /* TODO: && nodata_param NOT specified */
-	if (maptype == CELL_TYPE)
-	    G_warning(_("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);
-	else
-	    G_warning(_("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;
-}
-
 static void AttachMetadata(GDALDatasetH hDS, char **papszMetadataOptions)
 /* function taken from gdal_translate */
 {
@@ -373,7 +105,7 @@
 {
 
     struct GModule *module;
-    struct Flag *flag_l;
+    struct Flag *flag_l, *flag_c;
     struct Option *input, *format, *type, *output, *createopt, *metaopt,
 	*nodataopt;
 
@@ -390,7 +122,7 @@
 
     module = G_define_module();
     module->description =
-	_("Exports GRASS raster map into GDAL supported formats.");
+	_("Exports GRASS raster maps into GDAL supported formats.");
     module->keywords = _("raster, export");
 
     flag_l = G_define_flag();
@@ -398,6 +130,11 @@
     flag_l->description = _("List supported output formats");
     flag_l->guisection = _("Print");
 
+    flag_c = G_define_flag();
+    flag_c->key = 'c';
+    flag_c->label = _("Do not export long colortable");
+    flag_c->description = _("Only applicable to Byte or UInt16 data types.");
+
     input = G_define_standard_option(G_OPT_R_INPUT);
     input->required = NO;
     input->description = _("Name of raster map (or group) to export");
@@ -430,7 +167,7 @@
     type->type = TYPE_STRING;
     type->description = _("File type");
     type->options =
-	"Byte,Int16,UInt16,UInt32,Int32,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64";
+	"Byte,Int16,UInt16,Int32,UInt32,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64";
     type->required = NO;
 
     output = G_define_standard_option(G_OPT_R_OUTPUT);
@@ -442,8 +179,8 @@
     createopt = G_define_option();
     createopt->key = "createopt";
     createopt->type = TYPE_STRING;
-    createopt->description = _("Creation option to the output format driver. "
-			       "Multiple options may be listed.");
+    createopt->description =
+	_("Creation option(s) to pass to the output format driver");
     createopt->multiple = YES;
     createopt->required = NO;
 
@@ -463,9 +200,11 @@
     nodataopt->multiple = NO;
     nodataopt->required = NO;
 
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
+
 #ifndef __ALLOW_DYNAMIC_OPTIONS__
     /* Init GDAL */
     GDALAllRegister();
@@ -476,9 +215,8 @@
 	exit(EXIT_SUCCESS);
     }
 
-    if (!input->answer) {
+    if (!input->answer)
 	G_fatal_error(_("Required parameter <%s> not set"), input->key);
-    }
 
     /* Try to open input GRASS raster.. */
     mapset = G_find_cell2(input->answer, "");
@@ -633,7 +371,7 @@
 		    nodataval = (double)(short)0x8000;
 		}
 		else {
-		    datatype = GDT_Int32;	/* need to furthermore fine tune this? */
+		    datatype = GDT_Int32;  /* need to fine tune this more? */
 		    nodataval = (double)(int)0x80000000;
 		}
 	    }
@@ -645,6 +383,11 @@
 	nodataval = atof(nodataopt->answer);
     }
 
+/* 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" :
@@ -722,9 +465,9 @@
 						     ref.file[band].mapset),
 			      band + 1);
 	}
-	if (import_band
+	if (export_band
 	    (hCurrDS, band + 1, ref.file[band].name, ref.file[band].mapset,
-	     &cellhead, maptype, nodataval, nodataopt->key) < 0)
+	     &cellhead, maptype, nodataval, nodataopt->key, flag_c->answer) < 0)
 	    G_warning(_("Unable to export raster map <%s>"),
 		      ref.file[band].name);
     }
@@ -740,12 +483,12 @@
     }
 
     GDALClose(hDstDS);
+
     if (hMEMDS)
 	GDALClose(hMEMDS);
 
     CSLDestroy(papszOptions);
 
     G_done_msg(" ");
-
     exit(EXIT_SUCCESS);
 }



More information about the grass-commit mailing list