Geografic-Information-System/lib/raster/vrt.c

/*!
  \file lib/raster/vrt.c
  
  \brief Raster Library - virtual GRASS raster maps.
  
  (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 (strcmp(name, vname) == 0)
	    G_fatal_error(_("A virtual raster can not contain itself"));

	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);
}

/* must only be called by get_map_row_nomask() 
 * move to get_row.c as read_data_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;

	    /* recurse into get_map_row(), collect data for all tiles 
	     * a mask is applied to the collected data 
	     * after this function returns */
	    Rast_set_null_value(tmpbuf, rd_window->cols, data_type);
	    /* avoid Rast__check_for_auto_masking() */
	    tfd = Rast__open_old(p->name, p->mapset);
	    Rast_get_row_nomask(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;
}