[GRASS-SVN] r72762 - grass/trunk/lib/raster
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 4 05:25:30 PDT 2018
Author: mmetz
Date: 2018-06-04 05:25:30 -0700 (Mon, 04 Jun 2018)
New Revision: 72762
Added:
grass/trunk/lib/raster/vrt.c
Modified:
grass/trunk/lib/raster/R.h
grass/trunk/lib/raster/close.c
grass/trunk/lib/raster/get_row.c
grass/trunk/lib/raster/open.c
Log:
libraster: +GRASS virtual raster (VRT)
Modified: grass/trunk/lib/raster/R.h
===================================================================
--- grass/trunk/lib/raster/R.h 2018-06-04 12:15:27 UTC (rev 72761)
+++ grass/trunk/lib/raster/R.h 2018-06-04 12:25:30 UTC (rev 72762)
@@ -36,6 +36,21 @@
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 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 */
Modified: grass/trunk/lib/raster/close.c
===================================================================
--- grass/trunk/lib/raster/close.c 2018-06-04 12:15:27 UTC (rev 72761)
+++ grass/trunk/lib/raster/close.c 2018-06-04 12:25:30 UTC (rev 72762)
@@ -152,6 +152,8 @@
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);
Modified: grass/trunk/lib/raster/get_row.c
===================================================================
--- grass/trunk/lib/raster/get_row.c 2018-06-04 12:15:27 UTC (rev 72761)
+++ grass/trunk/lib/raster/get_row.c 2018-06-04 12:25:30 UTC (rev 72762)
@@ -561,6 +561,11 @@
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) {
Modified: grass/trunk/lib/raster/open.c
===================================================================
--- grass/trunk/lib/raster/open.c 2018-06-04 12:15:27 UTC (rev 72761)
+++ grass/trunk/lib/raster/open.c 2018-06-04 12:25:30 UTC (rev 72762)
@@ -165,6 +165,7 @@
struct Reclass reclass;
char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
struct GDAL_link *gdal;
+ struct R_vrt *vrt;
Rast__init();
@@ -268,6 +269,8 @@
}
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 @@
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,7 +319,8 @@
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() */
@@ -320,9 +327,12 @@
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 @@
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) {
Added: grass/trunk/lib/raster/vrt.c
===================================================================
--- grass/trunk/lib/raster/vrt.c (rev 0)
+++ grass/trunk/lib/raster/vrt.c 2018-06-04 12:25:30 UTC (rev 72762)
@@ -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;
+}
Property changes on: grass/trunk/lib/raster/vrt.c
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+text/x-csrc
\ No newline at end of property
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
More information about the grass-commit
mailing list