vrt.c

Go to the documentation of this file.

/*!
   \file lib/raster/vrt.c
 
   \brief Raster Library - virtual GRASS raster maps.
 
   SPDX-FileCopyrightText: 2010 GRASS Development Team
   SPDX-License-Identifier: GPL-2.0-or-later
 
   \author Markus Metz
 */
 
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
 
#include "R.h"
 
static 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.0;
                    while (east < p->cellhd.west)
                        east += 360.0;
                }
                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 R_vrt *vrt;
    struct tileinfo *ti;
    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);
 
    // R__.fileinfo can be reallocated by concurrent open/close calls,
    // so we extract the vrt pointer under the same lock that protects the I/O
    // below. This expects callers to open the raster once per thread before
    // entering any parallel code.
#pragma omp critical(raster_vrt_read)
    {
        vrt = R__.fileinfo[fd].vrt;
    }
    ti = vrt->tileinfo;
 
    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;
 
    /* parallelised reading of the real raster maps
     * constituting a GRASS virtual raster
     * causes IO read errors and segmentation faults:
     * enforce reading of the different rasters in only one thread
     * Serialize the entire open/read/close cycle for each tile
     * because they all access the non-thread-safe global R__.fileinfo array. */
#pragma omp critical(raster_vrt_read)
    for (i = 0; i < vrt->tlist->n_values; i++) {
        struct tileinfo *p = &ti[vrt->tlist->value[i]];
 
        /* open the tile only if it overlaps with the current reading
         * row and the current reading EW extents */
        if (p->cellhd.north > rows && p->cellhd.south <= rown &&
            p->cellhd.west < rd_window->east &&
            p->cellhd.east > rd_window->west) {
            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);
 
            /* 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)) {
                    memcpy(p1, p2, size);
                }
            }
            have_tile = 1;
        }
    }
    G_free(tmpbuf);
 
    return have_tile;
}