[GRASS-dev] some questions about future development
Glynn Clements
glynn at gclements.plus.com
Tue Aug 19 18:46:13 EDT 2008
Markus Neteler wrote:
> > BTW, a r.external would be a great addition.
>
> Yes, see:
> http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_5_4/gips/gip-0002.txt
>
> Still not implemented... next GSoC project?
It's not that much work.
I have attached a "first draft" patch which adds the essential support
to lib/gis. It has only been minimally tested, but seems to work.
To test:
1. Copy a raster map.
2. Export it to an image at its native resolution.
3. Create a cell_misc/<name>/gdal file consisting of two lines: the
filename followed by the band number.
4. Edit the cellhd/<name> file so that "compressed" is 0, and (if it's
a CELL map) "format" is 3.
5. Run "r.support -s <name>" to update the range.
Tomorrow, I'll look into creating an r.external module, adding null
support, etc.
--
Glynn Clements <glynn at gclements.plus.com>
-------------- next part --------------
Index: lib/gis/opencell.c
===================================================================
--- lib/gis/opencell.c (revision 32899)
+++ lib/gis/opencell.c (working copy)
@@ -159,6 +159,9 @@
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 @@
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 @@
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() */
Index: lib/gis/G.h
===================================================================
--- lib/gis/G.h (revision 32899)
+++ lib/gis/G.h (working copy)
@@ -44,6 +44,9 @@
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 */
Index: lib/gis/get_row.c
===================================================================
--- lib/gis/get_row.c (revision 32899)
+++ lib/gis/get_row.c (working copy)
@@ -198,12 +198,35 @@
/*--------------------------------------------------------------------------*/
+#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);
@@ -352,6 +375,36 @@
/*--------------------------------------------------------------------------*/
+static void cell_values_raw(const unsigned char *data,
+ const COLUMN_MAPPING *cmap, int nbytes,
+ void *cell, int n)
+{
+ unsigned char *c = cell;
+ const unsigned char *d;
+ COLUMN_MAPPING cmapold = 0;
+ int i;
+
+ for (i = 0; i < n; i++) {
+ if (!cmap[i]) {
+ memset(c + i*nbytes, 0, nbytes);
+ continue;
+ }
+
+ if (cmap[i] == cmapold) {
+ memcpy(c + i*nbytes, c + (i-1)*nbytes, nbytes);
+ continue;
+ }
+
+ d = data + (cmap[i] - 1) * nbytes;
+
+ memcpy(c + i*nbytes, d, nbytes);
+
+ cmapold = cmap[i];
+ }
+}
+
+/*--------------------------------------------------------------------------*/
+
/* transfer_to_cell_XY takes bytes from fcb->data, converts these bytes with
the appropriate procedure (e.g. XDR or byte reordering) into type X
values which are put into array G__.work_buf.
@@ -370,6 +423,11 @@
cell_values_int, cell_values_float, cell_values_double};
struct fileinfo *fcb = &G__.fileinfo[fd];
+#ifdef GDAL_LINK
+ if (fcb->gdal)
+ cell_values_raw(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);
Index: lib/gis/closecell.c
===================================================================
--- lib/gis/closecell.c (revision 32899)
+++ lib/gis/closecell.c (working copy)
@@ -133,6 +133,11 @@
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);
Index: lib/gis/Makefile
===================================================================
--- lib/gis/Makefile (revision 32899)
+++ lib/gis/Makefile (working copy)
@@ -1,7 +1,7 @@
MODULE_TOPDIR = ../..
LIB_NAME = $(GIS_LIBNAME)
-EXTRA_LIBS = $(XDRLIB) $(SOCKLIB) $(DATETIMELIB) $(INTLLIB) $(MATHLIB)
+EXTRA_LIBS = $(GDALLIBS) $(XDRLIB) $(SOCKLIB) $(DATETIMELIB) $(INTLLIB) $(MATHLIB)
DATASRC = ellipse.table datum.table datumtransform.table FIPS.code state27 state83 projections gui.tcl
include $(MODULE_TOPDIR)/include/Make/Lib.make
@@ -15,6 +15,7 @@
ifneq ($(USE_LARGEFILES),)
EXTRA_CFLAGS = -D_FILE_OFFSET_BITS=64
endif
+EXTRA_CFLAGS += -DGDAL_LINK
default: lib $(FMODE_OBJ) $(DATAFILES) $(COLORFILES)
Index: lib/gis/gdal.c
===================================================================
--- lib/gis/gdal.c (revision 0)
+++ lib/gis/gdal.c (revision 0)
@@ -0,0 +1,86 @@
+#include <grass/config.h>
+#include <grass/gis.h>
+#include <gdal.h>
+
+#ifdef GDAL_LINK
+struct GDAL_link *G_get_gdal_link(const char *name, const char *mapset)
+{
+ static int initialized;
+ char filename[GPATH_MAX];
+ char band_str[64];
+ int band_num;
+ GDALDatasetH data;
+ GDALRasterBandH band;
+ GDALDataType type;
+ struct GDAL_link *gdal;
+ RASTER_MAP_TYPE map_type;
+ FILE *fp;
+ int err = 0;
+
+ 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;
+
+ if (!G_getl2(filename, sizeof(filename), fp))
+ err = 1;
+
+ if (!G_getl2(band_str, sizeof(band_str), fp))
+ err = 1;
+
+ band_num = atoi(band_str);
+ if (!band_num)
+ err = 1;
+
+ fclose(fp);
+
+ if (err)
+ 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;
+ }
+
+ switch (map_type) {
+ case CELL_TYPE: type = GDT_Int32; break;
+ case FCELL_TYPE: type = GDT_Float32; break;
+ case DCELL_TYPE: type = GDT_Float64; break;
+ default: type = GDT_Float64; break;
+ }
+
+ gdal = G_calloc(1, sizeof(struct GDAL_link));
+
+ gdal->filename = G_store(filename);
+ gdal->band_num = band_num;
+ 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
Index: include/gis.h
===================================================================
--- include/gis.h (revision 32899)
+++ include/gis.h (working copy)
@@ -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,17 @@
int count;
};
+#ifdef GDAL_LINK
+struct GDAL_link
+{
+ char *filename;
+ int band_num;
+ GDALDatasetH data;
+ GDALRasterBandH band;
+ GDALDataType type;
+};
+#endif
+
/*============================== Prototypes ================================*/
/* Since there are so many prototypes for the gis library they are stored */
Index: include/gisdefs.h
===================================================================
--- include/gisdefs.h (revision 32899)
+++ include/gisdefs.h (working copy)
@@ -521,6 +521,11 @@
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);
More information about the grass-dev
mailing list