[GRASS-SVN] r68823 - in grass-addons/grass7/raster3d: . r3.what
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 30 19:27:24 PDT 2016
Author: annakrat
Date: 2016-06-30 19:27:24 -0700 (Thu, 30 Jun 2016)
New Revision: 68823
Added:
grass-addons/grass7/raster3d/r3.what/
grass-addons/grass7/raster3d/r3.what/Makefile
grass-addons/grass7/raster3d/r3.what/main.c
grass-addons/grass7/raster3d/r3.what/r3.what.html
Log:
r3.what: new module for querying 3D rasters similar to r.what
Added: grass-addons/grass7/raster3d/r3.what/Makefile
===================================================================
--- grass-addons/grass7/raster3d/r3.what/Makefile (rev 0)
+++ grass-addons/grass7/raster3d/r3.what/Makefile 2016-07-01 02:27:24 UTC (rev 68823)
@@ -0,0 +1,12 @@
+MODULE_TOPDIR = ../..
+
+PGM=r3.what
+
+LIBES = $(RASTER3DLIB) $(GISLIB) $(VECTORLIB)
+DEPENDENCIES = $(RASTER3DDEP) $(GISDEP) $(VECTORDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/grass7/raster3d/r3.what/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.what/main.c
===================================================================
--- grass-addons/grass7/raster3d/r3.what/main.c (rev 0)
+++ grass-addons/grass7/raster3d/r3.what/main.c 2016-07-01 02:27:24 UTC (rev 68823)
@@ -0,0 +1,292 @@
+
+/****************************************************************************
+ *
+ * MODULE: r3.what
+ *
+ * AUTHOR(S): Anna Petrasova (kratochanna gmail.com)
+ *
+ * PURPOSE: Queries 3D raster at specified 2D or 3D coordinates
+ *
+ * COPYRIGHT: (C) 2016 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.
+ *
+ *****************************************************************************/
+#include <grass/gis.h>
+#include <grass/raster3d.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+
+FILE *openAscii(char *file)
+{
+ FILE *fp;
+
+ if (file && file[0] != '-') {
+ fp = fopen(file, "w");
+ if (fp == NULL) {
+ perror(file);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ }
+ else
+ fp = stdout;
+
+ return fp;
+}
+
+
+void query(RASTER3D_Map * input_map, double east, double north, FILE * fp,
+ char *fs, RASTER3D_Region * region, int type, char *null_val)
+{
+
+ int x, y, depth;
+ DCELL dvalue;
+ FCELL fvalue;
+
+ Rast3d_location2coord(region, north, east, region->top, &x, &y, &depth);
+ fprintf(fp, "%f%s%f", east, fs, north);
+ for (depth = 0; depth < region->depths; depth++) {
+ if (type == FCELL_TYPE) {
+ Rast3d_get_value(input_map, x, y, depth, &fvalue, FCELL_TYPE);
+ if (Rast3d_is_null_value_num(&fvalue, FCELL_TYPE))
+ fprintf(fp, "%s%s", fs, null_val);
+ else
+ fprintf(fp, "%s%f", fs, fvalue);
+ }
+ else {
+ Rast3d_get_value(input_map, x, y, depth, &dvalue, DCELL_TYPE);
+ if (Rast3d_is_null_value_num(&dvalue, DCELL_TYPE))
+ fprintf(fp, "%s%s", fs, null_val);
+ else
+ fprintf(fp, "%s%f", fs, dvalue);
+ }
+ }
+ fprintf(fp, "\n");
+}
+
+void query3D(RASTER3D_Map * input_map, double east, double north, double top,
+ FILE * fp, char *fs, RASTER3D_Region * region, int type,
+ char *null_val)
+{
+
+ int x, y, depth;
+ DCELL dvalue;
+ FCELL fvalue;
+
+ Rast3d_location2coord(region, north, east, top, &x, &y, &depth);
+ fprintf(fp, "%f%s%f%s%f", east, fs, north, fs, top);
+ if (type == FCELL_TYPE) {
+ Rast3d_get_value(input_map, x, y, depth, &fvalue, FCELL_TYPE);
+ if (Rast3d_is_null_value_num(&fvalue, FCELL_TYPE))
+ fprintf(fp, "%s%s", fs, null_val);
+ else
+ fprintf(fp, "%s%f", fs, fvalue);
+ }
+ else {
+ Rast3d_get_value(input_map, x, y, depth, &dvalue, DCELL_TYPE);
+ if (Rast3d_is_null_value_num(&dvalue, DCELL_TYPE))
+ fprintf(fp, "%s%s", fs, null_val);
+ else
+ fprintf(fp, "%s%f", fs, dvalue);
+ }
+ fprintf(fp, "\n");
+}
+
+
+int main(int argc, char *argv[])
+{
+ struct Option *input, *output, *points_opt, *coords_opt, *coords3d_opt,
+ *fs_opt, *null_val;
+ struct Flag *mask, *z_flag;
+ struct GModule *module;
+ FILE *fp;
+ RASTER3D_Region region;
+ RASTER3D_Map *input_map;
+ struct Map_info Map;
+ struct line_pnts *Points;
+ int changemask;
+ int i;
+ double east, north, top;
+ int type, ltype;
+ char *fs;
+ int z;
+
+ /* Initialize GRASS */
+ G_gisinit(argv[0]);
+ module = G_define_module();
+ G_add_keyword(_("raster3d"));
+ G_add_keyword(_("query"));
+ G_add_keyword(_("voxel"));
+ module->description =
+ _("Queries 3D raster in specified 2D or 3D coordinates.");
+
+ input = G_define_standard_option(G_OPT_R3_INPUT);
+ input->key = "input";
+ input->description = _("Input 3D raster map");
+ input->guisection = _("Query");
+
+ coords_opt = G_define_standard_option(G_OPT_M_COORDS);
+ coords_opt->required = NO;
+ coords_opt->multiple = YES;
+ coords_opt->description = _("Query with 2D coordinates");
+ coords_opt->guisection = _("Query");
+
+ coords3d_opt = G_define_option();
+ coords3d_opt->key = "coordinates_3d";
+ coords3d_opt->type = TYPE_DOUBLE;
+ coords3d_opt->required = NO;
+ coords3d_opt->multiple = YES;
+ coords3d_opt->key_desc = "east,north,top";
+ coords3d_opt->gisprompt = "old,coords,coords";
+ coords3d_opt->description = _("Query with 3D coordinates");
+ coords3d_opt->guisection = _("Query");
+
+ points_opt = G_define_standard_option(G_OPT_V_MAP);
+ points_opt->key = "points";
+ points_opt->label = _("Name of vector points map for query");
+ points_opt->required = NO;
+ points_opt->guisection = _("Query");
+
+ output = G_define_standard_option(G_OPT_F_OUTPUT);
+ output->required = NO;
+ output->description = _("Name for output file");
+ output->guisection = _("Print");
+
+ fs_opt = G_define_standard_option(G_OPT_F_SEP);
+ fs_opt->guisection = _("Print");
+
+ null_val = G_define_standard_option(G_OPT_M_NULL_VALUE);
+ null_val->answer = "*";
+ null_val->guisection = _("Print");
+
+ mask = G_define_flag();
+ mask->key = 'm';
+ mask->description = _("Use 3D raster mask (if exists) with input map");
+
+ z_flag = G_define_flag();
+ z_flag->key = 'z';
+ z_flag->description = _("Ignore points Z values");
+
+ G_option_required(coords_opt, coords3d_opt, points_opt, NULL);
+ G_option_exclusive(coords_opt, points_opt, NULL);
+ G_option_exclusive(coords3d_opt, coords_opt, NULL);
+
+ /* Have GRASS get inputs */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ if (!G_find_raster3d(input->answer, ""))
+ Rast3d_fatal_error(_("3D raster map <%s> not found"), input->answer);
+
+ /* Initiate the default settings */
+ Rast3d_init_defaults();
+
+ /* Figure out the current region settings */
+ Rast3d_get_window(®ion);
+
+ /* Open the map and use XY cache mode */
+ input_map =
+ Rast3d_open_cell_old(input->answer,
+ G_find_raster3d(input->answer, ""), ®ion,
+ RASTER3D_TILE_SAME_AS_FILE,
+ RASTER3D_USE_CACHE_DEFAULT);
+
+ if (!input_map)
+ Rast3d_fatal_error(_("Unable to open 3D raster map <%s>"),
+ input->answer);
+
+ type = Rast3d_tile_type_map(input_map);
+
+ /* Open the output ascii file */
+ fp = openAscii(output->answer);
+
+ fs = G_option_to_separator(fs_opt);
+
+ z = z_flag->answer;
+
+ /*if requested set the mask on */
+ if (mask->answer) {
+ if (Rast3d_mask_file_exists()) {
+ changemask = 0;
+ if (Rast3d_mask_is_off(input_map)) {
+ Rast3d_mask_on(input_map);
+ changemask = 1;
+ }
+ }
+ }
+ /* open vector points map */
+ if (points_opt->answer) {
+ Vect_set_open_level(1); /* topology not required */
+ if (Vect_open_old(&Map, points_opt->answer, "") < 0)
+ G_fatal_error(_("Unable to open vector map <%s>"),
+ points_opt->answer);
+ Points = Vect_new_line_struct();
+ while (1) {
+ ltype = Vect_read_next_line(&Map, Points, NULL);
+ if (ltype == -1)
+ G_fatal_error(_("Unable to read vector map <%s>"),
+ Vect_get_full_name(&Map));
+ else if (ltype == -2)
+ break;
+ else if (!(ltype & GV_POINTS)) {
+ G_warning(_("Line is not point or centroid, skipped"));
+ }
+ else {
+ east = Points->x[0];
+ north = Points->y[0];
+ if (Vect_is_3d(&Map) && !z) {
+ top = Points->z[0];
+ query3D(input_map, east, north, top, fp, fs, ®ion,
+ type, null_val->answer);
+ }
+ else
+ query(input_map, east, north, fp, fs, ®ion, type,
+ null_val->answer);
+ }
+ }
+ }
+ else {
+ /* loop through coordinates */
+ if (coords3d_opt->answer) {
+ for (i = 0; coords3d_opt->answers[i] != NULL; i += 3) {
+ G_scan_easting(coords3d_opt->answers[i], &east,
+ G_projection());
+ G_scan_northing(coords3d_opt->answers[i + 1], &north,
+ G_projection());
+ top = atof(coords3d_opt->answers[i + 2]);
+ query3D(input_map, east, north, top, fp, fs, ®ion, type,
+ null_val->answer);
+ }
+ }
+ else {
+ for (i = 0; coords_opt->answers[i] != NULL; i += 2) {
+ G_scan_easting(coords_opt->answers[i], &east, G_projection());
+ G_scan_northing(coords_opt->answers[i + 1], &north,
+ G_projection());
+ query(input_map, east, north, fp, fs, ®ion, type,
+ null_val->answer);
+ }
+ }
+ }
+
+ /* We set the mask off, if it was off before */
+ if (mask->answer) {
+ if (Rast3d_mask_file_exists())
+ if (Rast3d_mask_is_on(input_map) && changemask)
+ Rast3d_mask_off(input_map);
+ }
+
+ /* Close files and exit */
+ if (!Rast3d_close(input_map))
+ Rast3d_fatal_error(_("Unable to close 3D raster map"));
+
+ if (output)
+ if (fclose(fp))
+ Rast3d_fatal_error(_("Unable to close new ASCII file"));
+
+ return 0;
+}
Property changes on: grass-addons/grass7/raster3d/r3.what/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster3d/r3.what/r3.what.html
===================================================================
--- grass-addons/grass7/raster3d/r3.what/r3.what.html (rev 0)
+++ grass-addons/grass7/raster3d/r3.what/r3.what.html 2016-07-01 02:27:24 UTC (rev 68823)
@@ -0,0 +1,58 @@
+<h2>DESCRIPTION</h2>
+
+<em>r3.what</em> queries 3D raster at specified coordinates. The <em>input</em>
+parameter is a valid 3D raster map in the current mapset search path.
+The <em>output</em> parameter is used to output results in a text file.
+If <em>output</em> is not specified or '-' is used, then standard output
+ (stdout) is used.
+
+<p>
+Coordinates can be specified directly using options <b>coordinates</b>
+as coordinate pairs (east, north), or using option <b>coordinates_3d</b>
+as east, north and top. If the coordinates are 2D, the result will be
+a list of values representing the vertical column of the 3D raster at that
+coordinates (from bottom to top). If 3D coordinates are provided only one value
+is returned. User can specify multiple coordinates.
+<p>
+Another option is to specify
+vector points map using option <b>points</b>. It can be either 2D or 3D.
+If a 3D vector points map is specified, one 3D raster value per point is printed.
+If that is not desired, flag <b>z</b> will ignore the z coordinate and
+the module returns values for the vertical column.
+
+<h2>Example</h2>
+We create a 3D raster where values depend on the depth.
+We query the 3D raster with two 2D coordinates:
+<div class="code"><pre>
+g.region s=0 n=100 w=0 e=100 b=0 t=100 res3=1 res=1 tbres=10 -p3
+r3.mapcalc "new = depth()"
+r3.what input=new at user1 coordinates=10,20,50,50
+</pre></div>
+
+<pre>
+10.000000|20.000000|1.000000|2.000000|3.000000|4.000000|5.000000|6.000000|7.000000|8.000000|9.000000|10.000000
+50.000000|50.000000|1.000000|2.000000|3.000000|4.000000|5.000000|6.000000|7.000000|8.000000|9.000000|10.000000
+</pre>
+Here we query the 3D raster with two 3D coordinates:
+<div class="code"><pre>
+r3.what input=new at user1 coordinates_3d=10,10,1,10,10,11
+</pre></div>
+
+<pre>
+10.000000|10.000000|1.000000|1.000000
+10.000000|10.000000|11.000000|2.000000
+</pre>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.what.html">r.what</a>,
+<a href="r3.out.ascii.html">r3.out.ascii</a>,
+<a href="v.what.rast3.html">v.what.rast3</a>,
+<a href="g.region.html">g.region</a>
+</em>
+
+<h2>AUTHORS</h2>
+Anna Petrasova, NCSU GeoForAll Lab
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster3d/r3.what/r3.what.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list