Add write support to GDAL-link, r.external.out

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

include/gisdefs.h

@@ -549,7 +549,9 @@ void G_fpreclass_perform_id(const struct FPReclass *, const CELL *, DCELL *,
 /* gdal.c */
 void G_init_gdal(void);
 struct GDAL_link *G_get_gdal_link(const char *, const char *);
+struct GDAL_link *G_create_gdal_link(const char *, RASTER_MAP_TYPE);
 void G_close_gdal_link(struct GDAL_link *);
+int G_close_gdal_write_link(struct GDAL_link *);
 
 /* geodesic.c */
 int G_begin_geodesic_equation(double, double, double, double);
@@ -932,7 +934,6 @@ int G_put_cellhd(const char *, struct Cell_head *);
 
 /* put_row.c */
 int G_put_map_row(int, const CELL *);
-int G__put_null_value_row(int, const char *);
 int G_put_raster_row(int, const void *, RASTER_MAP_TYPE);
 int G_put_c_raster_row(int, const CELL *);
 int G_put_f_raster_row(int, const FCELL *);

lib/gis/closecell.c

@@ -158,17 +158,178 @@ static int close_old(int fd)
     return 1;
 }
 
-static int close_new(int fd, int ok)
+static void write_support_files(int fd)
 {
     struct fileinfo *fcb = &G__.fileinfo[fd];
-    int stat;
     struct Categories cats;
     struct History hist;
-    char path[GPATH_MAX];
     CELL cell_min, cell_max;
-    int row, i, open_mode;
+    char path[GPATH_MAX];
+
+    /* remove color table */
+    G_remove_colors(fcb->name, "");
+
+    /* create a history file */
+    G_short_history(fcb->name, "raster", &hist);
+    G_write_history(fcb->name, &hist);
+
+    /* write the range */
+    if (fcb->map_type == CELL_TYPE) {
+	G_write_range(fcb->name, &fcb->range);
+	G__remove_fp_range(fcb->name);
+    }
+    /*NOTE: int range for floating point maps is not written out */
+    else {			/* if(fcb->map_type != CELL_TYPE) */
+
+	G_write_fp_range(fcb->name, &fcb->fp_range);
+	G_construct_default_range(&fcb->range);
+	/* this range will be used to add default rule to quant structure */
+    }
+
+    if (fcb->map_type != CELL_TYPE)
+	fcb->cellhd.format = -1;
+    else			/* CELL map */
+	fcb->cellhd.format = fcb->nbytes - 1;
+
+    /* write header file */
+    G_put_cellhd(fcb->name, &fcb->cellhd);
+
+    /* if map is floating point write the quant rules, otherwise remove f_quant */
+    if (fcb->map_type != CELL_TYPE) {
+	/* DEFAULT RANGE QUANT
+	   G_get_fp_range_min_max(&fcb->fp_range, &dcell_min, &dcell_max);
+	   if(!G_is_d_null_value(&dcell_min) && !G_is_d_null_value(&dcell_max))
+	   {
+	   G_get_range_min_max(&fcb->range, &cell_min, &cell_max);
+	   G_quant_add_rule(&fcb->quant, dcell_min, dcell_max, 
+	   cell_min, cell_max);
+	   }
+	*/
+	G_quant_round(&fcb->quant);
+	if (G_write_quant(fcb->name, fcb->mapset, &fcb->quant) < 0)
+	    G_warning(_("unable to write quant file!"));
+    }
+    else {
+	/* remove cell_misc/name/f_quant */
+	G__file_name_misc(path, "cell_misc", QUANT_FILE, fcb->name,
+			  fcb->mapset);
+	remove(path);
+    }
+
+    /* create empty cats file */
+    G_get_range_min_max(&fcb->range, &cell_min, &cell_max);
+    if (G_is_c_null_value(&cell_max))
+	cell_max = 0;
+    G_init_cats(cell_max, (char *)NULL, &cats);
+    G_write_cats(fcb->name, &cats);
+    G_free_cats(&cats);
+
+    /* write the histogram */
+    /* only works for integer maps */
+    if ((fcb->map_type == CELL_TYPE)
+	&& (fcb->want_histogram)) {
+	G_write_histogram_cs(fcb->name, &fcb->statf);
+	G_free_cell_stats(&fcb->statf);
+    }
+    else {
+	G_remove_histogram(fcb->name);
+    }
+}
+
+static int close_new_gdal(int fd, int ok)
+{
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+    char path[GPATH_MAX];
+    int stat = 1;
+
+    if (ok) {
+	int cell_fd;
+
+	G_debug(1, "close %s GDAL", fcb->name);
+
+	if (fcb->cur_row < fcb->cellhd.rows) {
+	    int row;
+
+	    G_zero_raster_buf(fcb->data, fcb->map_type);
+	    for (row = fcb->cur_row; row < fcb->cellhd.rows; row++)
+		G_put_raster_row(fd, fcb->data, fcb->map_type);
+	    G_free(fcb->data);
+	    fcb->data = NULL;
+	}
+
+	/* create path : full null file name */
+	G__make_mapset_element_misc("cell_misc", fcb->name);
+	G__file_name_misc(path, "cell_misc", NULL_FILE, fcb->name,
+			  G_mapset());
+	remove(path);
+
+	/* write 0-length cell file */
+	G__make_mapset_element("cell");
+	G__file_name(path, "cell", fcb->name, fcb->mapset);
+	cell_fd = creat(path, 0666);
+	close(cell_fd);
+
+	if (fcb->map_type != CELL_TYPE) {	/* floating point map */
+	    if (write_fp_format(fd) != 0) {
+		G_warning(_("Error writing floating point format file for map %s"),
+			  fcb->name);
+		stat = -1;
+	    }
+
+	    /* write 0-length fcell file */
+	    G__make_mapset_element("fcell");
+	    G__file_name(path, "fcell", fcb->name, fcb->mapset);
+	    cell_fd = creat(path, 0666);
+	    close(cell_fd);
+	}
+	else {
+	    /* remove fcell/name file */
+	    G__file_name(path, "fcell", fcb->name, fcb->mapset);
+	    remove(path);
+	    /* remove cell_misc/name/f_format */
+	    G__file_name_misc(path, "cell_misc", FORMAT_FILE, fcb->name,
+			      fcb->mapset);
+	    remove(path);
+	}
+
+	if (G_close_gdal_write_link(fcb->gdal) < 0)
+	    stat = -1;
+    }
+    else {
+	remove(fcb->gdal->filename);
+	G_close_gdal_link(fcb->gdal);
+    }
+
+    /* NOW CLOSE THE FILE DESCRIPTOR */
+    close(fd);
+    fcb->open_mode = -1;
+
+    if (fcb->data != NULL)
+	G_free(fcb->data);
+
+    if (ok)
+	write_support_files(fd);
+
+    G_free(fcb->name);
+    G_free(fcb->mapset);
+
+    if (fcb->map_type != CELL_TYPE)
+	G_quant_free(&fcb->quant);
+
+    return stat;
+}
+
+static int close_new(int fd, int ok)
+{
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+    int stat;
+    char path[GPATH_MAX];
+    int row, i;
     const char *CELL_DIR;
 
+    if (fcb->gdal)
+	return close_new_gdal(fd, ok);
+
     if (ok) {
 	switch (fcb->open_mode) {
 	case OPEN_NEW_COMPRESSED:
@@ -272,8 +433,6 @@ static int close_new(int fd, int ok)
     /* NOW CLOSE THE FILE DESCRIPTOR */
 
     close(fd);
-    /* remember open_mode */
-    open_mode = fcb->open_mode;
     fcb->open_mode = -1;
 
     if (fcb->data != NULL)
@@ -307,76 +466,8 @@ static int close_new(int fd, int ok)
 	G_free(fcb->temp_name);
     }
 
-    if (ok) {
-	/* remove color table */
-	G_remove_colors(fcb->name, "");
-
-	/* create a history file */
-	G_short_history(fcb->name, "raster", &hist);
-	G_write_history(fcb->name, &hist);
-
-	/* write the range */
-	if (fcb->map_type == CELL_TYPE) {
-	    G_write_range(fcb->name, &fcb->range);
-	    G__remove_fp_range(fcb->name);
-	}
-	/*NOTE: int range for floating point maps is not written out */
-	else {			/* if(fcb->map_type != CELL_TYPE) */
-
-	    G_write_fp_range(fcb->name, &fcb->fp_range);
-	    G_construct_default_range(&fcb->range);
-	    /* this range will be used to add default rule to quant structure */
-	}
-
-	if (fcb->map_type != CELL_TYPE)
-	    fcb->cellhd.format = -1;
-	else			/* CELL map */
-	    fcb->cellhd.format = fcb->nbytes - 1;
-
-	/* write header file */
-	G_put_cellhd(fcb->name, &fcb->cellhd);
-
-	/* if map is floating point write the quant rules, otherwise remove f_quant */
-	if (fcb->map_type != CELL_TYPE) {
-	    /* DEFAULT RANGE QUANT
-	       G_get_fp_range_min_max(&fcb->fp_range, &dcell_min, &dcell_max);
-	       if(!G_is_d_null_value(&dcell_min) && !G_is_d_null_value(&dcell_max))
-	       {
-	       G_get_range_min_max(&fcb->range, &cell_min, &cell_max);
-	       G_quant_add_rule(&fcb->quant, dcell_min, dcell_max, 
-	       cell_min, cell_max);
-	       }
-	     */
-	    G_quant_round(&fcb->quant);
-	    if (G_write_quant(fcb->name, fcb->mapset, &fcb->quant) < 0)
-		G_warning(_("unable to write quant file!"));
-	}
-	else {
-	    /* remove cell_misc/name/f_quant */
-	    G__file_name_misc(path, "cell_misc", QUANT_FILE, fcb->name,
-			      fcb->mapset);
-	    remove(path);
-	}
-
-	/* create empty cats file */
-	G_get_range_min_max(&fcb->range, &cell_min, &cell_max);
-	if (G_is_c_null_value(&cell_max))
-	    cell_max = 0;
-	G_init_cats(cell_max, (char *)NULL, &cats);
-	G_write_cats(fcb->name, &cats);
-	G_free_cats(&cats);
-
-	/* write the histogram */
-	/* only works for integer maps */
-	if ((fcb->map_type == CELL_TYPE)
-	    && (fcb->want_histogram)) {
-	    G_write_histogram_cs(fcb->name, &fcb->statf);
-	    G_free_cell_stats(&fcb->statf);
-	}
-	else {
-	    G_remove_histogram(fcb->name);
-	}
-    }				/* OK */
+    if (ok)
+	write_support_files(fd);
 
     G_free(fcb->name);
     G_free(fcb->mapset);

lib/gis/gdal.c

@@ -1,6 +1,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <unistd.h>
 #include <grass/config.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
@@ -31,6 +32,15 @@ static CPLErr CPL_STDCALL (*pGDALRasterIO)(
     int nDSXOff, int nDSYOff, int nDSXSize, int nDSYSize,
     void * pBuffer, int nBXSize, int nBYSize,GDALDataType eBDataType,
     int nPixelSpace, int nLineSpace);
+static GDALDriverH CPL_STDCALL (*pGDALGetDriverByName)(const char *);
+static const char * CPL_STDCALL (*pGDALGetMetadataItem)(GDALMajorObjectH, const char *, const char *);
+static GDALDatasetH CPL_STDCALL (*pGDALCreate)(GDALDriverH hDriver, const char *, int, int, int, GDALDataType, char **);
+static GDALDatasetH CPL_STDCALL (*pGDALCreateCopy)(GDALDriverH, const char *, GDALDatasetH, int, char **, GDALProgressFunc, void *);
+static CPLErr CPL_STDCALL (*pGDALSetRasterNoDataValue)(GDALRasterBandH, double);
+static CPLErr CPL_STDCALL (*pGDALSetGeoTransform)(GDALDatasetH, double *);
+static CPLErr CPL_STDCALL (*pGDALSetProjection)(GDALDatasetH, const char *);
+static const char * CPL_STDCALL (*pGDALGetDriverShortName)(GDALDriverH);
+static GDALDriverH CPL_STDCALL (*pGDALGetDatasetDriver)(GDALDatasetH);
 
 #if GDAL_DYNAMIC
 # if defined(__unix) && !defined(__unix__)
@@ -98,22 +108,40 @@ static void init_gdal(void)
 {
     load_library();
 
-    pGDALAllRegister   = get_symbol("GDALAllRegister");
-    pGDALOpen          = get_symbol("GDALOpen");
-    pGDALClose         = get_symbol("GDALClose");
-    pGDALGetRasterBand = get_symbol("GDALGetRasterBand");
-    pGDALRasterIO      = get_symbol("GDALRasterIO");
+    pGDALAllRegister          = get_symbol("GDALAllRegister");
+    pGDALOpen                 = get_symbol("GDALOpen");
+    pGDALClose                = get_symbol("GDALClose");
+    pGDALGetRasterBand        = get_symbol("GDALGetRasterBand");
+    pGDALRasterIO             = get_symbol("GDALRasterIO");
+    pGDALGetDriverByName      = get_symbol("GDALGetDriverByName");
+    pGDALGetMetadataItem      = get_symbol("GDALGetMetadataItem");
+    pGDALCreate               = get_symbol("GDALCreate");
+    pGDALCreateCopy           = get_symbol("GDALCreateCopy");
+    pGDALSetRasterNoDataValue = get_symbol("GDALSetRasterNoDataValue");
+    pGDALSetGeoTransform      = get_symbol("GDALSetGeoTransform");
+    pGDALSetProjection        = get_symbol("GDALSetProjection");
+    pGDALGetDriverShortName   = get_symbol("GDALGetDriverShortName");
+    pGDALGetDatasetDriver     = get_symbol("GDALGetDatasetDriver");
 }
 
 #else /* GDAL_DYNAMIC */
 
 static void init_gdal(void)
 {
-    pGDALAllRegister   = &GDALAllRegister;
-    pGDALOpen          = &GDALOpen;
-    pGDALClose         = &GDALClose;
-    pGDALGetRasterBand = &GDALGetRasterBand;
-    pGDALRasterIO      = &GDALRasterIO;
+    pGDALAllRegister          = &GDALAllRegister;
+    pGDALOpen                 = &GDALOpen;
+    pGDALClose                = &GDALClose;
+    pGDALGetRasterBand        = &GDALGetRasterBand;
+    pGDALRasterIO             = &GDALRasterIO;
+    pGDALGetDriverByName      = &GDALGetDriverByName;
+    pGDALGetMetadataItem      = &GDALGetMetadataItem;
+    pGDALCreate               = &GDALCreate;
+    pGDALCreateCopy           = &GDALCreateCopy;
+    pGDALSetRasterNoDataValue = &GDALSetRasterNoDataValue;
+    pGDALSetGeoTransform      = &GDALSetGeoTransform;
+    pGDALSetProjection        = &GDALSetProjection;
+    pGDALGetDriverShortName   = &GDALGetDriverShortName;
+    pGDALGetDatasetDriver     = &GDALGetDatasetDriver;
 }
 
 #endif /* GDAL_DYNAMIC */
@@ -246,6 +274,211 @@ struct GDAL_link *G_get_gdal_link(const char *name, const char *mapset)
     return gdal;
 }
 
+struct GDAL_Options {
+    const char *dir;
+    const char *ext;
+    const char *format;
+    char **options;
+};
+
+static struct state {
+    int initialized;
+    struct GDAL_Options opts;
+    struct Key_Value *projinfo, *projunits;
+    char *srswkt;
+} state;
+
+static struct state *st = &state;
+
+static void read_gdal_options(void)
+{
+    FILE *fp;
+    struct Key_Value *key_val;
+    const char *p;
+
+    fp = G_fopen_old("", "GDAL", G_mapset());
+    if (!fp)
+	G_fatal_error(_("Unable to open GDAL file"));
+    key_val = G_fread_key_value(fp);
+    fclose(fp);
+
+    p = G_find_key_value("directory", key_val);
+    if (!p)
+	p = "gdal";
+    if (*p == '/')
+	st->opts.dir = G_store(p);
+    else;
+    {
+	char path[GPATH_MAX];
+	G__file_name(path, p, "", G_mapset());
+	st->opts.dir = G_store(path);
+	if (access(path, 0) != 0)
+	    G__make_mapset_element(p);
+    }
+
+    p = G_find_key_value("extension", key_val);
+    st->opts.ext = G_store(p ? p : "");
+
+    p = G_find_key_value("format", key_val);
+    st->opts.format = G_store(p ? p : "GTiff");
+
+    p = G_find_key_value("options", key_val);
+    st->opts.options = p ? G_tokenize(p, ",") : NULL;
+
+    G_free_key_value(key_val);
+}
+
+struct GDAL_link *G_create_gdal_link(const char *name, RASTER_MAP_TYPE map_type)
+{
+#ifdef GDAL_LINK
+    char path[GPATH_MAX];
+    GDALDriverH driver;
+    double transform[6];
+    struct GDAL_link *gdal;
+    FILE *fp;
+    struct Key_Value *key_val;
+    char buf[32];
+
+
+    G_init_gdal();
+
+    if (!G_is_initialized(&st->initialized)) {
+	read_gdal_options();
+	st->projinfo = G_get_projinfo();
+	st->projunits = G_get_projunits();
+#if 0
+	/* We cannot use GPJ_grass_to_wkt() here because that would create a
+	   circular dependency between libgis and libgproj */
+	if (st->projinfo && st->projunits)
+	    st->srswkt = GPJ_grass_to_wkt(st->projinfo, st->projunits);
+#endif
+	G_initialize_done(&st->initialized);
+    }
+
+    gdal = G_calloc(1, sizeof(struct GDAL_link));
+
+    sprintf(path, "%s/%s%s", st->opts.dir, name, st->opts.ext);
+    gdal->filename = G_store(path);
+    gdal->band_num = 1;
+    gdal->hflip = 0;
+    gdal->vflip = 0;
+
+    switch (map_type) {
+    case CELL_TYPE:
+	switch (G__.nbytes) {
+	case 1:
+	    gdal->type = GDT_Byte;
+	    gdal->null_val = (DCELL) 0xFF;
+	    break;
+	case 2:
+	    gdal->type = GDT_UInt16;
+	    gdal->null_val = (DCELL) 0xFFFF;
+	    break;
+	case 3:
+	case 4:
+	    gdal->type = GDT_Int32;
+	    gdal->null_val = (DCELL) 0x80000000U;
+	    break;
+	}
+	break;
+    case FCELL_TYPE:
+	gdal->type = GDT_Float32;
+	G_set_d_null_value(&gdal->null_val, 1);
+	break;
+    case DCELL_TYPE:
+	gdal->type = GDT_Float64;
+	G_set_d_null_value(&gdal->null_val, 1);
+	break;
+    default:
+	G_fatal_error(_("Invalid map type <%d>"), map_type);
+	break;
+    }
+
+    driver = (*pGDALGetDriverByName)(st->opts.format);
+    if (!driver)
+	G_fatal_error(_("Unable to get <%s> driver"), st->opts.format);
+
+    /* Does driver support GDALCreate ? */
+    if ((*pGDALGetMetadataItem)(driver, GDAL_DCAP_CREATE, NULL))
+    {
+	gdal->data = (*pGDALCreate)(driver, gdal->filename, G__.window.cols, G__.window.rows,
+				1, gdal->type, st->opts.options);
+	if (!gdal->data)
+	    G_fatal_error(_("Unable to create <%s> dataset using <%s> driver"),
+			  name, st->opts.format);
+    }
+    /* If not - create MEM driver for intermediate dataset. 
+     * Check if raster can be created at all (with GDALCreateCopy) */
+    else if ((*pGDALGetMetadataItem)(driver, GDAL_DCAP_CREATECOPY, NULL)) {
+	GDALDriverH mem_driver;
+
+	G_message(_("Driver <%s> does not support direct writing. "
+		    "Using MEM driver for intermediate dataset."),
+		  st->opts.format);
+
+	mem_driver = (*pGDALGetDriverByName)("MEM");
+	if (!mem_driver)
+	    G_fatal_error(_("Unable to get in-memory raster driver"));
+
+	gdal->data = (*pGDALCreate)(mem_driver, "", G__.window.cols, G__.window.rows,
+				    1, gdal->type, st->opts.options);
+	if (!gdal->data)
+	    G_fatal_error(_("Unable to create <%s> dataset using memory driver"),
+			  name);
+    }
+    else
+	G_fatal_error(_("Driver <%s> does not support creating rasters"),
+		      st->opts.format);
+
+    gdal->band = (*pGDALGetRasterBand)(gdal->data, gdal->band_num);
+
+    (*pGDALSetRasterNoDataValue)(gdal->band, gdal->null_val);
+
+    /* Set Geo Transform  */
+    transform[0] = G__.window.west;
+    transform[1] = G__.window.ew_res;
+    transform[2] = 0.0;
+    transform[3] = G__.window.north;
+    transform[4] = 0.0;
+    transform[5] = -G__.window.ns_res;
+
+    if ((*pGDALSetGeoTransform)(gdal->data, transform) >= CE_Failure)
+	G_warning(_("Unable to set geo transform"));
+
+    if (st->srswkt)
+	if ((*pGDALSetProjection)(gdal->data, st->srswkt) == CE_Failure)
+	    G_warning(_("Unable to set projection"));
+
+    fp = G_fopen_new_misc("cell_misc", "gdal", name);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell_misc/%s/gdal file"), name);
+
+    key_val = G_create_key_value();
+
+    G_set_key_value("file", gdal->filename, key_val);
+
+    sprintf(buf, "%d", gdal->band_num);
+    G_set_key_value("band", buf, key_val);
+
+    sprintf(buf, "%.22g", gdal->null_val);
+    G_set_key_value("null", buf, key_val);
+
+    sprintf(buf, "%d", gdal->type);
+    G_set_key_value("type", buf, key_val);
+
+    if (G_fwrite_key_value(fp, key_val) < 0)
+	G_fatal_error(_("Error writing cell_misc/%s/gdal file"), name);
+
+    G_free_key_value(key_val);
+
+    fclose(fp);
+
+    return gdal;
+#else
+    return NULL;
+#endif
+}
+
 void G_close_gdal_link(struct GDAL_link *gdal)
 {
 #ifdef GDAL_LINK
@@ -255,6 +488,33 @@ void G_close_gdal_link(struct GDAL_link *gdal)
     G_free(gdal);
 }
 
+int G_close_gdal_write_link(struct GDAL_link *gdal)
+{
+    int stat = 1;
+#ifdef GDAL_LINK
+    GDALDriverH src_drv = (*pGDALGetDatasetDriver)(gdal->data);
+
+    if (G_strcasecmp((*pGDALGetDriverShortName)(src_drv), "MEM") == 0) {
+	GDALDriverH dst_drv = (*pGDALGetDriverByName)(st->opts.format);
+	GDALDatasetH dst = (*pGDALCreateCopy)(dst_drv, gdal->filename, gdal->data, FALSE,
+					      st->opts.options, NULL, NULL);
+	if (!dst) {
+	    G_warning(_("Unable to create output file <%s> using driver <%s>"),
+		      gdal->filename, st->opts.format);
+	    stat = -1;
+	}
+	(*pGDALClose)(dst);
+    }
+
+    (*pGDALClose)(gdal->data);
+
+#endif
+    G_free(gdal->filename);
+    G_free(gdal);
+
+    return stat;
+}
+
 #ifdef GDAL_LINK
 CPLErr G_gdal_raster_IO(
     GDALRasterBandH band, GDALRWFlag rw_flag,

lib/gis/opencell.c

@@ -487,6 +487,70 @@ int G_open_fp_cell_new_uncompressed(const char *name)
     return G__open_raster_new(name, OPEN_NEW_UNCOMPRESSED, G__.fp_type);
 }
 
+#ifdef HAVE_GDAL
+static int G__open_raster_new_gdal(char *map, char *mapset, RASTER_MAP_TYPE map_type)
+{
+    int fd;
+    struct fileinfo *fcb;
+    int i;
+
+    /* dummy descriptor to reserve the fileinfo slot */
+    fd = open("/dev/null", O_RDONLY);
+    if (fd < 0)
+	return -1;
+
+    fcb = new_fileinfo(fd);
+
+    /* mark closed */
+    fcb->map_type = map_type;
+    fcb->open_mode = -1;
+
+    fcb->gdal = G_create_gdal_link(map, map_type);
+    if (!fcb->gdal)
+	return -1;
+
+    fcb->cellhd = G__.window;
+    fcb->cellhd.compressed = 0;
+    fcb->nbytes = G_raster_size(fcb->map_type);
+    /* for writing fcb->data is allocated to be G__.window.cols * 
+       sizeof(CELL or DCELL or FCELL)  */
+    fcb->data = G_calloc(G__.window.cols, fcb->nbytes);
+
+    fcb->name = map;
+    fcb->mapset = mapset;
+    fcb->cur_row = 0;
+
+    fcb->row_ptr = NULL;
+    fcb->temp_name = NULL;
+    fcb->null_temp_name = NULL;
+    fcb->null_cur_row = 0;
+    fcb->min_null_row = 0;
+    for (i = 0; i < NULL_ROWS_INMEM; i++)
+	fcb->NULL_ROWS[i] = NULL;
+
+    if (fcb->map_type != CELL_TYPE)
+	G_quant_init(&(fcb->quant));
+
+    /* init cell stats */
+    /* now works only for int maps */
+    if (fcb->map_type == CELL_TYPE)
+	if ((fcb->want_histogram = G__.want_histogram))
+	    G_init_cell_stats(&fcb->statf);
+
+    /* init range and if map is double/float init d/f_range */
+    G_init_range(&fcb->range);
+
+    if (fcb->map_type != CELL_TYPE)
+	G_init_fp_range(&fcb->fp_range);
+
+    /* mark file as open for write */
+    fcb->open_mode = OPEN_NEW_UNCOMPRESSED;
+    fcb->io_error = 0;
+
+    return fd;
+}
+#endif /* HAVE_GDAL */
+
 static int G__open_raster_new(const char *name, int open_mode,
 			      RASTER_MAP_TYPE map_type)
 {
@@ -532,6 +596,11 @@ static int G__open_raster_new(const char *name, int open_mode,
     /* make sure window is set */
     G__init_window();
 
+#ifdef HAVE_GDAL
+    if (G_find_file2("", "GDAL", G_mapset()))
+	return G__open_raster_new_gdal(map, mapset, map_type);
+#endif
+
     /* open a tempfile name */
     tempname = G_tempfile();
     fd = creat(tempname, 0666);
@@ -553,6 +622,7 @@ static int G__open_raster_new(const char *name, int open_mode,
     /* mark closed */
     fcb->map_type = map_type;
     fcb->open_mode = -1;
+    fcb->gdal = NULL;
 
     /* for writing fcb->data is allocated to be G__.window.cols * 
        sizeof(CELL or DCELL or FCELL)  */

lib/gis/put_row.c

@@ -59,7 +59,7 @@
  *
  ***********************************************************************
  *
- *  G__put_null_value_row(fd, buf)
+ *  put_null_value_row(fd, buf)
  *      int fd                  File descriptor where data is to be written
  *      char *buf               Buffer holding null data
  *
@@ -107,10 +107,13 @@ static int put_raster_row(int fd, const void *buf, RASTER_MAP_TYPE data_type,
 
 /*--------------------------------------------------------------------------*/
 
-int G__put_null_value_row(int fd, const char *buf)
+static int put_null_value_row(int fd, const char *buf)
 {
     struct fileinfo *fcb = &G__.fileinfo[fd];
 
+    if (fcb->gdal)
+	G_fatal_error(_("GDAL output doesn't support writing null rows separately"));
+
     switch (put_null_data(fd, buf, fcb->null_cur_row)) {
     case -1:
 	return -1;
@@ -562,6 +565,58 @@ static int put_data(int fd, char *null_buf, const CELL *cell,
 
 /*--------------------------------------------------------------------------*/
 
+static int put_data_gdal(int fd, const void *rast, int row, int n,
+			 int zeros_r_nulls, RASTER_MAP_TYPE map_type)
+{
+#ifdef HAVE_GDAL
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+    int size = G_raster_size(map_type);
+    DCELL null_val = fcb->gdal->null_val;
+    const void *src;
+    void *work_buf, *dst;
+    GDALDataType datatype;
+    CPLErr err;
+    int i;
+
+    if (row < 0 || row >= fcb->cellhd.rows)
+	return 0;
+
+    if (n <= 0)
+	return 0;
+
+    work_buf = G__alloca(n * size);
+
+    switch (map_type) {
+    case CELL_TYPE:	datatype = GDT_Int32;	break;
+    case FCELL_TYPE:	datatype = GDT_Float32;	break;
+    case DCELL_TYPE:	datatype = GDT_Float64;	break;
+    }
+
+    src = rast;
+    dst = work_buf;
+
+    for (i = 0; i < n; i++) {
+	if (G_is_null_value(src, map_type) || zeros_r_nulls && !*(CELL *)src)
+	    G_set_raster_value_d(dst, null_val, map_type);
+	else
+	    memcpy(dst, src, size);
+	src = G_incr_void_ptr(src, size);
+	dst = G_incr_void_ptr(dst, size);
+    }
+
+    err = G_gdal_raster_IO(fcb->gdal->band, GF_Write, 0, row, n, 1, work_buf,
+			   n, 1, datatype, 0, 0);
+
+    G__freea(work_buf);
+
+    return err == CE_None ? 1 : -1;
+#else
+    return -1;
+#endif
+}
+
+/*--------------------------------------------------------------------------*/
+
 /*--------------------------------------------------------------------------*/
 
 /*--------------------------------------------------------------------------*/
@@ -570,6 +625,11 @@ static int put_raster_data(int fd, char *null_buf, const void *rast,
 			   int row, int n,
 			   int zeros_r_nulls, RASTER_MAP_TYPE map_type)
 {
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+
+    if (fcb->gdal)
+	return put_data_gdal(fd, rast, row, n, zeros_r_nulls, map_type);
+
     return (map_type == CELL_TYPE)
 	? put_data(fd, null_buf, rast, row, n, zeros_r_nulls)
 	: put_fp_data(fd, null_buf, rast, row, n, map_type);
@@ -808,7 +868,10 @@ static int put_raster_row(int fd, const void *buf, RASTER_MAP_TYPE data_type,
     fcb->cur_row++;
 
     /* write the null row for the data row */
-    stat = G__put_null_value_row(fd, null_buf);
+    if (fcb->gdal)
+	stat = 0;
+    else
+	stat = put_null_value_row(fd, null_buf);
     G__freea(null_buf);
     return stat;
 }

raster/Makefile

@@ -21,6 +21,7 @@ SUBDIRS = \
 	r.distance \
 	r.drain \
 	r.external \
+	r.external.out \
 	r.fill.dir \
 	r.flow \
 	r.grow.distance \

raster/r.external.out/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM  = r.external.out
+
+LIBES= $(GISLIB) $(GDALLIBS) $(GMATHLIB)
+DEPENDENCIES = $(GISDEP) $(GMATHDEP)
+EXTRA_INC = $(PROJINC) $(GDALCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+ifneq ($(USE_GDAL),)
+default: cmd
+endif

raster/r.external.out/main.c

@@ -0,0 +1,242 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.external.out
+ *               
+ * AUTHOR(S):    Glynn Clements, based on r.out.gdal
+ *
+ * PURPOSE:      Make GRASS write raster maps utilizing the GDAL library.
+ *
+ * COPYRIGHT:    (C) 2008 by Glynn Clements and 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 <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/imagery.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+
+#include "gdal.h"
+
+static void list_formats(void)
+{
+    /* -------------------------------------------------------------------- */
+    /*      List supported formats and exit.                                */
+    /*         code from GDAL 1.2.5  gcore/gdal_misc.cpp                    */
+    /*         Copyright (c) 1999, Frank Warmerdam                          */
+    /* -------------------------------------------------------------------- */
+    int iDr;
+
+    fprintf(stdout, _("Supported Formats:\n"));
+    for (iDr = 0; iDr < GDALGetDriverCount(); iDr++) {
+	GDALDriverH hDriver = GDALGetDriver(iDr);
+	const char *pszRWFlag;
+
+	if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL))
+	    pszRWFlag = "rw+";
+	else if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, NULL))
+	    pszRWFlag = "rw";
+	else
+	    continue;
+
+	fprintf(stdout, "  %s (%s): %s\n",
+		GDALGetDriverShortName(hDriver),
+		pszRWFlag, GDALGetDriverLongName(hDriver));
+    }
+}
+
+static char *format_list(void)
+{
+    char *buf, *p;
+    int len = 0;
+    int first = 1;
+    int i;
+
+    for (i = 0; i < GDALGetDriverCount(); i++) {
+	GDALDriverH driver = GDALGetDriver(i);
+
+	if (!GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL) &&
+	    !GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL))
+	    continue;
+
+	len += strlen(GDALGetDriverShortName(driver)) + 1;
+    }
+
+    buf = G_malloc(len);
+    p = buf;
+
+    for (i = 0; i < GDALGetDriverCount(); i++) {
+	GDALDriverH driver = GDALGetDriver(i);
+	const char *name;
+
+	if (!GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL) &&
+	    !GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL))
+	    continue;
+
+	if (first)
+	    first = 0;
+	else
+	    *p++ = ',';
+
+	name = GDALGetDriverShortName(driver);
+	strcpy(p, name);
+	p += strlen(name);
+    }
+    *p++ = '\0';
+
+    return buf;
+}
+
+static void check_format(const char *format)
+{
+    GDALDriverH driver = GDALGetDriverByName(format);
+
+    if (!driver)
+	G_fatal_error(_("Format <%s> not supported"), format);
+
+    if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, NULL))
+	return;
+
+    if (GDALGetMetadataItem(driver, GDAL_DCAP_CREATECOPY, NULL)) {
+	G_warning(_("Format <%s> does not support direct write"), format);
+	return;
+    }
+
+    G_fatal_error(_("Format <%s> does not support writing"), format);
+}
+
+static void make_link(const char *dir, const char *ext,
+		      const char *format, char **options)
+{
+    struct Key_Value *key_val = G_create_key_value();
+    char *opt_str = NULL;
+    FILE *fp;
+
+    if (options && *options) {
+	int n_opts = 0, opt_len = 0, i;
+	char *p;
+
+	for (i = 0; options[i]; i++) {
+	    n_opts++;
+	    opt_len += strlen(options[i]) + 1;
+	}
+
+	opt_str = G_malloc(opt_len);
+	p = opt_str;
+
+	for (i = 0; i < n_opts; i++) {
+	    if (i > 0)
+		*p++ = ',';
+	    strcpy(p, options[i]);
+	    p += strlen(options[i]);
+	}
+	*p++ = '\0';
+    }
+
+    if (ext && ext[0] != '.') {
+	char *p;
+	G_asprintf(&p, ".%s", ext);
+	ext = p;
+    }
+
+    if (dir)
+	G_set_key_value("directory", dir, key_val);
+    if (ext)
+	G_set_key_value("extension", ext, key_val);
+    if (format)
+	G_set_key_value("format", format, key_val);
+    if (opt_str)
+	G_set_key_value("options", opt_str, key_val);
+
+    fp = G_fopen_new("", "GDAL");
+    if (!fp)
+	G_fatal_error(_("Unable to create GDAL file"));
+
+    if (G_fwrite_key_value(fp, key_val) < 0)
+	G_fatal_error(_("Error writing GDAL file"));
+
+    fclose(fp);
+}
+
+int main(int argc, char *argv[])
+{
+    struct GModule *module;
+    struct {
+	struct Option *dir, *ext, *format, *opts;
+    } parm;
+    struct Flag *flag_f, *flag_r;
+
+    G_gisinit(argv[0]);
+
+    GDALAllRegister();
+
+    module = G_define_module();
+    module->keywords = _("raster, import");
+    module->description =
+	_("Link GDAL supported raster file to a binary raster map layer.");
+
+    parm.dir = G_define_option();
+    parm.dir->key = "directory";
+    parm.dir->description = _("Name of output directory");
+    parm.dir->required = NO;
+    parm.dir->type = TYPE_STRING;
+
+    parm.ext = G_define_option();
+    parm.ext->key = "extension";
+    parm.ext->description = _("Extension for output files");
+    parm.ext->required = NO;
+    parm.ext->type = TYPE_STRING;
+
+    parm.format = G_define_option();
+    parm.format->key = "format";
+    parm.format->description = _("Format of output files");
+    parm.format->required = NO;
+    parm.format->type = TYPE_STRING;
+    parm.format->options = format_list();
+
+    parm.opts = G_define_option();
+    parm.opts->key = "options";
+    parm.opts->description = _("Creation options");
+    parm.opts->required = NO;
+    parm.opts->multiple = YES;
+    parm.opts->type = TYPE_STRING;
+
+    flag_f = G_define_flag();
+    flag_f->key = 'f';
+    flag_f->description = _("List supported formats and exit");
+    flag_f->guisection = _("Print");
+
+    flag_r = G_define_flag();
+    flag_r->key = 'r';
+    flag_r->description = _("Cease using GDAL and revert to native output");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    if (flag_f->answer) {
+	list_formats();
+	exit(EXIT_SUCCESS);
+    }
+
+    if (flag_r->answer) {
+	G_remove("", "GDAL");
+	exit(EXIT_SUCCESS);
+    }
+
+    if (parm.format->answer)
+	check_format(parm.format->answer);
+
+    make_link(parm.dir->answer, parm.ext->answer, parm.format->answer,
+	      parm.opts->answers);
+
+    exit(EXIT_SUCCESS);
+}
+

raster/r.external.out/r.external.out.html

@@ -0,0 +1,22 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.external.out</em> instructs GRASS to write raster maps as data
+files using GDAL.
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="r.in.gdal.html">r.in.gdal</a>,
+<a href="r.out.gdal.html">r.in.gdal</a>,
+<a href="r.external.html">r.external</a>
+</em>
+
+<h2>REFERENCES</h2>
+
+GDAL Pages: <a href="http://www.gdal.org">http://www.gdal.org/</a><br>
+
+<h2>AUTHOR</h2>
+
+Glynn Clements
+
+<p>
+<i>Last changed: $Date: 2008-08-15 08:16:42 +0200 (Fri, 15 Aug 2008) $</i>