libraster: +GRASS virtual raster (VRT)

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

lib/raster/R.h

@@ -36,6 +36,21 @@ extern CPLErr Rast_gdal_raster_IO(GDALRasterBandH, GDALRWFlag,
 				  void *, int, int, GDALDataType, int, int);
 #endif
 
+struct tileinfo		/* Information for tiles */
+{
+    char *name;			/* Name of open file            */
+    char *mapset;		/* Mapset of open file          */
+    struct Cell_head cellhd;	/* Cell header                  */
+    struct ilist *clist;	/* columns inside current region */
+};
+
+struct R_vrt
+{
+    int tilecount;
+    struct tileinfo *tileinfo;
+    struct ilist *tlist;
+};
+
 struct fileinfo			/* Information for opened cell files */
 {
     int open_mode;		/* see defines below            */
@@ -67,6 +82,7 @@ struct fileinfo			/* Information for opened cell files */
     struct GDAL_link *gdal;
     int data_fd;		/* Raster data fd               */
     off_t *null_row_ptr;	/* Null file row addresses      */
+    struct R_vrt *vrt;
 };
 
 struct R__			/*  Structure of library globals */

lib/raster/close.c

@@ -152,6 +152,8 @@ static int close_old(int fd)
 
     if (fcb->gdal)
 	Rast_close_gdal_link(fcb->gdal);
+    if (fcb->vrt)
+	Rast_close_vrt(fcb->vrt);
 
     if (fcb->null_bits)
 	G_free(fcb->null_bits);

lib/raster/get_row.c

@@ -561,6 +561,11 @@ static int get_map_row_nomask(int fd, void *rast, int row,
     int r;
     int row_status;
 
+    /* is this the best place to read a vrt row, or
+     * call Rast_get_vrt_row() earlier ? */
+    if (fcb->vrt)
+	return Rast_get_vrt_row(fd, rast, row, data_type);
+
     row_status = compute_window_row(fd, row, &r);
 
     if (!row_status) {

lib/raster/open.c

@@ -165,6 +165,7 @@ int Rast__open_old(const char *name, const char *mapset)
     struct Reclass reclass;
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
     struct GDAL_link *gdal;
+    struct R_vrt *vrt;
 
     Rast__init();
 
@@ -268,6 +269,8 @@ int Rast__open_old(const char *name, const char *mapset)
     }
 
     gdal = Rast_get_gdal_link(r_name, r_mapset);
+    vrt = Rast_get_vrt(r_name, r_mapset);
+    cell_fd = -1;
     if (gdal) {
 #ifdef HAVE_GDAL
 	cell_fd = -1;
@@ -276,6 +279,9 @@ int Rast__open_old(const char *name, const char *mapset)
 		      r_name, r_mapset);
 #endif
     }
+    else if (vrt) {
+	cell_fd = -1;
+    }
     else {
 	/* now actually open file for reading */
 	cell_fd = G_open_old(cell_dir, r_name, r_mapset);
@@ -313,16 +319,20 @@ int Rast__open_old(const char *name, const char *mapset)
 	fcb->reclass = reclass;
 
     fcb->gdal = gdal;
-    if (!gdal)
+    fcb->vrt = vrt;
+    if (!gdal && !vrt) {
 	/* check for compressed data format, making initial reads if necessary */
 	if (Rast__check_format(fd) < 0) {
 	    close(cell_fd);	/* warning issued by check_format() */
 	    G_fatal_error(_("Error reading format for <%s@%s>"),
 			  r_name, r_mapset);
 	}
+    }
 
-    /* create the mapping from cell file to window */
-    Rast__create_window_mapping(fd);
+    if (!vrt) {
+	/* create the mapping from cell file to window */
+	Rast__create_window_mapping(fd);
+    }
 
     /*
      * allocate the data buffer
@@ -349,7 +359,7 @@ int Rast__open_old(const char *name, const char *mapset)
     fcb->nbytes = MAP_NBYTES;
     fcb->null_row_ptr = NULL;
 
-    if (!gdal) {
+    if (!gdal && !vrt) {
 	/* First, check for compressed null file */
 	fcb->null_fd = G_open_old_misc("cell_misc", NULL_FILE, r_name, r_mapset);
 	if (fcb->null_fd < 0) {

lib/raster/vrt.c

@@ -0,0 +1,226 @@
+/*!
+  \file lib/raster/gdal.c
+  
+  \brief Raster Library - Utilization of GDAL library.
+  
+  (C) 2010 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.
+  
+  \author Markus Metz
+*/
+
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+
+#include "R.h"
+
+int cmp_wnd(const void *a, const void *b)
+{
+    struct Cell_head *cellhda = &((struct tileinfo *) a)->cellhd;
+    struct Cell_head *cellhdb = &((struct tileinfo *) b)->cellhd;
+
+    /* sort from descending N to S, then ascending from W to E */
+    if (cellhda->south > cellhdb->south)
+	return -1;
+    if (cellhda->south < cellhdb->south)
+	return 1;
+    if (cellhda->north > cellhdb->north)
+	return -1;
+    if (cellhda->north < cellhdb->north)
+	return 1;
+    if (cellhda->west < cellhdb->west)
+	return -1;
+    if (cellhda->west > cellhdb->west)
+	return 1;
+    if (cellhda->east < cellhdb->east)
+	return -1;
+    if (cellhda->east > cellhdb->east)
+	return 1;
+
+    return 0;
+}
+
+struct R_vrt *Rast_get_vrt(const char *vname, const char *vmapset)
+{
+    FILE *fp;
+    int talloc, tilecount;
+    struct tileinfo *ti;
+    struct R_vrt *vrt;
+    struct Cell_head *rd_window = &R__.rd_window;
+    struct ilist *tlist;
+
+    tilecount = 0;
+    ti = NULL;
+
+    if (!G_find_raster2(vname, vmapset))
+	return NULL;
+
+    fp = G_fopen_old_misc("cell_misc", "vrt", vname, vmapset);
+    if (!fp)
+	return NULL;
+
+    tlist = G_new_ilist();
+    talloc = 0;
+    while (1) {
+	char buf[GNAME_MAX];
+	char *name;
+	const char *mapset;
+	struct tileinfo *p;
+
+	if (!G_getl2(buf, sizeof(buf), fp))
+	    break;
+
+	/* Ignore empty lines */
+	if (!*buf)
+	    continue;
+
+	name = buf;
+	if ((mapset = G_find_raster(name, "")) == NULL)
+	    G_fatal_error(_("Tile raster map <%s> not found"), name);
+
+	if (tilecount >= talloc) {
+	    talloc += 100;
+	    ti = G_realloc(ti, talloc * sizeof(struct tileinfo));
+	}
+	p = &ti[tilecount];
+
+	p->name = G_store(name);
+	p->mapset = G_store(mapset);
+	Rast_get_cellhd(p->name, p->mapset, &(p->cellhd));
+	p->clist = NULL;
+
+	if (rd_window->proj == PROJECTION_LL) {
+	    while (p->cellhd.west >= rd_window->east) {
+		p->cellhd.west -= 360.0;
+		p->cellhd.east -= 360.0;
+	    }
+	    while (p->cellhd.east <= rd_window->west) {
+		p->cellhd.west += 360.0;
+		p->cellhd.east += 360.0;
+	    }
+	}
+
+	if (p->cellhd.north > rd_window->south && 
+	    p->cellhd.south <= rd_window->north &&
+	    p->cellhd.west < rd_window->east && 
+	    p->cellhd.east >= rd_window->west) {
+	    
+	    int col;
+	    double east;
+	    
+	    G_ilist_add(tlist, tilecount);
+	    
+	    p->clist = G_new_ilist();
+	    for (col = 0; col < rd_window->cols; col++) {
+		east = rd_window->west + rd_window->ew_res * (col + 0.5);
+		
+		if (rd_window->proj == PROJECTION_LL) {
+		    while (east > p->cellhd.east)
+			east -= 360;
+		    while (east < p->cellhd.west)
+			east += 360;
+		}
+		if (east >= p->cellhd.west && east < p->cellhd.east)
+		    G_ilist_add(p->clist, col);
+	    }
+	}
+	tilecount++;
+    }
+
+    if (tilecount > 1)
+	qsort(ti, tilecount, sizeof(struct tileinfo), cmp_wnd);
+
+    fclose(fp);
+    
+    vrt = G_calloc(1, sizeof(struct R_vrt));
+    vrt->tilecount = tilecount;
+    vrt->tileinfo = ti;
+    vrt->tlist = tlist;
+
+    return vrt;
+}
+
+void Rast_close_vrt(struct R_vrt *vrt)
+{
+    int i;
+    
+    for (i = 0; i < vrt->tilecount; i++) {
+	struct tileinfo *p;
+
+	p = &(vrt->tileinfo[i]);
+	
+	G_free(p->name);
+	G_free(p->mapset);
+	if (p->clist)
+	    G_free_ilist(p->clist);
+    }
+    G_free(vrt->tileinfo);
+    G_free_ilist(vrt->tlist);
+    G_free(vrt);
+}
+
+int Rast_get_vrt_row(int fd, void *buf, int row, RASTER_MAP_TYPE data_type)
+{
+    struct fileinfo *fcb = &R__.fileinfo[fd];
+    struct R_vrt *vrt = fcb->vrt;
+    struct tileinfo *ti = vrt->tileinfo;
+    struct Cell_head *rd_window = &R__.rd_window;
+    double rown, rows;
+    int i, j;
+    int have_tile;
+    void *tmpbuf;
+    size_t size = Rast_cell_size(data_type);
+
+    rown = rd_window->north - rd_window->ns_res * row;
+    rows = rd_window->north - rd_window->ns_res * (row + 1);
+
+    Rast_set_null_value(buf, rd_window->cols, data_type);
+    tmpbuf = Rast_allocate_input_buf(data_type);
+    have_tile = 0;
+
+    for (i = 0; i < vrt->tlist->n_values; i++) {
+	struct tileinfo *p = &ti[vrt->tlist->value[i]];
+
+	if (p->cellhd.north > rows && p->cellhd.south <= rown) {
+	    int tfd;
+	    void *p1, *p2;
+
+	    Rast_set_null_value(tmpbuf, rd_window->cols, data_type);
+	    tfd = Rast_open_old(p->name, p->mapset);
+	    Rast_get_row(tfd, tmpbuf, row, data_type);
+	    Rast_unopen(tfd);
+	    
+	    p1 = buf;
+	    p2 = tmpbuf;
+	    /* restrict to start and end col ? */
+	    for (j = 0; j < p->clist->n_values; j++) {
+		p1 = (unsigned char *)buf + size * p->clist->value[j];
+		p2 = (unsigned char *)tmpbuf + size * p->clist->value[j];
+		
+		if (!Rast_is_null_value(p2, data_type)) {
+		    switch (data_type) {
+		    case CELL_TYPE:
+			*(CELL *) p1 = *(CELL *) p2;
+			break;
+		    case FCELL_TYPE:
+			*(FCELL *) p1 = *(FCELL *) p2;
+			break;
+		    case DCELL_TYPE:
+			*(DCELL *) p1 = *(DCELL *) p2;
+			break;
+		    default:
+			break;
+		    }
+		}
+	    }
+	    have_tile = 1;
+	}
+    }
+    G_free(tmpbuf);
+
+    return have_tile;
+}