[GRASS-SVN] r72763 - in grass/trunk/raster: . r.buildvrt
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 4 05:26:10 PDT 2018
Author: mmetz
Date: 2018-06-04 05:26:10 -0700 (Mon, 04 Jun 2018)
New Revision: 72763
Added:
grass/trunk/raster/r.buildvrt/
grass/trunk/raster/r.buildvrt/r.buildvrt.html
Removed:
grass/trunk/raster/r.buildvrt/list.c
grass/trunk/raster/r.buildvrt/proj.c
grass/trunk/raster/r.buildvrt/r.external.html
grass/trunk/raster/r.buildvrt/window.c
Modified:
grass/trunk/raster/r.buildvrt/Makefile
grass/trunk/raster/r.buildvrt/link.c
grass/trunk/raster/r.buildvrt/main.c
grass/trunk/raster/r.buildvrt/proto.h
Log:
r.buildvrt: new module to build a GRASS virtual raster (VRT)
Modified: grass/trunk/raster/r.buildvrt/Makefile
===================================================================
--- grass/trunk/raster/r.external/Makefile 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/Makefile 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,13 +1,10 @@
MODULE_TOPDIR = ../..
-PGM = r.external
+PGM = r.buildvrt
-LIBES = $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(GDALLIBS) $(MATHLIB) $(IMAGERYLIB)
-DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP) $(IMAGERYDEP)
-EXTRA_INC = $(PROJINC) $(GDALCFLAGS)
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
-ifneq ($(USE_GDAL),)
default: cmd
-endif
Modified: grass/trunk/raster/r.buildvrt/link.c
===================================================================
--- grass/trunk/raster/r.external/link.c 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/link.c 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,99 +1,11 @@
+#include <stdio.h>
#include <grass/gis.h>
#include <grass/glocale.h>
-#include <gdal.h>
-
#include "proto.h"
-void query_band(GDALRasterBandH hBand, const char *output,
- struct Cell_head *cellhd, struct band_info *info)
+void make_cell(const char *output, int maptype)
{
- int bGotMin, bGotMax;
-
- info->gdal_type = GDALGetRasterDataType(hBand);
-
- info->null_val = GDALGetRasterNoDataValue(hBand, &info->has_null);
-
- cellhd->compressed = 0;
-
- switch (info->gdal_type) {
- case GDT_Float32:
- info->data_type = FCELL_TYPE;
- cellhd->format = -1;
- break;
-
- case GDT_Float64:
- info->data_type = DCELL_TYPE;
- cellhd->format = -1;
- break;
-
- case GDT_Byte:
- info->data_type = CELL_TYPE;
- cellhd->format = 0;
- break;
-
- case GDT_Int16:
- case GDT_UInt16:
- info->data_type = CELL_TYPE;
- cellhd->format = 1;
- break;
-
- case GDT_Int32:
- case GDT_UInt32:
- info->data_type = CELL_TYPE;
- cellhd->format = 3;
- break;
-
- default:
- G_fatal_error(_("Complex types not supported"));
- break;
- }
-
- info->range[0] = GDALGetRasterMinimum(hBand, &bGotMin);
- info->range[1] = GDALGetRasterMaximum(hBand, &bGotMax);
- if(!(bGotMin && bGotMax))
- GDALComputeRasterMinMax(hBand, 0, info->range);
-
- Rast_init_colors(&info->colors);
-
- if (GDALGetRasterColorTable(hBand) != NULL) {
- GDALColorTableH hCT;
- int count, i;
-
- G_verbose_message(_("Copying color table for %s"), output);
-
- hCT = GDALGetRasterColorTable(hBand);
- count = GDALGetColorEntryCount(hCT);
-
- for (i = 0; i < count; i++) {
- GDALColorEntry sEntry;
-
- GDALGetColorEntryAsRGB(hCT, i, &sEntry);
- if (sEntry.c4 == 0)
- continue;
-
- Rast_set_c_color(i, sEntry.c1, sEntry.c2, sEntry.c3, &info->colors);
- }
- }
- else {
- if (info->gdal_type == GDT_Byte) {
- /* set full 0..255 range to grey scale: */
- G_verbose_message(_("Setting grey color table for <%s> (full 8bit range)"),
- output);
- Rast_make_grey_scale_colors(&info->colors, 0, 255);
- }
- else {
- /* set data range to grey scale: */
- G_verbose_message(_("Setting grey color table for <%s> (data range)"),
- output);
- Rast_make_grey_scale_colors(&info->colors,
- (int) info->range[0], (int) info->range[1]);
- }
- }
-}
-
-void make_cell(const char *output, const struct band_info *info)
-{
FILE *fp;
fp = G_fopen_new("cell", output);
@@ -102,7 +14,7 @@
fclose(fp);
- if (info->data_type == CELL_TYPE)
+ if (maptype == CELL_TYPE)
return;
fp = G_fopen_new("fcell", output);
@@ -112,57 +24,37 @@
fclose(fp);
}
-void make_link(const char *input, const char *output, int band,
- const struct band_info *info, int flip)
+void make_link(const struct input *inputs, int num_inputs,
+ const char *output)
{
- struct Key_Value *key_val = G_create_key_value();
- char null_str[256], type_str[8], band_str[8];
+ int i;
FILE *fp;
- sprintf(band_str, "%d", band);
+ fp = G_fopen_new_misc("cell_misc", "vrt", output);
+ if (!fp)
+ G_fatal_error(_("Unable to create cell_misc/%s/vrt file"), output);
- if (info->has_null) {
- if (info->data_type == CELL_TYPE)
- sprintf(null_str, "%d", (int) info->null_val);
- else
- sprintf(null_str, "%.22g", info->null_val);
+ for (i = 0; i < num_inputs; i++) {
+ const struct input *p = &inputs[i];
+
+ fprintf(fp, "%s@%s\n", p->name, p->mapset);
}
- else
- strcpy(null_str, "none");
- sprintf(type_str, "%d", info->gdal_type);
-
- G_set_key_value("file", input, key_val);
- G_set_key_value("band", band_str, key_val);
- G_set_key_value("null", null_str, key_val);
- G_set_key_value("type", type_str, key_val);
- if (flip & FLIP_H)
- G_set_key_value("hflip", "yes", key_val);
- if (flip & FLIP_V)
- G_set_key_value("vflip", "yes", key_val);
-
- fp = G_fopen_new_misc("cell_misc", "gdal", output);
- if (!fp)
- G_fatal_error(_("Unable to create cell_misc/%s/gdal file"), output);
-
- if (G_fwrite_key_value(fp, key_val) < 0)
- G_fatal_error(_("Error writing cell_misc/%s/gdal file"), output);
-
fclose(fp);
}
-void write_fp_format(const char *output, const struct band_info *info)
+void write_fp_format(const char *output, int maptype)
{
struct Key_Value *key_val;
const char *type;
FILE *fp;
- if (info->data_type == CELL_TYPE)
+ if (maptype == CELL_TYPE)
return;
key_val = G_create_key_value();
- type = (info->data_type == FCELL_TYPE)
+ type = (maptype == FCELL_TYPE)
? "float"
: "double";
G_set_key_value("type", type, key_val);
@@ -191,35 +83,39 @@
Rast_write_quant(output, G_mapset(), &quant);
}
-void create_map(const char *input, int band, const char *output,
- struct Cell_head *cellhd, struct band_info *info,
- const char *title, int flip)
+void create_map(const struct input *inputs, int num_inputs, const char *output,
+ struct Cell_head *cellhd, int maptype, DCELL dmin, DCELL dmax,
+ int have_stats, struct R_stats *ostats, const char *title)
{
struct History history;
struct Categories cats;
+ struct Colors colors;
Rast_put_cellhd(output, cellhd);
- make_cell(output, info);
+ make_cell(output, maptype);
- make_link(input, output, band, info, flip);
+ make_link(inputs, num_inputs, output);
- if (info->data_type == CELL_TYPE) {
+ if (maptype == CELL_TYPE) {
struct Range range;
- range.min = (CELL)info->range[0];
- range.max = (CELL)info->range[1];
+ range.min = (CELL)dmin;
+ range.max = (CELL)dmax;
range.first_time = 0;
Rast_write_range(output, &range);
}
else {
struct FPRange fprange;
- fprange.min = info->range[0];
- fprange.max = info->range[1];
+ fprange.min = dmin;
+ fprange.max = dmax;
fprange.first_time = 0;
Rast_write_fp_range(output, &fprange);
- write_fp_format(output, info);
+ write_fp_format(output, maptype);
write_fp_quant(output);
}
+ G_remove_misc("cell_misc", "stats", output);
+ if (have_stats)
+ Rast_write_rstats(output, ostats);
G_verbose_message(_("Creating support files for %s"), output);
Rast_short_history(output, "raster", &history);
@@ -226,7 +122,8 @@
Rast_command_history(&history);
Rast_write_history(output, &history);
- Rast_write_colors(output, G_mapset(), &info->colors);
+ if (Rast_read_colors(inputs[0].name, inputs[0].mapset, &colors) == 1)
+ Rast_write_colors(output, G_mapset(), &colors);
Rast_init_cats(NULL, &cats);
Rast_write_cats((char *)output, &cats);
Deleted: grass/trunk/raster/r.buildvrt/list.c
===================================================================
--- grass/trunk/raster/r.external/list.c 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/list.c 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,93 +0,0 @@
-#include <stdio.h>
-
-#include <grass/gis.h>
-#include <grass/gprojects.h>
-#include <grass/glocale.h>
-
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-#include <gdal.h>
-
-#include "proto.h"
-
-void list_formats(void)
-{
- /* -------------------------------------------------------------------- */
- /* List supported formats and exit. */
- /* code from GDAL 1.2.5 gcore/gdal_misc.cpp */
- /* Copyright (c) 1999, Frank Warmerdam */
- /* -------------------------------------------------------------------- */
- int iDr;
-
- G_message(_("Supported formats:"));
- for (iDr = 0; iDr < GDALGetDriverCount(); iDr++) {
- GDALDriverH hDriver = GDALGetDriver(iDr);
- const char *pszRWFlag;
-
-#ifdef GDAL_DCAP_RASTER
- /* Starting with GDAL 2.0, vector drivers can also be returned */
- /* Only keep raster drivers */
- if (!GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, NULL))
- continue;
-#endif
-
- if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL))
- pszRWFlag = "rw+";
- else if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, NULL))
- pszRWFlag = "rw";
- else
- pszRWFlag = "ro";
-
- fprintf(stdout, " %s (%s): %s\n",
- GDALGetDriverShortName(hDriver),
- pszRWFlag, GDALGetDriverLongName(hDriver));
- }
-}
-
-void list_bands(struct Cell_head *cellhd, GDALDatasetH hDS)
-{
- struct Cell_head loc_wind;
- struct Key_Value *proj_info = NULL, *proj_units = NULL;
- struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
- int n_bands, i_band, proj_same;
- GDALRasterBandH hBand;
- GDALDataType gdal_type;
-
- if (GPJ_wkt_to_grass(cellhd, &proj_info,
- &proj_units, GDALGetProjectionRef(hDS), 0) < 0) {
- proj_same = 0;
- }
- else {
-
- G_get_default_window(&loc_wind);
- if (loc_wind.proj != PROJECTION_XY) {
- loc_proj_info = G_get_projinfo();
- loc_proj_units = G_get_projunits();
- }
-
-
- if (loc_wind.proj != cellhd->proj ||
- (G_compare_projections
- (loc_proj_info, loc_proj_units, proj_info, proj_units)) < 0) {
- proj_same = 0;
- }
- else {
- proj_same = 1;
- }
-
- }
-
- n_bands = GDALGetRasterCount(hDS);
-
- for (i_band = 1; i_band <= n_bands; i_band++) {
-
- hBand = GDALGetRasterBand(hDS, i_band);
- gdal_type = GDALGetRasterDataType(hBand);
-
- fprintf(stdout, "%d,%s,%d\n", i_band, GDALGetDataTypeName(gdal_type),
- proj_same);
-
- }
-}
-
Modified: grass/trunk/raster/r.buildvrt/main.c
===================================================================
--- grass/trunk/raster/r.external/main.c 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/main.c 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,14 +1,14 @@
/****************************************************************************
*
- * MODULE: r.external
+ * MODULE: r.buildvrt
*
- * AUTHOR(S): Glynn Clements, based on r.in.gdal
- * List GDAL layers by Martin Landa <landa.martin gmail.com> 8/2011
+ * AUTHOR(S): Markus Metz, based on r.external
*
- * PURPOSE: Link raster map into GRASS utilizing the GDAL library.
+ * PURPOSE: Build a VRT (Virtual Raster) that is a mosaic of the
+ * list of input raster maps.
*
- * COPYRIGHT: (C) 2008-2015 by Glynn Clements and the GRASS Development Team
+ * COPYRIGHT: (C) 2018 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
@@ -17,71 +17,83 @@
*****************************************************************************/
#include <stdlib.h>
-#include <unistd.h>
+#include <stdio.h>
#include <math.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
-#include <grass/imagery.h>
#include <grass/glocale.h>
-#include <gdal.h>
-#include <ogr_srs_api.h>
-#include <cpl_conv.h>
-
#include "proto.h"
+int cmp_wnd(const void *a, const void *b)
+{
+ struct Cell_head *cellhda = &((struct input *) a)->cellhd;
+ struct Cell_head *cellhdb = &((struct input *) 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;
+}
+
int main(int argc, char *argv[])
{
- const char *input, *source, *output;
+ const char *output;
char *title;
struct Cell_head cellhd;
- GDALDatasetH hDS;
- GDALRasterBandH hBand;
struct GModule *module;
struct {
- struct Option *input, *source, *output, *band, *title;
+ struct Option *input, *file, *output, *title;
} parm;
- struct {
- struct Flag *o, *j, *f, *e, *h, *v, *t, *a;
- } flag;
- int min_band, max_band, band;
- struct band_info info;
- int flip;
- struct Ref reference;
+ int i, j, num_inputs;
+ struct input *inputs = NULL;
+ char nsresstr[1024], ewresstr[1024];
+ int maptype;
+ struct FPRange fprange;
+ DCELL dmin, dmax;
+ int have_stats;
+ struct R_stats rstats, ostats;
G_gisinit(argv[0]);
module = G_define_module();
G_add_keyword(_("raster"));
- G_add_keyword(_("import"));
- G_add_keyword(_("external"));
+ G_add_keyword(_("mosaic"));
+ G_add_keyword(_("virtual raster"));
module->description =
- _("Links GDAL supported raster data as a pseudo GRASS raster map.");
+ _("Build a VRT (Virtual Raster) from the list of input raster maps.");
- parm.input = G_define_standard_option(G_OPT_F_INPUT);
- parm.input->description = _("Name of raster file to be linked");
+ parm.input = G_define_standard_option(G_OPT_R_INPUTS);
+ parm.input->description = _("Name of input raster files");
parm.input->required = NO;
parm.input->guisection = _("Input");
- parm.source = G_define_option();
- parm.source->key = "source";
- parm.source->description = _("Name of non-file GDAL data source");
- parm.source->required = NO;
- parm.source->type = TYPE_STRING;
- parm.source->key_desc = "name";
- parm.source->guisection = _("Input");
-
+ parm.file = G_define_standard_option(G_OPT_F_INPUT);
+ parm.file->key = "file";
+ parm.file->description = _("Input file with one raster map name per line");
+ parm.file->required = NO;
+ parm.file->guisection = _("Input");
+
parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.output->guisection = _("Output");
- parm.band = G_define_option();
- parm.band->key = "band";
- parm.band->type = TYPE_INTEGER;
- parm.band->required = NO;
- parm.band->description = _("Band to select (default is all bands)");
- parm.band->guisection = _("Input");
-
parm.title = G_define_option();
parm.title->key = "title";
parm.title->key_desc = "phrase";
@@ -88,179 +100,211 @@
parm.title->type = TYPE_STRING;
parm.title->required = NO;
parm.title->description = _("Title for resultant raster map");
- parm.title->guisection = _("Metadata");
+ parm.title->guisection = _("Output");
- flag.f = G_define_flag();
- flag.f->key = 'f';
- flag.f->description = _("List supported formats and exit");
- flag.f->guisection = _("Print");
- flag.f->suppress_required = YES;
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
- flag.o = G_define_flag();
- flag.o->key = 'o';
- flag.o->label =
- _("Override projection check (use current location's projection)");
- flag.o->description =
- _("Assume that the dataset has same projection as the current location");
+ if (parm.input->answer && parm.file->answer)
+ G_fatal_error(_("%s= and %s= are mutually exclusive"),
+ parm.input->key, parm.file->key);
+
+ if (!parm.input->answer && !parm.file->answer)
+ G_fatal_error(_("Please specify %s= or %s="),
+ parm.input->key, parm.file->key);
- flag.j = G_define_flag();
- flag.j->key = 'j';
- flag.j->description =
- _("Perform projection check only and exit");
- flag.j->suppress_required = YES;
+ /* read input maps from file */
+ if (parm.file->answer) {
+ FILE *in;
+ int max_inputs;
- flag.e = G_define_flag();
- flag.e->key = 'e';
- flag.e->label = _("Extend region extents based on new dataset");
- flag.e->description = _("Also updates the default region if in the PERMANENT mapset");
+ if (strcmp(parm.file->answer, "-") == 0)
+ in = stdin;
+ else {
+ in = fopen(parm.file->answer, "r");
+ if (!in)
+ G_fatal_error(_("Unable to open input file <%s>"), parm.file->answer);
+ }
+
+ num_inputs = 0;
+ max_inputs = 0;
- flag.a = G_define_flag();
- flag.a->key = 'a';
- flag.a->label = _("Auto-adjustment for lat/lon");
- flag.a->description = _("Attempt to fix small precision errors in resolution and extents");
+ for (;;) {
+ char buf[GNAME_MAX];
+ char *name;
+ const char *mapset;
+ struct input *p;
- flag.h = G_define_flag();
- flag.h->key = 'h';
- flag.h->description = _("Flip horizontally");
+ if (!G_getl2(buf, sizeof(buf), in))
+ break;
- flag.v = G_define_flag();
- flag.v->key = 'v';
- flag.v->description = _("Flip vertically");
+ /* Ignore empty lines */
+ if (!*buf)
+ continue;
- flag.t = G_define_flag();
- flag.t->key = 't';
- flag.t->label =
- _("List available bands including band type in dataset and exit");
- flag.t->description = _("Format: band number,type,projection check");
- flag.t->guisection = _("Print");
- flag.t->suppress_required = YES;
+ name = buf;
+ if ((mapset = G_find_raster(name, "")) == NULL)
+ G_fatal_error(_("Input raster map <%s> not found"), name);
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ Rast_read_fp_range(name, mapset, &fprange);
+ dmin = fprange.min;
+ if (Rast_is_d_null_value(&dmin)) {
+ G_verbose_message(_("Input map <%s@%s> is all NULL, skipping"),
+ name, mapset);
+ continue;
+ }
- GDALAllRegister();
+ if (num_inputs >= max_inputs) {
+ max_inputs += 100;
+ inputs = G_realloc(inputs, max_inputs * sizeof(struct input));
+ }
+ p = &inputs[num_inputs++];
- if (flag.f->answer) {
- list_formats();
- exit(EXIT_SUCCESS);
+ p->name = G_store(name);
+ p->mapset = G_store(mapset);
+ p->maptype = Rast_map_type(p->name, p->mapset);
+ Rast_get_cellhd(p->name, p->mapset, &(p->cellhd));
+ }
+ fclose(in);
+
+ if (num_inputs < 1)
+ G_fatal_error(_("No raster map name found in input file"));
+
+ if (num_inputs == 1)
+ G_fatal_error(_("Only one raster map name found in input file"));
}
+ else {
+ for (i = 0; parm.input->answers[i]; i++)
+ ;
+ num_inputs = i;
- input = parm.input->answer;
- source = parm.source->answer;
- output = parm.output->answer;
+ if (num_inputs < 1)
+ G_fatal_error(_("Raster map not found"));
- flip = 0;
- if (flag.h->answer)
- flip |= FLIP_H;
- if (flag.v->answer)
- flip |= FLIP_V;
+ if (num_inputs == 1)
+ G_fatal_error(_("Only one raster map name found"));
- if (parm.title->answer) {
- title = G_store(parm.title->answer);
- G_strip(title);
- }
- else
- title = NULL;
+ inputs = G_malloc(num_inputs * sizeof(struct input));
- if (!input && !source)
- G_fatal_error(_("%s= or %s= must be given"),
- parm.input->key, parm.source->key);
+ j = 0;
+ for (i = 0; i < num_inputs; i++) {
+ char *name;
+ const char *mapset;
+ struct input *p = &inputs[i];
- if (input && source)
- G_fatal_error(_("%s= and %s= are mutually exclusive"),
- parm.input->key, parm.source->key);
-
- if (input && !G_is_absolute_path(input)) {
- char path[GPATH_MAX], *cwd;
- cwd = CPLGetCurrentDir();
- if (!cwd)
- G_fatal_error(_("Unable to get current working directory"));
-
- G_snprintf(path, GPATH_MAX, "%s%c%s", cwd, HOST_DIRSEP, input);
- input = G_store(path);
- CPLFree(cwd);
- }
+ name = parm.input->answers[i];
+ if ((mapset = G_find_raster(name, "")) == NULL)
+ G_fatal_error(_("Input raster map <%s> not found"), name);
- if (!input)
- input = source;
+ Rast_read_fp_range(name, mapset, &fprange);
+ dmin = fprange.min;
+ if (Rast_is_d_null_value(&dmin)) {
+ G_verbose_message(_("Input map <%s@%s> is all NULL, skipping"),
+ name, mapset);
+ continue;
+ }
- hDS = GDALOpen(input, GA_ReadOnly);
- if (hDS == NULL)
- return 1;
+ p = &inputs[j++];
- setup_window(&cellhd, hDS, &flip);
-
- if (flag.t->answer) {
- list_bands(&cellhd, hDS);
- /* close the GDALDataset to avoid segfault in libgdal */
- GDALClose(hDS);
- exit(EXIT_SUCCESS);
+ p->name = G_store(name);
+ p->mapset = G_store(mapset);
+ p->maptype = Rast_map_type(p->name, p->mapset);
+ Rast_get_cellhd(p->name, p->mapset, &(p->cellhd));
+ }
+ num_inputs = j;
}
- check_projection(&cellhd, hDS, NULL, 0, flag.o->answer, flag.j->answer);
+ qsort(inputs, num_inputs, sizeof(struct input), cmp_wnd);
- if (flag.a->answer && cellhd.proj == PROJECTION_LL) {
- G_adjust_Cell_head(&cellhd, 1, 1);
- G_adjust_window_ll(&cellhd);
+ /* check resolution and maptype of input maps */
+ cellhd = inputs[0].cellhd;
+ cellhd.compressed = 0;
+ G_format_resolution(cellhd.ns_res, nsresstr, G_projection());
+ G_format_resolution(cellhd.ew_res, ewresstr, G_projection());
+ maptype = inputs[0].maptype;
+
+ Rast_set_d_null_value(&dmin, 1);
+ Rast_set_d_null_value(&dmax, 1);
+ if (Rast_read_fp_range(inputs[0].name, inputs[0].mapset, &fprange) == 1) {
+ dmin = fprange.min;
+ dmax = fprange.max;
}
-
- Rast_set_window(&cellhd);
-
- if (parm.band->answer)
- min_band = max_band = atoi(parm.band->answer);
+ Rast_set_d_null_value(&(ostats.sum), 1);
+ Rast_set_d_null_value(&(ostats.sumsq), 1);
+ ostats.count = 0;
+ have_stats = 1;
+ if (Rast_read_rstats(inputs[0].name, inputs[0].mapset, &rstats) == 1) {
+ ostats.sum = rstats.sum;
+ ostats.sumsq = rstats.sumsq;
+ ostats.count = rstats.count;
+ }
else
- min_band = 1, max_band = GDALGetRasterCount(hDS);
+ have_stats = 0;
- G_verbose_message(_("Proceeding with import..."));
+ for (i = 1; i < num_inputs; i++) {
+ char tnsresstr[1024], tewresstr[1024];
+ int tmaptype;
+ struct input *p = &inputs[i];
- if (max_band > min_band) {
- if (I_find_group(output) == 1)
- G_warning(_("Imagery group <%s> already exists and will be overwritten."), output);
- I_init_group_ref(&reference);
- }
+ G_format_resolution(p->cellhd.ns_res, tnsresstr, G_projection());
+ G_format_resolution(p->cellhd.ew_res, tewresstr, G_projection());
+ tmaptype = p->maptype;
- for (band = min_band; band <= max_band; band++) {
- char *output2, *title2 = NULL;
+ if (tmaptype != maptype)
+ G_warning(_("Input maptypes are different"));
+ if (strcmp(nsresstr, tnsresstr) != 0)
+ G_warning(_("Input ns resolutions are different"));
+ if (strcmp(ewresstr, tewresstr) != 0)
+ G_warning(_("Input ns resolutions are different"));
- G_message(_("Reading band %d of %d..."),
- band, GDALGetRasterCount( hDS ));
+ if (cellhd.north < p->cellhd.north)
+ cellhd.north = p->cellhd.north;
+ if (cellhd.south > p->cellhd.south)
+ cellhd.south = p->cellhd.south;
+ if (cellhd.east < p->cellhd.east)
+ cellhd.east = p->cellhd.east;
+ if (cellhd.west > p->cellhd.west)
+ cellhd.west = p->cellhd.west;
- hBand = GDALGetRasterBand(hDS, band);
- if (!hBand)
- G_fatal_error(_("Selected band (%d) does not exist"), band);
-
- if (max_band > min_band) {
- G_asprintf(&output2, "%s.%d", output, band);
- if (title)
- G_asprintf(&title2, "%s (band %d)", title, band);
- G_debug(1, "Adding raster map <%s> to group <%s>", output2, output);
- I_add_file_to_group_ref(output2, G_mapset(), &reference);
+ if (Rast_read_fp_range(p->name, p->mapset, &fprange) == 1) {
+ if (Rast_is_d_null_value(&dmin)) {
+ dmin = fprange.min;
+ dmax = fprange.max;
+ }
+ else {
+ if (dmin > fprange.min)
+ dmin = fprange.min;
+ if (dmax < fprange.max)
+ dmax = fprange.max;
+ }
}
- else {
- output2 = G_store(output);
- if (title)
- title2 = G_store(title);
+ if (have_stats &&
+ Rast_read_rstats(inputs[0].name, inputs[0].mapset, &rstats) == 1) {
+ ostats.sum += rstats.sum;
+ ostats.sumsq += rstats.sumsq;
+ ostats.count += rstats.count;
}
+ else
+ have_stats = 0;
+ }
- query_band(hBand, output2, &cellhd, &info);
- create_map(input, band, output2, &cellhd, &info, title, flip);
+ G_adjust_Cell_head(&cellhd, 0, 0);
- G_free(output2);
- G_free(title2);
- }
+ if (maptype == CELL_TYPE)
+ cellhd.format = 3;
+ else
+ cellhd.format = -1;
- /* close the GDALDataset to avoid segfault in libgdal */
- GDALClose(hDS);
+ output = parm.output->answer;
- if (flag.e->answer)
- update_default_window(&cellhd);
-
- /* Create the imagery group if multiple bands are imported */
- if (max_band > min_band) {
- I_put_group_ref(output, &reference);
- I_put_group(output);
- G_message(_("Imagery group <%s> created"), output);
+ title = NULL;
+ if (parm.title->answer) {
+ title = G_store(parm.title->answer);
+ G_strip(title);
}
+ create_map(inputs, num_inputs, output, &cellhd, maptype, dmin, dmax,
+ have_stats, &ostats, title);
+
exit(EXIT_SUCCESS);
}
Deleted: grass/trunk/raster/r.buildvrt/proj.c
===================================================================
--- grass/trunk/raster/r.external/proj.c 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/proj.c 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,290 +0,0 @@
-#include <grass/gis.h>
-#include <grass/gprojects.h>
-#include <grass/glocale.h>
-
-#include <gdal.h>
-#include <ogr_srs_api.h>
-
-/* keep in sync with r.in.gdal, v.in.ogr, v.external */
-void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
- char *outloc, int create_only, int override,
- int check_only)
-{
- struct Cell_head loc_wind;
- struct Key_Value *proj_info, *proj_units, *proj_epsg;
- struct Key_Value *loc_proj_info, *loc_proj_units;
- const char *wkt;
- char error_msg[8096];
- int proj_trouble;
-
- /* -------------------------------------------------------------------- */
- /* Fetch the projection in GRASS form. */
- /* -------------------------------------------------------------------- */
- proj_info = NULL;
- proj_units = NULL;
- proj_epsg = NULL;
- loc_proj_info = NULL;
- loc_proj_units = NULL;
-
- wkt = GDALGetProjectionRef(hDS);
- /* proj_trouble:
- * 0: valid srs
- * 1: no srs, default to xy
- * 2: unreadable srs, default to xy
- */
-
- /* Projection only required for checking so convert non-interactively */
- proj_trouble = 0;
- if (wkt && *wkt) {
- OGRSpatialReferenceH hSRS;
-
- hSRS = OSRNewSpatialReference(wkt);
- if (hSRS != NULL)
- GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0);
-
- if (!hSRS || (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS))) {
- G_important_message(_("Input contains an invalid SRS. "
- "WKT definition:\n%s"), wkt);
-
- proj_trouble = 2;
- }
- else{
- const char *authkey, *authname, *authcode;
-
- if (OSRIsProjected(hSRS))
- authkey = "PROJCS";
- else /* is geographic */
- authkey = "GEOGCS";
-
- authname = OSRGetAuthorityName(hSRS, authkey);
- if (authname && *authname && strcmp(authname, "EPSG") == 0) {
- authcode = OSRGetAuthorityCode(hSRS, authkey);
- if (authcode && *authcode) {
- G_debug(1, "found EPSG:%s", authcode);
- proj_epsg = G_create_key_value();
- G_set_key_value("epsg", authcode, proj_epsg);
- }
- }
- }
- if (hSRS)
- OSRDestroySpatialReference(hSRS);
- }
- else {
- G_important_message(_("No projection information available"));
- cellhd->proj = PROJECTION_XY;
- cellhd->zone = 0;
- proj_trouble = 1;
- }
-
- /* -------------------------------------------------------------------- */
- /* Do we need to create a new location? */
- /* -------------------------------------------------------------------- */
- if (outloc != NULL) {
- /* do not create a xy location if an existing SRS was unreadable */
- if (proj_trouble == 2) {
- G_fatal_error(_("Unable to convert input map projection to GRASS "
- "format; cannot create new location."));
- }
- else {
- if (0 != G_make_location_epsg(outloc, cellhd, proj_info,
- proj_units, proj_epsg)) {
- G_fatal_error(_("Unable to create new location <%s>"),
- outloc);
- }
- G_message(_("Location <%s> created"), outloc);
-
- G_unset_window(); /* new location, projection, and window */
- G_get_window(cellhd);
- }
-
- /* If create only, clean up and exit here */
- if (create_only) {
- GDALClose(hDS);
- exit(EXIT_SUCCESS);
- }
- }
- else {
- int err = 0;
- void (*msg_fn)(const char *, ...);
-
- if (check_only && override) {
- /* can't check when over-riding check */
- override = 0;
- }
-
- if (proj_trouble == 2) {
- strcpy(error_msg, _("Unable to convert input map projection information "
- "to GRASS format."));
- if (override) {
- msg_fn = G_warning;
- }
- else {
- msg_fn = G_fatal_error;
- }
- msg_fn(error_msg);
- if (!override) {
- exit(EXIT_FAILURE);
- }
- }
-
- /* -------------------------------------------------------------------- */
- /* Does the projection of the current location match the */
- /* dataset? */
- /* -------------------------------------------------------------------- */
- G_get_default_window(&loc_wind);
- /* fetch LOCATION PROJ info */
- if (loc_wind.proj != PROJECTION_XY) {
- loc_proj_info = G_get_projinfo();
- loc_proj_units = G_get_projunits();
- }
-
- if (override) {
- cellhd->proj = loc_wind.proj;
- cellhd->zone = loc_wind.zone;
- G_message(_("Over-riding projection check"));
- }
- else if (loc_wind.proj != cellhd->proj
- || (err =
- G_compare_projections(loc_proj_info, loc_proj_units,
- proj_info, proj_units)) != 1) {
- int i_value;
-
- strcpy(error_msg,
- _("Projection of dataset does not"
- " appear to match current location.\n\n"));
-
- /* TODO: output this info sorted by key: */
- if (loc_wind.proj != cellhd->proj || err != -2) {
- /* error in proj_info */
- if (loc_proj_info != NULL) {
- strcat(error_msg, _("Location PROJ_INFO is:\n"));
- for (i_value = 0; i_value < loc_proj_info->nitems;
- i_value++)
- sprintf(error_msg + strlen(error_msg), "%s: %s\n",
- loc_proj_info->key[i_value],
- loc_proj_info->value[i_value]);
- strcat(error_msg, "\n");
- }
-
- if (proj_info != NULL) {
- strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
- for (i_value = 0; i_value < proj_info->nitems; i_value++)
- sprintf(error_msg + strlen(error_msg), "%s: %s\n",
- proj_info->key[i_value],
- proj_info->value[i_value]);
- }
- else {
- strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
- if (cellhd->proj == PROJECTION_XY)
- sprintf(error_msg + strlen(error_msg),
- "Dataset proj = %d (unreferenced/unknown)\n",
- cellhd->proj);
- else if (cellhd->proj == PROJECTION_LL)
- sprintf(error_msg + strlen(error_msg),
- "Dataset proj = %d (lat/long)\n",
- cellhd->proj);
- else if (cellhd->proj == PROJECTION_UTM)
- sprintf(error_msg + strlen(error_msg),
- "Dataset proj = %d (UTM), zone = %d\n",
- cellhd->proj, cellhd->zone);
- else
- sprintf(error_msg + strlen(error_msg),
- "Dataset proj = %d (unknown), zone = %d\n",
- cellhd->proj, cellhd->zone);
- }
- if (loc_wind.proj != cellhd->proj) {
- strcat(error_msg, "\nERROR: proj\n");
- }
- else {
- strcat(error_msg, "\nERROR: ");
- switch (err) {
- case -1:
- strcat(error_msg, "proj\n");
- break;
- case -2:
- strcat(error_msg, "units\n");
- break;
- case -3:
- strcat(error_msg, "datum\n");
- break;
- case -4:
- strcat(error_msg, "ellps, a, es\n");
- break;
- case -5:
- strcat(error_msg, "zone\n");
- break;
- case -6:
- strcat(error_msg, "south\n");
- break;
- case -7:
- strcat(error_msg, "x_0\n");
- break;
- case -8:
- strcat(error_msg, "y_0\n");
- break;
- case -9:
- strcat(error_msg, "lon_0\n");
- break;
- case -10:
- strcat(error_msg, "lat_0\n");
- break;
- case -11:
- strcat(error_msg, "lat_1, lat2\n");
- break;
- }
- }
- }
- else {
- /* error in proj_units */
- if (loc_proj_units != NULL) {
- strcat(error_msg, "Location PROJ_UNITS is:\n");
- for (i_value = 0; i_value < loc_proj_units->nitems;
- i_value++)
- sprintf(error_msg + strlen(error_msg), "%s: %s\n",
- loc_proj_units->key[i_value],
- loc_proj_units->value[i_value]);
- strcat(error_msg, "\n");
- }
-
- if (proj_units != NULL) {
- strcat(error_msg, "Dataset PROJ_UNITS is:\n");
- for (i_value = 0; i_value < proj_units->nitems; i_value++)
- sprintf(error_msg + strlen(error_msg), "%s: %s\n",
- proj_units->key[i_value],
- proj_units->value[i_value]);
- }
- }
- if (!check_only) {
- strcat(error_msg,
- _("\nIn case of no significant differences in the projection definitions,"
- " use the -o flag to ignore them and use"
- " current location definition.\n"));
- strcat(error_msg,
- _("Consider generating a new location from the input dataset using "
- "the 'location' parameter.\n"));
- }
-
- if (check_only)
- msg_fn = G_warning;
- else
- msg_fn = G_fatal_error;
- msg_fn("%s", error_msg);
- if (check_only) {
- GDALClose(hDS);
- exit(EXIT_FAILURE);
- }
- }
- else {
- if (check_only)
- msg_fn = G_message;
- else
- msg_fn = G_verbose_message;
- msg_fn(_("Projection of input dataset and current location "
- "appear to match"));
- if (check_only) {
- GDALClose(hDS);
- exit(EXIT_SUCCESS);
- }
- }
- }
-}
Modified: grass/trunk/raster/r.buildvrt/proto.h
===================================================================
--- grass/trunk/raster/r.external/proto.h 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/proto.h 2018-06-04 12:26:10 UTC (rev 72763)
@@ -5,41 +5,19 @@
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
-struct band_info
+struct input
{
- RASTER_MAP_TYPE data_type;
- GDALDataType gdal_type;
- int has_null;
- double null_val;
- double range[2];
- struct Colors colors;
+ const char *name;
+ const char *mapset;
+ int maptype;
+ struct Cell_head cellhd;
};
-enum flip {
- FLIP_H = 1,
- FLIP_V = 2,
-};
-
/* link.c */
-void query_band(GDALRasterBandH, const char *,
- struct Cell_head *, struct band_info *);
-void make_cell(const char *, const struct band_info *);
-void make_link(const char *, const char *, int,
- const struct band_info *, int);
-void write_fp_format(const char *, const struct band_info *);
+void make_cell(const char *, int);
+void make_link(const struct input *, int, const char *);
+void write_fp_format(const char *, int);
void write_fp_quant(const char *);
-void create_map(const char *, int, const char *,
- struct Cell_head *, struct band_info *,
- const char *, int);
-
-/* list.c */
-void list_layers(FILE *, const char *);
-void list_formats(void);
-void list_bands(struct Cell_head *, GDALDatasetH);
-
-/* proj.c */
-void check_projection(struct Cell_head *, GDALDatasetH, char *, int, int, int);
-
-/* window.c */
-void setup_window(struct Cell_head *, GDALDatasetH, int *);
-void update_default_window(struct Cell_head *);
+void create_map(const struct input *, int, const char *,
+ struct Cell_head *, int, DCELL, DCELL,
+ int, struct R_stats *, const char *);
Copied: grass/trunk/raster/r.buildvrt/r.buildvrt.html (from rev 72749, grass/trunk/raster/r.external/r.external.html)
===================================================================
--- grass/trunk/raster/r.buildvrt/r.buildvrt.html (rev 0)
+++ grass/trunk/raster/r.buildvrt/r.buildvrt.html 2018-06-04 12:26:10 UTC (rev 72763)
@@ -0,0 +1,53 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.buildvrt</em> builds a virtual raster (VRT) that is a mosaic of
+the list of input raster maps. The purpose of such a VRT is to provide
+fast access to small subsets of the VRT, also with multiple simultaneous
+read requests.
+
+<h2>NOTES</h2>
+
+<em>r.buildvrt</em> creates a list of raster maps that can be
+located in different mapsets. The ouput is a read-only link to
+the original raster maps which is only valid if the original raster
+maps remain at the originally indicated mapset. A VRT can also be built
+from raster maps registered with <em>r.external</em>.
+
+<p>
+Reading the whole VRT is slower than reading the equivalent single
+raster map. Only reading small parts of the VRT provides a performance
+benefit.
+
+
+<h2>EXAMPLES</h2>
+
+<h3>VRT from a DEM in the North Carolina sample dataset</h3>
+
+<div class="code"><pre>
+# set the region
+g.region raster=elev_state_500m -p
+# higher resolution
+g.region res=50 -p
+# resample the DEM to 50 meter
+r.resamp.interp input=elev_state_500m output=elev_state_50m method=bilinear
+# create tiles
+r.tile input=elev_state_50m output=elev_state_50m_tile_ width=1000 height=1000 overlap=0
+# dump list of tiles to a file
+g.list type=raster pattern=elev_state_50m_tile_* > tilelist
+# build a vrt
+r.buildvrt file=tilelist output=elev_state_50m_vrt
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.tile.html">r.tile</a>,
+<a href="r.patch.html">r.patch</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Markus Metz
+
+<p>
+<i>Last changed: $Date$</i>
Deleted: grass/trunk/raster/r.buildvrt/r.external.html
===================================================================
--- grass/trunk/raster/r.external/r.external.html 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/r.external.html 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,86 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.external</em> allows a user to link a GDAL supported raster file to a binary
-raster map layer, from any GDAL supported raster map format, with an optional
-title. The file is not imported but just registered as GRASS raster map.
-
-<h2>NOTES</h2>
-
-In essence, <em>r.external</em> creates a read-only link to the
-original dataset which is only valid if the original dataset remains
-at the originally indicated directory and filename.
-
-<h2>NULL data handling</h2>
-
-GDAL-linked (<em>r.external</em>) maps do not have or use a NULL
-bitmap, hence <em>r.null</em> cannot manipulate them directly. Here
-NULL cells are those whose value matches the value reported by the
-GDALGetRasterNoDataValue() function.
-
-To apply the GDAL-linked the user need to either create a MASK (e.g.
-with <em>r.mask</em>) and then "apply" it using e.g. <em>r.resample</em>,
-or use <em>r.mapcalc</em> to create a copy with the appropriate categories
-changed to NULL (if() condition).
-
-<h2>EXAMPLES</h2>
-
-<h3>RGB Orthophoto from GeoTIFF</h3>
-
-<div class="code"><pre>
-# import of all channels (each channel will become a GRASS raster map):
-r.external input=/home/user/data/maps/059100.tif output=ortho
-g.region raster=ortho.3 -p
-d.rgb r=ortho.1 g=ortho.2 b=ortho.3
-r.composite r=ortho.1 g=ortho.2 b=ortho.3 output=ortho.rgb
-</pre></div>
-
-<h3>Processing workflow without data import and export</h3>
-
-External raster maps to be processed can be directly linked using <em>r.external</em>;
-likewise, results can be written out to standard raster formats with
-<em>r.external.out</em> (GDAL supported formats):
-
-<div class="code"><pre>
-# register GeoTIFF file to be used in current mapset:
-r.external input=terra_lst1km20030314.LST_Day.tif output=modis_celsius
-
-# define output directory for files resulting from GRASS calculation:
-r.external.out directory=$HOME/gisoutput/ format="GTiff"
-
-# perform GRASS calculation (here: extract pixels > 20 deg C)
-# this stores the output map directly as GeoTIFF:
-r.mapcalc "warm.tif = if(modis_celsius > 20.0, modis_celsius, null() )"
-
-# cease GDAL output connection and turn back to write GRASS raster files:
-r.external.out -r
-
-# now use the resulting file elsewhere
-gdalinfo $HOME/gisoutput/warm.tif
-</pre></div>
-
-<h2>REFERENCES</h2>
-
-GDAL Pages: <a href="http://www.gdal.org">http://www.gdal.org/</a><br>
-
-<h2>SEE ALSO</h2>
-
-<em>
-<a href="r.import.html">r.import</a>,
-<a href="r.in.gdal.html">r.in.gdal</a>,
-<a href="r.external.out.html">r.external.out</a>
-</em>
-
-<p>
-<em>
-<a href="v.import.html">v.import</a>,
-<a href="v.in.ogr.html">v.in.ogr</a>,
-<a href="v.external.html">v.external</a>,
-<a href="v.external.out.html">v.external.out</a>
-</em>
-
-<h2>AUTHOR</h2>
-
-Glynn Clements
-
-<p>
-<i>Last changed: $Date$</i>
Deleted: grass/trunk/raster/r.buildvrt/window.c
===================================================================
--- grass/trunk/raster/r.external/window.c 2018-05-28 21:42:31 UTC (rev 72749)
+++ grass/trunk/raster/r.buildvrt/window.c 2018-06-04 12:26:10 UTC (rev 72763)
@@ -1,97 +0,0 @@
-#include <grass/gis.h>
-#include <grass/glocale.h>
-
-#include <gdal.h>
-
-#include "proto.h"
-
-void setup_window(struct Cell_head *cellhd, GDALDatasetH hDS, int *flip)
-{
- /* -------------------------------------------------------------------- */
- /* Set up the window representing the data we have. */
- /* -------------------------------------------------------------------- */
-
- double adfGeoTransform[6];
-
- cellhd->rows = GDALGetRasterYSize(hDS);
- cellhd->rows3 = GDALGetRasterYSize(hDS);
- cellhd->cols = GDALGetRasterXSize(hDS);
- cellhd->cols3 = GDALGetRasterXSize(hDS);
-
- if (GDALGetGeoTransform(hDS, adfGeoTransform) == CE_None) {
- if (adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0)
- G_fatal_error(_("Input raster map is rotated - cannot import. "
- "You may use 'gdalwarp' to transform the map to North-up."));
- if (adfGeoTransform[1] <= 0.0) {
- G_message(_("Applying horizontal flip"));
- *flip |= FLIP_H;
- }
- if (adfGeoTransform[5] >= 0.0) {
- G_message(_("Applying vertical flip"));
- *flip |= FLIP_V;
- }
-
- cellhd->north = adfGeoTransform[3];
- cellhd->ns_res = fabs(adfGeoTransform[5]);
- cellhd->ns_res3 = fabs(adfGeoTransform[5]);
- cellhd->south = cellhd->north - cellhd->ns_res * cellhd->rows;
- cellhd->west = adfGeoTransform[0];
- cellhd->ew_res = fabs(adfGeoTransform[1]);
- cellhd->ew_res3 = fabs(adfGeoTransform[1]);
- cellhd->east = cellhd->west + cellhd->cols * cellhd->ew_res;
- cellhd->top = 1.;
- cellhd->bottom = 0.;
- cellhd->tb_res = 1.;
- cellhd->depths = 1;
- }
- else {
- cellhd->north = cellhd->rows;
- cellhd->south = 0.0;
- cellhd->ns_res = 1.0;
- cellhd->ns_res3 = 1.0;
- cellhd->west = 0.0;
- cellhd->east = cellhd->cols;
- cellhd->ew_res = 1.0;
- cellhd->ew_res3 = 1.0;
- cellhd->top = 1.;
- cellhd->bottom = 0.;
- cellhd->tb_res = 1.;
- cellhd->depths = 1;
- }
-}
-
-void update_default_window(struct Cell_head *cellhd)
-{
- /* -------------------------------------------------------------------- */
- /* Extend current window based on dataset. */
- /* -------------------------------------------------------------------- */
-
- struct Cell_head cur_wind;
-
- if (strcmp(G_mapset(), "PERMANENT") == 0)
- /* fixme: expand WIND and DEFAULT_WIND independently. (currently
- WIND gets forgotten and DEFAULT_WIND is expanded for both) */
- G_get_default_window(&cur_wind);
- else
- G_get_window(&cur_wind);
-
- cur_wind.north = MAX(cur_wind.north, cellhd->north);
- cur_wind.south = MIN(cur_wind.south, cellhd->south);
- cur_wind.west = MIN(cur_wind.west, cellhd->west);
- cur_wind.east = MAX(cur_wind.east, cellhd->east);
-
- cur_wind.rows = (int)ceil((cur_wind.north - cur_wind.south)
- / cur_wind.ns_res);
- cur_wind.south = cur_wind.north - cur_wind.rows * cur_wind.ns_res;
-
- cur_wind.cols = (int)ceil((cur_wind.east - cur_wind.west)
- / cur_wind.ew_res);
- cur_wind.east = cur_wind.west + cur_wind.cols * cur_wind.ew_res;
-
- if (strcmp(G_mapset(), "PERMANENT") == 0) {
- G_put_element_window(&cur_wind, "", "DEFAULT_WIND");
- G_message(_("Default region for this location updated"));
- }
- G_put_window(&cur_wind);
- G_message(_("Region for the current mapset updated"));
-}
More information about the grass-commit
mailing list