Add r.external and supporting infrastructure

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

include/gis.h

@@ -28,6 +28,10 @@
 #include <grass/config.h>
 #include <grass/datetime.h>
 
+#ifdef GDAL_LINK
+#include <gdal.h>
+#endif
+
 /*=========================== Constants/Defines ============================*/
 
 #if !defined __GNUC__ || __GNUC__ < 2
@@ -623,6 +627,18 @@ struct TimeStamp
     int count;
 };
 
+#ifdef GDAL_LINK
+struct GDAL_link
+{
+    char *filename;
+    int band_num;
+    DCELL null_val;
+    GDALDatasetH data;
+    GDALRasterBandH band;
+    GDALDataType type;
+};
+#endif
+
 /*============================== Prototypes ================================*/
 
 /* Since there are so many prototypes for the gis library they are stored */

include/gisdefs.h

@@ -521,6 +521,11 @@ void G_fpreclass_perform_if(const struct FPReclass *, const CELL *, FCELL *,
 			    int);
 void G_fpreclass_perform_id(const struct FPReclass *, const CELL *, DCELL *,
 			    int);
+/* gdal.c */
+#ifdef GDAL_LINK
+struct GDAL_link *G_get_gdal_link(const char *, const char *);
+void G_close_gdal_link(struct GDAL_link *);
+#endif
 
 /* geodesic.c */
 int G_begin_geodesic_equation(double, double, double, double);

lib/gis/G.h

@@ -44,6 +44,9 @@ struct fileinfo			/* Information for opened cell files */
     unsigned char *null_work_buf;	/* data buffer for reading null rows    */
     int min_null_row;		/* Minimum row null row number in memory */
     struct Quant quant;
+#ifdef GDAL_LINK
+    struct GDAL_link *gdal;
+#endif
 };
 
 struct G__			/*  Structure of library globals */

lib/gis/Makefile

@@ -15,6 +15,11 @@ COLORFILES := $(patsubst colors/%,$(ETC)/colors/%,$(COLORSRC))
 ifneq ($(USE_LARGEFILES),)
 	EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
 endif
+ifneq ($(GDAL_LINK),)
+EXTRA_CFLAGS += -DGDAL_LINK
+EXTRA_INC += $(PROJINC) $(GDALCFLAGS)
+EXTRA_LIBS += $(GDALLIBS)
+endif
 
 default: lib $(FMODE_OBJ) $(DATAFILES) $(COLORFILES)
 

lib/gis/closecell.c

@@ -133,6 +133,11 @@ static int close_old(int fd)
        This is obsolete since now the mask_bus is always allocated
      */
 
+#ifdef GDAL_LINK
+    if (fcb->gdal)
+	G_close_gdal_link(fcb->gdal);
+#endif
+
     for (i = 0; i < NULL_ROWS_INMEM; i++)
 	G_free(fcb->NULL_ROWS[i]);
     G_free(fcb->null_work_buf);

lib/gis/gdal.c

@@ -0,0 +1,119 @@
+#include <grass/config.h>
+#include <grass/gis.h>
+
+#ifdef GDAL_LINK
+
+#include <gdal.h>
+
+struct GDAL_link *G_get_gdal_link(const char *name, const char *mapset)
+{
+    static int initialized;
+    const char *filename;
+    int band_num;
+    GDALDatasetH data;
+    GDALRasterBandH band;
+    GDALDataType type;
+    struct GDAL_link *gdal;
+    RASTER_MAP_TYPE map_type, req_type;
+    FILE *fp;
+    struct Key_Value *key_val;
+    const char *p;
+    DCELL null_val;
+
+    if (!G_find_cell2(name, mapset))
+	return NULL;
+
+    map_type = G_raster_map_type(name, mapset);
+    if (map_type < 0)
+	return NULL;
+
+    fp = G_fopen_old_misc("cell_misc", "gdal", name, mapset);
+    if (!fp)
+	return NULL;
+    key_val = G_fread_key_value(fp);
+    fclose(fp);
+
+    if (!key_val)
+	return NULL;
+
+    filename = G_find_key_value("file", key_val);
+    if (!filename)
+	return NULL;
+
+    p = G_find_key_value("band", key_val);
+    if (!p)
+	return NULL;
+    band_num = atoi(p);
+    if (!band_num)
+	return NULL;
+
+    p = G_find_key_value("null", key_val);
+    if (!p)
+	return NULL;
+    if (strcmp(p, "none") == 0)
+	G_set_d_null_value(&null_val, 1);
+    else
+	null_val = atof(p);
+
+    p = G_find_key_value("type", key_val);
+    if (!p)
+	return NULL;
+    type = atoi(p);
+
+
+    switch (type) {
+    case GDT_Byte:
+    case GDT_Int16:
+    case GDT_UInt16:
+    case GDT_Int32:
+    case GDT_UInt32:
+	req_type = CELL_TYPE;
+	break;
+    case GDT_Float32:
+	req_type = FCELL_TYPE;
+	break;
+    case GDT_Float64:
+	req_type = DCELL_TYPE;
+	break;
+    default:
+	return NULL;
+    }
+
+    if (req_type != map_type)
+	return NULL;
+
+    if (!initialized) {
+	GDALAllRegister();
+	initialized = 1;
+    }
+
+    data = GDALOpen(filename, GA_ReadOnly);
+    if (!data)
+	return NULL;
+
+    band = GDALGetRasterBand(data, band_num);
+    if (!band) {
+	GDALClose(data);
+	return NULL;
+    }
+
+    gdal = G_calloc(1, sizeof(struct GDAL_link));
+
+    gdal->filename = G_store(filename);
+    gdal->band_num = band_num;
+    gdal->null_val = null_val;
+    gdal->data = data;
+    gdal->band = band;
+    gdal->type = type;
+
+    return gdal;
+}
+
+void G_close_gdal_link(struct GDAL_link *gdal)
+{
+    GDALClose(gdal->data);
+    G_free(gdal->filename);
+    G_free(gdal);
+}
+
+#endif

lib/gis/get_row.c

@@ -198,12 +198,35 @@ static int read_data_uncompressed(int fd, int row, unsigned char *data_buf,
 
 /*--------------------------------------------------------------------------*/
 
+#ifdef GDAL_LINK
+static int read_data_gdal(int fd, int row, unsigned char *data_buf, int *nbytes)
+{
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+    CPLErr err;
+
+    *nbytes = fcb->nbytes;
+
+    err = GDALRasterIO(
+	fcb->gdal->band, GF_Read, 0, row, fcb->cellhd.cols, 1, data_buf,
+	fcb->cellhd.cols, 1, fcb->gdal->type, 0, 0);
+
+    return err == CE_None ? 0 : -1;
+}
+#endif
+
+/*--------------------------------------------------------------------------*/
+
 /* Actually read a row of data in */
 
 static int read_data(int fd, int row, unsigned char *data_buf, int *nbytes)
 {
     struct fileinfo *fcb = &G__.fileinfo[fd];
 
+#ifdef GDAL_LINK
+    if (fcb->gdal)
+	return read_data_gdal(fd, row, data_buf, nbytes);
+#endif
+
     if (!fcb->cellhd.compressed)
 	return read_data_uncompressed(fd, row, data_buf, nbytes);
 
@@ -350,6 +373,104 @@ static void cell_values_double(int fd, const unsigned char *data,
 
 /*--------------------------------------------------------------------------*/
 
+#ifdef GDAL_LINK
+
+/*--------------------------------------------------------------------------*/
+
+static void gdal_values_int(int fd, const unsigned char *data,
+			    const COLUMN_MAPPING *cmap, int nbytes,
+			    CELL *cell, int n)
+{
+    struct fileinfo *fcb = &G__.fileinfo[fd];
+    const unsigned char *d;
+    COLUMN_MAPPING cmapold = 0;
+    int i;
+
+    for (i = 0; i < n; i++) {
+	if (!cmap[i]) {
+	    cell[i] = 0;
+	    continue;
+	}
+
+	if (cmap[i] == cmapold) {
+	    cell[i] = cell[i-1];
+	    continue;
+	}
+
+	d = data + (cmap[i] - 1) * nbytes;
+
+	switch (fcb->gdal->type) {
+	case GDT_Byte:	    cell[i] = *(GByte   *)d;	break;
+	case GDT_Int16:	    cell[i] = *(GInt16  *)d;	break;
+	case GDT_UInt16:    cell[i] = *(GUInt16 *)d;	break;
+	case GDT_Int32:	    cell[i] = *(GInt32  *)d;	break;
+	case GDT_UInt32:    cell[i] = *(GUInt32 *)d;	break;
+	default:
+	    /* shouldn't happen */
+	    G_set_c_null_value(&cell[i], 1);
+	    break;
+	}
+
+	cmapold = cmap[i];
+    }
+}
+
+/*--------------------------------------------------------------------------*/
+
+static void gdal_values_float(int fd, const float *data,
+			      const COLUMN_MAPPING *cmap, int nbytes,
+			      FCELL *cell, int n)
+{
+    COLUMN_MAPPING cmapold = 0;
+    int i;
+
+    for (i = 0; i < n; i++) {
+	if (!cmap[i]) {
+	    cell[i] = 0;
+	    continue;
+	}
+
+	if (cmap[i] == cmapold) {
+	    cell[i] = cell[i-1];
+	    continue;
+	}
+
+	cell[i] = data[cmap[i] - 1];
+
+	cmapold = cmap[i];
+    }
+}
+
+/*--------------------------------------------------------------------------*/
+
+static void gdal_values_double(int fd, const double *data,
+			       const COLUMN_MAPPING *cmap, int nbytes,
+			       DCELL *cell, int n)
+{
+    COLUMN_MAPPING cmapold = 0;
+    int i;
+
+    for (i = 0; i < n; i++) {
+	if (!cmap[i]) {
+	    cell[i] = 0;
+	    continue;
+	}
+
+	if (cmap[i] == cmapold) {
+	    cell[i] = cell[i-1];
+	    continue;
+	}
+
+	cell[i] = data[cmap[i] - 1];
+
+	cmapold = cmap[i];
+    }
+}
+
+/*--------------------------------------------------------------------------*/
+
+#endif
+
 /*--------------------------------------------------------------------------*/
 
 /* transfer_to_cell_XY takes bytes from fcb->data, converts these bytes with
@@ -368,8 +489,19 @@ static void transfer_to_cell_XX(int fd, void *cell)
 {
     static void (*cell_values_type[3]) () = {
     cell_values_int, cell_values_float, cell_values_double};
+#ifdef GDAL_LINK
+    static void (*gdal_values_type[3]) () = {
+    gdal_values_int, gdal_values_float, gdal_values_double};
+#endif
     struct fileinfo *fcb = &G__.fileinfo[fd];
 
+#ifdef GDAL_LINK
+    if (fcb->gdal)
+    (gdal_values_type[fcb->map_type]) (fd, fcb->data, fcb->col_map,
+				       fcb->cur_nbytes, cell,
+				       G__.window.cols);
+    else
+#endif
     (cell_values_type[fcb->map_type]) (fd, fcb->data, fcb->col_map,
 				       fcb->cur_nbytes, cell,
 				       G__.window.cols);

lib/gis/opencell.c

@@ -159,6 +159,9 @@ int G__open_cell_old(const char *name, const char *mapset)
     RASTER_MAP_TYPE MAP_TYPE;
     struct Reclass reclass;
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
+#ifdef GDAL_LINK
+    struct GDAL_link *gdal;
+#endif
 
     /* make sure window is set    */
     G__init_window();
@@ -252,6 +255,13 @@ int G__open_cell_old(const char *name, const char *mapset)
 	MAP_NBYTES = CELL_nbytes;
     }
 
+#ifdef GDAL_LINK
+    gdal = G_get_gdal_link(name, mapset);
+    if (gdal)
+	/* dummy descriptor to reserve the fileinfo slot */
+	fd = open("/dev/null", O_RDONLY);
+    else
+#endif
     /* now actually open file for reading */
     fd = G_open_old(cell_dir, r_name, r_mapset);
     if (fd < 0)
@@ -287,6 +297,11 @@ int G__open_cell_old(const char *name, const char *mapset)
     if ((fcb->reclass_flag = reclass_flag))
 	G_copy(&fcb->reclass, &reclass, sizeof(reclass));
 
+#ifdef GDAL_LINK
+    if (gdal)
+	fcb->gdal = gdal;
+    else
+#endif
     /* check for compressed data format, making initial reads if necessary */
     if (G__check_format(fd) < 0) {
 	close(fd);		/* warning issued by check_format() */

raster/r.external/Makefile

@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM  = r.external
+
+LIBES= $(GPROJLIB) $(GISLIB) $(IMAGERYLIB) $(GDALLIBS) $(GMATHLIB)
+DEPENDENCIES = $(GPROJDEP) $(GISDEP) $(IMAGERYDEP) $(GMATHDEP)
+EXTRA_INC = $(PROJINC) $(GDALCFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.external/main.c

@@ -0,0 +1,566 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.in.gdal
+ *               
+ * AUTHOR(S):    Frank Warmerdam (copyright of this file)
+ *               Added optional GCP transformation: Markus Neteler 10/2001
+ *
+ * PURPOSE:      Imports many GIS/image formats into GRASS utilizing the GDAL
+ *               library.
+ *
+ * COPYRIGHT:    (C) 2001 by Frank Warmerdam
+ *
+ *               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"
+#include "ogr_srs_api.h"
+
+#ifndef MIN
+#  define MIN(a,b)      ((a<b) ? a : b)
+#endif
+#ifndef MAX
+#  define MAX(a,b)      ((a>b) ? a : b)
+#endif
+
+struct band_info
+{
+    RASTER_MAP_TYPE data_type;
+    GDALDataType gdal_type;
+    int has_null;
+    double null_val;
+    double range[2];
+    struct Colors colors;
+};
+
+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
+	    pszRWFlag = "ro";
+
+	fprintf(stdout, "  %s (%s): %s\n",
+		GDALGetDriverShortName(hDriver),
+		pszRWFlag, GDALGetDriverLongName(hDriver));
+    }
+}
+
+static void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS, int override)
+{
+    struct Cell_head loc_wind;
+    struct Key_Value *proj_info = NULL, *proj_units = NULL;
+    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
+    int projcomp_error = 0;
+    char error_msg[8096];
+
+    /* Projection only required for checking so convert non-interactively */
+    if (GPJ_wkt_to_grass(cellhd, &proj_info,
+			 &proj_units, GDALGetProjectionRef(hDS), 0) < 0)
+	G_warning(_("Unable to convert input raster map projection information to "
+		    "GRASS format for checking"));
+    else {
+	/* -------------------------------------------------------------------- */
+	/*      Does the projection of the current location match the           */
+	/*      dataset?                                                        */
+	/* -------------------------------------------------------------------- */
+	G_get_window(&loc_wind);
+	if (loc_wind.proj != PROJECTION_XY) {
+	    loc_proj_info = G_get_projinfo();
+	    loc_proj_units = G_get_projunits();
+	}
+
+	if (override) {
+	    cellhd->proj = loc_wind.proj;
+	    cellhd->zone = loc_wind.zone;
+	    G_warning(_("Over-riding projection check"));
+	}
+	else if (loc_wind.proj != cellhd->proj ||
+		 (projcomp_error = G_compare_projections(
+		      loc_proj_info, loc_proj_units, proj_info, proj_units)) < 0) {
+	    int i_value;
+
+	    strcpy(error_msg,
+		   _("Projection of dataset does not"
+		     " appear to match current location.\n\n"));
+
+	    /* TODO: output this info sorted by key: */
+	    if (loc_proj_info != NULL) {
+		strcat(error_msg, _("Location PROJ_INFO is:\n"));
+		for (i_value = 0;
+		     loc_proj_info != NULL &&
+			 i_value < loc_proj_info->nitems; i_value++)
+		    sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+			    loc_proj_info->key[i_value],
+			    loc_proj_info->value[i_value]);
+		strcat(error_msg, "\n");
+	    }
+
+	    if (proj_info != NULL) {
+		strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
+		for (i_value = 0;
+		     proj_info != NULL && i_value < proj_info->nitems;
+		     i_value++)
+		    sprintf(error_msg + strlen(error_msg), "%s: %s\n",
+			    proj_info->key[i_value],
+			    proj_info->value[i_value]);
+		strcat(error_msg, "\nERROR: ");
+		switch (projcomp_error) {
+		case -1:
+		    strcat(error_msg, "proj\n");
+		    break;
+		case -2:
+		    strcat(error_msg, "units\n");
+		    break;
+		case -3:
+		    strcat(error_msg, "datum\n");
+		    break;
+		case -4:
+		    strcat(error_msg, "ellps\n");
+		    break;
+		case -5:
+		    strcat(error_msg, "zone\n");
+		    break;
+		}
+	    }
+	    else {
+		strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
+		if (cellhd->proj == PROJECTION_XY)
+		    sprintf(error_msg + strlen(error_msg),
+			    "cellhd.proj = %d (unreferenced/unknown)\n",
+			    cellhd->proj);
+		else if (cellhd->proj == PROJECTION_LL)
+		    sprintf(error_msg + strlen(error_msg),
+			    "cellhd.proj = %d (lat/long)\n", cellhd->proj);
+		else if (cellhd->proj == PROJECTION_UTM)
+		    sprintf(error_msg + strlen(error_msg),
+			    "cellhd.proj = %d (UTM), zone = %d\n",
+			    cellhd->proj, cellhd->zone);
+		else if (cellhd->proj == PROJECTION_SP)
+		    sprintf(error_msg + strlen(error_msg),
+			    "cellhd.proj = %d (State Plane), zone = %d\n",
+			    cellhd->proj, cellhd->zone);
+		else
+		    sprintf(error_msg + strlen(error_msg),
+			    "cellhd.proj = %d (unknown), zone = %d\n",
+			    cellhd->proj, cellhd->zone);
+	    }
+	    strcat(error_msg,
+		   _("\nYou can use the -o flag to r.in.gdal to override this check and "
+		     "use the location definition for the dataset.\n"));
+	    strcat(error_msg,
+		   _("Consider generating a new location from the input dataset using "
+		     "the 'location' parameter.\n"));
+	    G_fatal_error(error_msg);
+	}
+	else {
+	    G_message(_("Projection of input dataset and current location "
+			"appear to match"));
+	}
+    }
+}
+
+static void setup_window(struct Cell_head *cellhd, GDALDatasetH hDS)
+{
+    /* -------------------------------------------------------------------- */
+    /*      Set up the window representing the data we have.                */
+    /* -------------------------------------------------------------------- */
+
+    double adfGeoTransform[6];
+
+    cellhd->rows = GDALGetRasterYSize(hDS);
+    cellhd->rows3 = GDALGetRasterYSize(hDS);
+    cellhd->cols = GDALGetRasterXSize(hDS);
+    cellhd->cols3 = GDALGetRasterXSize(hDS);
+
+    if (GDALGetGeoTransform(hDS, adfGeoTransform) == CE_None &&
+	adfGeoTransform[5] < 0.0) {
+	if (adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0)
+	    G_fatal_error(_("Input raster map is rotated - cannot import. "
+			    "You may use 'gdalwarp' to transform the map to North-up."));
+
+	cellhd->north = adfGeoTransform[3];
+	cellhd->ns_res = fabs(adfGeoTransform[5]);
+	cellhd->ns_res3 = fabs(adfGeoTransform[5]);
+	cellhd->south = cellhd->north - cellhd->ns_res * cellhd->rows;
+	cellhd->west = adfGeoTransform[0];
+	cellhd->ew_res = adfGeoTransform[1];
+	cellhd->ew_res3 = adfGeoTransform[1];
+	cellhd->east = cellhd->west + cellhd->cols * cellhd->ew_res;
+	cellhd->top = 1.;
+	cellhd->bottom = 0.;
+	cellhd->tb_res = 1.;
+	cellhd->depths = 1;
+    }
+    else {
+	cellhd->north = cellhd->rows;
+	cellhd->south = 0.0;
+	cellhd->ns_res = 1.0;
+	cellhd->ns_res3 = 1.0;
+	cellhd->west = 0.0;
+	cellhd->east = cellhd->cols;
+	cellhd->ew_res = 1.0;
+	cellhd->ew_res3 = 1.0;
+	cellhd->top = 1.;
+	cellhd->bottom = 0.;
+	cellhd->tb_res = 1.;
+	cellhd->depths = 1;
+    }
+}
+
+static void update_default_window(struct Cell_head *cellhd)
+{
+    /* -------------------------------------------------------------------- */
+    /*      Extend current window based on dataset.                         */
+    /* -------------------------------------------------------------------- */
+
+    struct Cell_head def_wind;
+
+    G_get_default_window(&def_wind);
+
+    def_wind.north = MAX(def_wind.north, cellhd->north);
+    def_wind.south = MIN(def_wind.south, cellhd->south);
+    def_wind.west = MIN(def_wind.west, cellhd->west);
+    def_wind.east = MAX(def_wind.east, cellhd->east);
+
+    def_wind.rows = (int)ceil((def_wind.north - def_wind.south)
+			      / def_wind.ns_res);
+    def_wind.south = def_wind.north - def_wind.rows * def_wind.ns_res;
+
+    def_wind.cols = (int)ceil((def_wind.east - def_wind.west)
+			      / def_wind.ew_res);
+    def_wind.east = def_wind.west + def_wind.cols * def_wind.ew_res;
+
+    G__put_window(&def_wind, "../PERMANENT", "DEFAULT_WIND");
+}
+
+static void query_band(GDALRasterBandH hBand, const char *output,
+		       struct Cell_head *cellhd, struct band_info *info)
+{
+    int bGotMin, bGotMax;
+
+    info->gdal_type = GDALGetRasterDataType(hBand);
+
+    info->null_val = GDALGetRasterNoDataValue(hBand, &info->has_null);
+
+    cellhd->compressed = 0;
+
+    switch (info->gdal_type) {
+    case GDT_Float32:
+	info->data_type = FCELL_TYPE;
+	cellhd->format = -1;
+	break;
+
+    case GDT_Float64:
+	info->data_type = DCELL_TYPE;
+	cellhd->format = -1;
+	break;
+
+    case GDT_Byte:
+	info->data_type = CELL_TYPE;
+	cellhd->format = 0;
+	break;
+
+    case GDT_Int16:
+    case GDT_UInt16:
+	info->data_type = CELL_TYPE;
+	cellhd->format = 1;
+	break;
+
+    case GDT_Int32:
+    case GDT_UInt32:
+	info->data_type = CELL_TYPE;
+	cellhd->format = 3;
+	break;
+
+    default:
+	G_fatal_error(_("Complex types not supported"));
+	break;
+    }
+
+    info->range[0] = GDALGetRasterMinimum(hBand, &bGotMin);
+    info->range[1] = GDALGetRasterMaximum(hBand, &bGotMax);
+    if(!(bGotMin && bGotMax))
+	GDALComputeRasterMinMax(hBand, TRUE, info->range);
+
+    G_init_colors(&info->colors);
+
+    if (GDALGetRasterColorTable(hBand) != NULL) {
+	GDALColorTableH hCT;
+	int count, i;
+
+	G_verbose_message(_("Copying color table for %s"), output);
+
+	hCT = GDALGetRasterColorTable(hBand);
+	count = GDALGetColorEntryCount(hCT);
+
+	for (i = 0; i < count; i++) {
+	    GDALColorEntry sEntry;
+
+	    GDALGetColorEntryAsRGB(hCT, i, &sEntry);
+	    if (sEntry.c4 == 0)
+		continue;
+
+	    G_set_color(i, sEntry.c1, sEntry.c2, sEntry.c3, &info->colors);
+	}
+    }
+    else {
+	if (info->gdal_type == GDT_Byte) {
+	    /* set full 0..255 range to grey scale: */
+	    G_verbose_message(_("Setting grey color table for <%s> (full 8bit range)"),
+			      output);
+	    G_make_grey_scale_colors(&info->colors, 0, 255);
+	}
+	else  {
+	    /* set data range to grey scale: */
+	    G_verbose_message(_("Setting grey color table for <%s> (data range)"),
+			      output);
+	    G_make_grey_scale_colors(&info->colors,
+				     (int) info->range[0], (int) info->range[1]);
+	}
+    }
+}
+
+static void make_cell(const char *output, const struct band_info *info)
+{
+    FILE *fp;
+
+    fp = G_fopen_new("cell", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell/%s file"), output);
+
+    fclose(fp);
+
+    if (info->data_type == CELL_TYPE)
+	return;
+
+    fp = G_fopen_new("fcell", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create fcell/%s file"), output);
+
+    fclose(fp);
+}
+
+static void make_link(const char *input, const char *output, int band,
+		      const struct band_info *info)
+{
+    struct Key_Value *key_val = G_create_key_value();
+    char null_str[256], type_str[8], band_str[8];
+    FILE *fp;
+
+    sprintf(band_str, "%d", band);
+
+    if (info->has_null) {
+	if (info->data_type == CELL_TYPE)
+	    sprintf(null_str, "%d", (int) info->null_val);
+	else
+	    sprintf(null_str, "%.22g", info->null_val);
+    }
+    else
+	strcpy(null_str, "none");
+
+    sprintf(type_str, "%d", info->gdal_type);
+
+    G_set_key_value("file", input, key_val);
+    G_set_key_value("band", band_str, key_val);
+    G_set_key_value("null", null_str, key_val);
+    G_set_key_value("type", type_str, key_val);
+
+    fp = G_fopen_new_misc("cell_misc", "gdal", output);
+    if (!fp)
+	G_fatal_error(_("Unable to create cell_misc/%s/gdal file"), output);
+
+    if (G_fwrite_key_value(fp, key_val) < 0)
+	G_fatal_error(_("Error writing cell_misc/%s/gdal file"), output);
+
+    fclose(fp);
+}
+
+static void create_map(const char *input, int band, const char *output,
+		       struct Cell_head *cellhd, struct band_info *info,
+		       const char *title)
+{
+    struct History history;
+
+    G_put_cellhd(output, cellhd);
+
+    make_cell(output, info);
+
+    make_link(input, output, band, info);
+
+    if (info->data_type == CELL_TYPE) {
+	struct Range range;
+	range.min = (CELL)info->range[0];
+	range.max = (CELL)info->range[1];
+	G_write_range(output, &range);
+    }
+    else {
+	struct FPRange fprange;
+	fprange.min = info->range[0];
+	fprange.max = info->range[1];
+	G_write_fp_range(output, &fprange);
+    }
+
+    G_verbose_message(_("Creating support files for %s"), output);
+    G_short_history(output, "raster", &history);
+    G_command_history(&history);
+    G_write_history(output, &history);
+
+    G_write_colors(output, G_mapset(), &info->colors);
+
+    if (title)
+	G_put_cell_title(output, title);
+
+    G_message(_("<%s> created"), output);
+
+    return;
+}
+
+int main(int argc, char *argv[])
+{
+    char *input;
+    char *output;
+    char *title;
+    struct Cell_head cellhd;
+    GDALDatasetH hDS;
+    GDALRasterBandH hBand;
+    struct GModule *module;
+    struct {
+	struct Option *input, *output, *band, *title;
+    } parm;
+    struct Flag *flag_o, *flag_f, *flag_e;
+    int band;
+    struct band_info info;
+
+    G_gisinit(argv[0]);
+
+    module = G_define_module();
+    module->keywords = _("raster, import");
+    module->description =
+	_("Link GDAL supported raster file to a binary raster map layer.");
+
+    parm.input = G_define_standard_option(G_OPT_F_INPUT);
+    parm.input->description = _("Raster file to be linked");
+    parm.input->required = NO;	/* not required because of -f flag */
+    parm.input->guisection = _("Required");
+
+    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+    parm.output->required = NO;	/* not required because of -f flag */
+    parm.output->guisection = _("Required");
+
+    parm.band = G_define_option();
+    parm.band->key = "band";
+    parm.band->type = TYPE_INTEGER;
+    parm.band->required = NO;
+    parm.band->description = _("Band to select");
+    parm.band->answer = "1";
+
+    parm.title = G_define_option();
+    parm.title->key = "title";
+    parm.title->key_desc = "\"phrase\"";
+    parm.title->type = TYPE_STRING;
+    parm.title->required = NO;
+    parm.title->description = _("Title for resultant raster map");
+    parm.title->guisection = _("Metadata");
+
+    flag_o = G_define_flag();
+    flag_o->key = 'o';
+    flag_o->description =
+	_("Override projection (use location's projection)");
+
+    flag_e = G_define_flag();
+    flag_e->key = 'e';
+    flag_e->description = _("Extend location extents based on new dataset");
+
+    flag_f = G_define_flag();
+    flag_f->key = 'f';
+    flag_f->description = _("List supported formats and exit");
+    flag_f->guisection = _("Print");
+
+    if (G_parser(argc, argv))
+	exit(EXIT_FAILURE);
+
+    GDALAllRegister();
+
+    if (flag_f->answer) {
+	list_formats();
+	exit(EXIT_SUCCESS);
+    }
+
+    input = parm.input->answer;
+
+    output = parm.output->answer;
+
+    if (parm.title->answer) {
+	title = G_store(parm.title->answer);
+	G_strip(title);
+    }
+    else
+	title = NULL;
+
+    if (parm.band->answer)
+	band = atoi(parm.band->answer);
+    else
+	band = 1;
+
+    if (!input)
+	G_fatal_error(_("Name for input raster map not specified"));
+
+    if (!output)
+	G_fatal_error(_("Name for output raster map not specified"));
+
+    hDS = GDALOpen(input, GA_ReadOnly);
+    if (hDS == NULL)
+	return 1;
+
+    setup_window(&cellhd, hDS);
+
+    check_projection(&cellhd, hDS, flag_o->answer);
+
+    G_verbose_message(_("Proceeding with import..."));
+
+    hBand = GDALGetRasterBand(hDS, band);
+    if (!hBand)
+	G_fatal_error(_("Selected band (%d) does not exist"), band);
+
+    if (G_set_window(&cellhd) < 0)
+	G_fatal_error(_("Unable to set window"));
+
+    query_band(hBand, output, &cellhd, &info);
+    create_map(input, band, output, &cellhd, &info, title);
+
+    if (flag_e->answer)
+	update_default_window(&cellhd);
+
+    G_done_msg(" ");
+
+    exit(EXIT_SUCCESS);
+}
+

raster/r.external/r.external.html

@@ -0,0 +1 @@
+<!-- this needs to be written -->