1
0
Fork 0

Initial checkin of new GRASS modules

This commit is contained in:
Ralf Gerlich 2010-05-06 09:04:51 +02:00
parent 1ad2b2e0e9
commit 74be0ac8b0
7 changed files with 604 additions and 0 deletions

3
gismodules/v.contri/.gitignore vendored Normal file
View file

@ -0,0 +1,3 @@
triangle.c
triangle.h
triangle.zip

297
gismodules/v.contri/main.c Normal file
View file

@ -0,0 +1,297 @@
/****************************************************************
*
* MODULE: v.contri
*
* AUTHOR(S): Ralf Gerlich
*
* PURPOSE: Create a conforming delauney triangulation from a
* vector map
*
* COPYRIGHT: (C) 2010 by Ralf Gerlich
*
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*
****************************************************************/
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/vector.h>
/* #define SINGLE */
#ifdef SINGLE
#define REAL float
#else /* not SINGLE */
#define REAL double
#endif /* not SINGLE */
#include "triangle.h"
void build_segments(struct triangulateio* tri, struct Map_info* map)
{
plus_t node_idx, line_idx;
const plus_t node_num = Vect_get_num_nodes(map);
const plus_t line_num = Vect_get_num_lines(map);
struct line_pnts *pnts = Vect_new_line_struct();
int index, point_idx, segment_idx;
int n1, n2, a1, a2;
/* Collect points and segments */
/* The nodes are shared between edges, but the other points are
* individual, so we place the nodes at the start for easy indexing
* and add the other points lazily.
*
* FIXME: actually, centroids are also nodes, but should not be
* included as vertices.
*/
tri->numberofpoints = node_num;
tri->numberofsegments = 0;
tri->numberofpointattributes = 1; // We store the z-coordinate there
for ( line_idx = 1; line_idx <= line_num; line_idx++ ) {
Vect_reset_line( pnts );
const int type = Vect_read_line( map, pnts, NULL, line_idx );
if ( type == GV_POINT ) {
tri->numberofpoints++;
} else if ( type == GV_BOUNDARY ) {
tri->numberofpoints+=pnts->n_points-2;
tri->numberofsegments+=pnts->n_points-1;
}
}
tri->pointlist = (REAL*) malloc( sizeof(REAL[2])*tri->numberofpoints );
tri->pointattributelist = (REAL*) malloc( sizeof(REAL)*tri->numberofpoints );
tri->segmentlist = (int*) malloc( sizeof(int[2])*tri->numberofsegments );
tri->segmentmarkerlist = (int*) malloc( sizeof(int)*tri->numberofsegments );
/* Add the nodes first */
for ( node_idx = 1; node_idx <= node_num; node_idx++ ) {
Vect_get_node_coor( map, node_idx,
&tri->pointlist[ 2*node_idx - 2 ],
&tri->pointlist[ 2*node_idx - 1 ],
&tri->pointattributelist[ node_idx - 1] );
}
point_idx = node_num;
segment_idx = 0;
for ( line_idx = 1; line_idx <= line_num; line_idx++ ) {
Vect_reset_line( pnts );
const int type = Vect_read_line( map, pnts, NULL, line_idx );
if ( type == GV_POINT ) {
tri->pointlist[ 2*point_idx + 0 ] = pnts->x[0];
tri->pointlist[ 2*point_idx + 1 ] = pnts->y[0];
tri->pointattributelist[ point_idx ] = pnts->z[0];
point_idx++;
} else if ( type == GV_BOUNDARY ) {
Vect_get_line_nodes( map, line_idx, &n1, &n2 );
Vect_get_line_areas( map, line_idx, &a1, &a2 );
int last_index = n1-1;
int segmarker = (a1==0 || a2==0 ? 1 : 0);
for ( index = 1; index < pnts->n_points-1; index++ ) {
tri->pointlist[ 2*point_idx + 0 ] = pnts->x[index];
tri->pointlist[ 2*point_idx + 1 ] = pnts->y[index];
tri->pointattributelist[ point_idx ] = pnts->y[index];
tri->segmentlist[ 2 * segment_idx + 0 ] = last_index;
tri->segmentlist[ 2 * segment_idx + 1 ] = point_idx;
tri->segmentmarkerlist[ segment_idx ] = segmarker;
segment_idx++;
last_index = point_idx;
point_idx++;
}
tri->segmentlist[ 2 * segment_idx + 0 ] = last_index;
tri->segmentlist[ 2 * segment_idx + 1 ] = n2-1;
tri->segmentmarkerlist[ segment_idx ] = segmarker;
segment_idx++;
}
}
Vect_destroy_line_struct( pnts );
}
void build_regions(struct triangulateio* tri, struct Map_info* map)
{
plus_t area_idx;
const plus_t area_num = Vect_get_num_areas(map);
struct line_pnts *pnts = Vect_new_line_struct();
int region_idx, hole_idx;
tri->numberofregions = 0;
tri->numberofholes = 0;
for ( area_idx = 0; area_idx <= area_num; area_idx++ ) {
if ( !Vect_area_alive( map, area_idx) ) {
continue;
}
if ( Vect_get_area_centroid( map, area_idx ) > 0 ) {
tri->numberofregions++;
} else {
tri->numberofholes++;
}
}
tri->regionlist = (REAL*) malloc(sizeof(REAL[4]) * tri->numberofregions);
tri->holelist = (REAL*) malloc(sizeof(REAL[2]) * tri->numberofholes);
region_idx = 0;
hole_idx = 0;
for ( area_idx = 0; area_idx <= area_num; area_idx++ ) {
if ( !Vect_area_alive( map, area_idx) ) {
continue;
}
const int centroid_idx = Vect_get_area_centroid( map, area_idx );
if ( centroid_idx > 0 ) {
Vect_reset_line( pnts );
Vect_read_line( map, pnts, NULL, centroid_idx );
tri->regionlist[ 4*region_idx + 0 ] = pnts->x[0];
tri->regionlist[ 4*region_idx + 1 ] = pnts->y[0];
tri->regionlist[ 4*region_idx + 2 ] = centroid_idx;
tri->regionlist[ 4*region_idx + 3 ] = -1.0; // maximum area, unused
region_idx++;
} else {
Vect_get_point_in_area( map,
area_idx,
&tri->holelist[ 2*hole_idx + 0],
&tri->holelist[ 2*hole_idx + 1] );
hole_idx++;
}
}
Vect_destroy_line_struct( pnts );
}
void build_triangulateio(struct triangulateio* tri, struct Map_info* map)
{
build_segments(tri, map);
build_regions(tri, map);
}
void build_outputvector(struct Map_info* newmap,
struct triangulateio* tri,
struct Map_info* oldmap)
{
int index;
struct line_pnts* points = Vect_new_line_struct();
struct line_cats* cats = Vect_new_cats_struct();
/* Write out the segments */
for ( index = 0; index < tri->numberofedges; index++ ) {
int * const edge = &tri->edgelist[ 2 * index ];
Vect_reset_line( points );
Vect_append_point( points,
tri->pointlist[ 2 * edge[0] + 0 ],
tri->pointlist[ 2 * edge[0] + 1 ],
tri->pointattributelist[ edge[0] ]);
Vect_append_point( points,
tri->pointlist[ 2 * edge[1] + 0 ],
tri->pointlist[ 2 * edge[1] + 1 ],
tri->pointattributelist[ edge[1] ]);
Vect_write_line( newmap, GV_BOUNDARY, points, cats );
}
/* Write out the region markers */
for ( index = 0; index < tri->numberoftriangles; index++ ) {
int * const triele = &tri->trianglelist[ index * tri->numberofcorners ];
REAL x, y, z;
x = ( tri->pointlist[ 2 * triele[0] + 0 ] +
tri->pointlist[ 2 * triele[1] + 0 ] +
tri->pointlist[ 2 * triele[2] + 0 ] ) / 3.0;
y = ( tri->pointlist[ 2 * triele[0] + 1 ] +
tri->pointlist[ 2 * triele[1] + 1 ] +
tri->pointlist[ 2 * triele[2] + 1 ] ) / 3.0;
z = ( tri->pointattributelist[ triele[0] ] +
tri->pointattributelist[ triele[1] ] +
tri->pointattributelist[ triele[2] ] ) / 3.0;
Vect_reset_line( points );
Vect_reset_cats( cats );
Vect_read_line( oldmap,
NULL,
cats,
(plus_t) tri->triangleattributelist[ index * tri->numberoftriangleattributes ]);
Vect_append_point( points, x, y, z );
Vect_write_line( newmap, GV_CENTROID, points, cats );
}
}
void build_outputholes(struct Map_info* newmap,
struct triangulateio* tri)
{
int index;
struct line_pnts* points = Vect_new_line_struct();
struct line_cats* cats = Vect_new_cats_struct();
/* Write out the hole markers */
for ( index = 0; index < tri->numberofholes; index++ ) {
REAL * const hole = &tri->holelist[ 2 * index ];
Vect_reset_line( points );
Vect_append_point( points, hole[0], hole[1], 0.0 );
Vect_write_line( newmap, GV_POINT, points, cats );
}
}
int main(int argc, char *argv[])
{
struct GModule *module; /* GRASS module for parsing arguments */
struct Option *old, *new;
struct Map_info oldmap, newmap;
G_gisinit(argv[0]);
/* initialize module */
module = G_define_module();
module->description = _("Create a conforming delauney triangulation from a vector");
/* Define the different options as defined in gis.h */
old = G_define_standard_option(G_OPT_V_INPUT);
new = G_define_standard_option(G_OPT_V_OUTPUT);
/* options and flags parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
if ( Vect_open_old( &oldmap, old->answer, NULL ) < 2 ) {
G_fatal_error("Unable to open vector map \"%s\"",
old->answer);
}
if ( Vect_open_new( &newmap, new->answer, Vect_is_3d(&oldmap) )!=1 ) {
G_fatal_error("Unable to create vector map \"%s\"",
new->answer);
}
Vect_copy_head_data( &oldmap, &newmap );
Vect_copy_tables( &oldmap, &newmap, 0 );
struct triangulateio in, out;
build_triangulateio(&in, &oldmap);
out.pointlist = (REAL*) NULL;
out.pointattributelist = (REAL*) NULL;
out.pointmarkerlist = (int*) NULL;
out.trianglelist = (int*) NULL;
out.triangleattributelist = (REAL*) NULL;
out.neighborlist = (int*) NULL;
out.segmentlist = (int*) NULL;
out.segmentmarkerlist = (int*) NULL;
out.edgelist = (int*) NULL;
out.edgemarkerlist = (int*) NULL;
triangulate("pzAejC", &in, &out, NULL);
build_outputvector(&newmap, &out, &oldmap);
//build_outputholes(&newmap, &in);
Vect_build(&newmap);
Vect_close(&oldmap);
Vect_close(&newmap);
/* Don't forget to report to caller sucessful end of data processing :) */
exit(EXIT_SUCCESS);
}

View file

@ -0,0 +1,21 @@
<h2>DESCRIPTION</h2>
<em>v.example</em> is an example vector module that does something like
labeling all vectors with value 1. A new map is written instead of updating
the input map.
<h2>EXAMPLE</h2>
<div class="code"><pre>
v.example input=map output=newmap
</pre></div>
<H2>SEE ALSO</H2>
<EM><A HREF="http://grass.itc.it/devel/index.php#prog">GRASS Programmer's Manual</A></EM>
<h2>AUTHOR</h2>
Radim Blazek, ITC-irst, Trento, Italy
<p><i>Last changed: $Date: 2008-08-15 08:16:42 +0200 (Fr, 15. Aug 2008) $</i>

182
gismodules/v.flatten/main.c Normal file
View file

@ -0,0 +1,182 @@
/****************************************************************
*
* MODULE: v.contri
*
* AUTHOR(S): Ralf Gerlich
*
* PURPOSE: Flatten a vector map according to several operators
*
* COPYRIGHT: (C) 2010 by Ralf Gerlich
*
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*
****************************************************************/
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/dbmi.h>
#include <grass/vector.h>
void apply_setz_operator( struct Map_info* map, double z, int* cats, int ncats, int layer)
{
const plus_t area_count = Vect_get_num_areas( map );
plus_t area_idx;
struct line_cats *lcats = Vect_new_cats_struct();
for (area_idx = 0; area_idx <= area_count; area_idx++) {
if ( !Vect_area_alive( map, area_idx ) ) {
/* Skip dead areas */
continue;
}
if ( ncats >= 0 ) {
// TODO: check whether the area has the given categories
}
// TODO: apply setz
}
}
void apply_mean_operator( struct Map_info* map, int* cats, int ncats, int layer )
{
if (ncats > 0) {
/* We have to greedily collect adjacent selected areas */
// TODO
} else if (ncats < 0) {
/* We have to greedily collect adjacent areas */
// TODO
}
}
void apply_slope_operator( struct Map_info* map, double slope, int* cats, int ncats, int layer)
{
if (ncats > 0) {
/* We apply slope to individual areas */
// TODO
} else if (ncats < 0) {
/* We apply slope to individual areas */
// TODO
}
}
int main(int argc, char *argv[])
{
struct GModule *module; /* GRASS module for parsing arguments */
struct Option *old, *new;
struct Option *whereopt;
struct Option *operator, *value, *slope;
struct Option *fieldopt;
int field;
struct Map_info oldmap, newmap;
struct field_info* fieldinfo;
dbDriver* driver;
int i, ncats, *cats;
G_gisinit(argv[0]);
/* initialize module */
module = G_define_module();
module->description = _("Create a conforming delauney triangulation from a vector");
/* Define the different options as defined in gis.h */
old = G_define_standard_option(G_OPT_V_INPUT);
whereopt = G_define_standard_option(G_OPT_DB_WHERE);
fieldopt = G_define_standard_option(G_OPT_V_FIELD);
new = G_define_standard_option(G_OPT_V_OUTPUT);
operator = G_define_option();
operator->key = "operator";
operator->description = _("Flattening operator to apply");
operator->type = TYPE_STRING;
operator->required = YES;
operator->options = "setz,mean,slope";
value = G_define_option();
value->key = "value";
value->description = _("The target value for the 'setz' operator");
value->type = TYPE_DOUBLE;
value->required = NO;
value->answer = "0";
slope = G_define_option();
slope->key = "slope";
slope->description = _("The slope value for the 'slope' operator");
slope->type = TYPE_DOUBLE;
slope->required = NO;
slope->answer = "0.4";
/* options and flags parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
if ( Vect_open_old( &oldmap, old->answer, NULL ) < 2 ) {
G_fatal_error("Unable to open vector map \"%s\"",
old->answer);
}
if ( Vect_open_new( &newmap, new->answer, Vect_is_3d(&oldmap) )!=1 ) {
G_fatal_error("Unable to create vector map \"%s\"",
new->answer);
}
Vect_copy_head_data( &oldmap, &newmap );
Vect_copy_tables( &oldmap, &newmap, 0 );
Vect_copy_map_lines( &oldmap, &newmap );
Vect_close(&oldmap);
if (whereopt->answer!=NULL) {
field = Vect_get_field_number( &newmap, fieldopt->answer );
fieldinfo = Vect_get_field( &newmap, field );
if (!fieldinfo) {
G_fatal_error(_("Database connection not defined for layer <%s>"),
fieldopt->answer);
}
G_debug(1, "Loading categories from table <%s>", fieldinfo->table);
driver = db_start_driver_open_database(fieldinfo->driver, fieldinfo->database);
if (driver == NULL) {
G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
fieldinfo->database, fieldinfo->driver);
}
ncats = db_select_int(driver, fieldinfo->table, fieldinfo->key, whereopt->answer,
&cats);
if (ncats == -1) {
G_fatal_error(_("Unable select records from table <%s>"), fieldinfo->table);
}
G_message(_("%d categories loaded from table <%s>"), ncats,
fieldinfo->table);
db_close_database(driver);
db_shutdown_driver(driver);
} else {
field = -1;
cats = NULL;
ncats = -1;
}
if (!strcmp(operator->answer, "setz")) {
apply_setz_operator( &newmap, atof(value->answer), cats, ncats, field );
} else if (!strcmp(operator->answer, "mean")) {
apply_mean_operator( &newmap, cats, ncats, field );
} else if (!strcmp(operator->answer, "slope")) {
apply_slope_operator( &newmap, atof(slope->answer), cats, ncats, field );
}
Vect_build(&newmap);
if (cats != NULL) {
G_free(cats);
}
Vect_close(&newmap);
/* Don't forget to report to caller sucessful end of data processing :) */
exit(EXIT_SUCCESS);
}

View file

@ -0,0 +1,21 @@
<h2>DESCRIPTION</h2>
<em>v.example</em> is an example vector module that does something like
labeling all vectors with value 1. A new map is written instead of updating
the input map.
<h2>EXAMPLE</h2>
<div class="code"><pre>
v.example input=map output=newmap
</pre></div>
<H2>SEE ALSO</H2>
<EM><A HREF="http://grass.itc.it/devel/index.php#prog">GRASS Programmer's Manual</A></EM>
<h2>AUTHOR</h2>
Radim Blazek, ITC-irst, Trento, Italy
<p><i>Last changed: $Date: 2008-08-15 08:16:42 +0200 (Fr, 15. Aug 2008) $</i>

View file

@ -0,0 +1,59 @@
/****************************************************************
*
* MODULE: v.contri
*
* AUTHOR(S): Ralf Gerlich
*
* PURPOSE: Export a TIN as FlightGear .btg.gz-file
*
* COPYRIGHT: (C) 2010 by Ralf Gerlich
*
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*
****************************************************************/
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/vector.h>
int main(int argc, char *argv[])
{
struct GModule *module; /* GRASS module for parsing arguments */
struct Option *old, *new;
struct Map_info oldmap;
G_gisinit(argv[0]);
/* initialize module */
module = G_define_module();
module->description = _("Export a TIN as FlightGear .btg.gz");
/* Define the different options as defined in gis.h */
old = G_define_standard_option(G_OPT_V_INPUT);
new = G_define_standard_option(G_OPT_F_OUTPUT);
/* options and flags parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
if ( Vect_open_old( &oldmap, old->answer, NULL ) < 2 ) {
G_fatal_error("Unable to open vector map \"%s\"",
old->answer);
}
// TODO: stripify triangles by categories
// TODO: add singular triangles
// TODO: create texture coordinates
// TODO: write out btg.gz
Vect_close(&oldmap);
/* Don't forget to report to caller sucessful end of data processing :) */
exit(EXIT_SUCCESS);
}

View file

@ -0,0 +1,21 @@
<h2>DESCRIPTION</h2>
<em>v.example</em> is an example vector module that does something like
labeling all vectors with value 1. A new map is written instead of updating
the input map.
<h2>EXAMPLE</h2>
<div class="code"><pre>
v.example input=map output=newmap
</pre></div>
<H2>SEE ALSO</H2>
<EM><A HREF="http://grass.itc.it/devel/index.php#prog">GRASS Programmer's Manual</A></EM>
<h2>AUTHOR</h2>
Radim Blazek, ITC-irst, Trento, Italy
<p><i>Last changed: $Date: 2008-08-15 08:16:42 +0200 (Fr, 15. Aug 2008) $</i>