GRASS Programming Howto
- temporary solution here, eventually it becomes a doygenized page in SVN.
- See also GRASS 6 Programmer's manual
Introduction
The reason for starting this Howto is the fact that it's very hard to get a good overview on the Grass code base, and documentation is sparse, lacking or not clear on many places. The aim of this Howto therefore, is to do something about that, so it will be easier to contribute to, and better document the Grass code base. I encourage everybody to submit changes, correct me where I'm wrong and add information on the places where I haven't been able to do so. For a start I've chosen to follow the outline of the Grass Programmer's manual with some changes where needed. So let's start (26-01-2003) and make it more fun/ease to program for Grass.
Greetings
J.c.m van der Kwast (aka Sjors) Bsc Tropical Agriculture, Plant Production. The Netherlands
Tools
If you are new to programming I suggest that you start looking at a few tutorials on c, which are widely available on the internet. Also try to find somebody who has some experience and get them to help you on you're way (I luckily had a friend who could help me get started, thanks C.W. de Gier). Also have a look at the Programmer's manuals (link: PDF 60/61 http://grass.itc.it/grass60/manuals/).
To get an impression of how the program/modules work there's a need to help you see which functions are called when, from which place and what it is they pass etc. To do this you'll have to us a debugger like DDD and GDB, here's a link to some debugging hints. This could be a separate chapter with some screen shots. We'll see how to do this.
For ease I also make a distinction in the sources for Grass5.03/53 which doesn't have debug info coded in, and Grass6 which has. To set it up for grass60 you have to give this command in the Grass shell:
g.gisenv set=DEBUG=#number
where #number stands for the level. There are three levels, 1, 3, 5. Play around with these to see the differences in output. If you want output directed to a file it's necessary to edit you're .bashrc or .bash_profile and place a line like this:
export GRASS_DEBUG__FILE=/$Path/debug.txt or any name you like
You have to logout and login again to make it work. Then create a an empty file with vi or a command like:
$touch /$Path/debug.txt
To view some of the ouput while running Grass57 you can start a new terminal or xterm and use this command:
$tail -f /$Path/debug.txt
This file can get very big fast, so remove it once in a while.
Emacs
I use Emacs as my main editor and viewer of the code. To get a overview, it's best to begin with making a TAGS file, which will make an overview of keywords (eg. function names, defines) and point to the places where it's used. To set this up use a command like this in the main directory of the Grass source code:
$/Path_to_grass_source_code/find . -name '*.[ch]' -print |etags -
Now you can start Emacs with "$emacs &" from you're shell. For a
complete reference read the Emacs manual, here are some of the
commands to search the TAGS file. I will use the Emacs notation,
mainly C-, push ctrl key, M- push meta or alt key.
For DEFINES and Functions, this will point to the code of the DEFINE or function:
'M-.' then type search string etc
For references of functions, this will point to any file where a function is used:
'M-x tags-search' then type search string etc
Now get handy with the tools and notice how incredibly complex the Grass sources are. This however is you're main source of learning, viewing and understanding what already has been coded.
Once you have the code (new one or fixed one)
Now if you decide to add or change something do this according to the rules explained in SUBMITTING.txt (link: http://www.geog.uni-hannover.de/grass/grass5/source/SUBMITTING) and take a look at "GNU Coding Standards" (link: http://www.gnu.org/prep/standards.html). Some other hints and needs for the Grass programmer:
- Document you're new code well, and in browsing other code and gaining valuable insights, also document this either here, but especially in the code itself.
- As there is no DEBUG info code for the main source we must consider this, like has been done in Grass57, and use this abundantly in new code.
- As I understand the idea is to use doxygen to extract the comments and documentation out of the source. So learn how to use this when submitting you're findings or changes.
Existing API Documentation
In order to get quick documentation on functions structures while programming there are tools designed for extracting documentation directly from the source code so that documentation may be auto-generated at every update.
GRASS6 html documentation is produced by:
make htmldocs
or
make pdfdocs
this will create the API documentation in the doxygenhtml directory just below the GRASS6 source tree.
GIS programming in General
In GIS, from a user viewpoint, there's two main aspects in general to deal with. The geographical aspect (coordinates/location) and data we want tie to this position, in order to calculate, combine, interpolate (data and/or coordinates), or whatever it is you want to do with these two aspects. Then there are two main ways of presentation of these aspects, in the form of raster and vector maps. In Grass there's an addition to that in the form of sites. In my understanding this will be changed to either a cell file for raster (already a fact), a definition of point structure for vectors (in Grass57). So in the end this will be not be there anymore as a separate structure.
In Grass it's known that now and in the past the structure of these two aspects were always somewhat confusing and somewhat limited. Mainly we had for both raster and vector only two possibilities to tie data to a location, being the cat and label. Externally there's always been the database connection, direct or via odbc. Now with the new vector capabilities of grass57, and the change in structure, there's the possibility to add more data then cat and label to a map. Everybody being very excited with these new features, it's to be expected that these will have to be implemented for raster maps as well.
Now for the programming part, we'll all have to think about how we can provide a structure which will ease all the tasks needed to do this. In general we should consider most probably that we have the main Grass program environment, which will be responsible for things like display, i/o handling etc, and modules for specific tasks (eg digitize new map).
In the Grass sources there's a great need for getting more structure, better defined api's, in general more according to proper coding standards. If you have followed the discussions lately on this subject on the developer's mailing list (November 2003/Jan 2004), I'm not the only one (also my friend helping me out to get started) who noticed the lack of consistency and proper coding standards in some parts of the Grass sources. Also a good overview of how all these modules and drivers are connected/interacting is missing. That's one of the reason's of course for starting on this Programming Howto. When finished a diagram will be placed here to explain this. later more on this.
Libraries
To gain a insight here look at the directory in $Path_to_grass_source_code/lib/, where all the library code is placed. A division is made according to the several distinctions in tasks these libraries perform. Here's a list of the main libraries so far:
GIS library (GisLib) Raster library (RasterLib) Vector library (VectLib) Imagery library (ImgLib) Raster Graphics Library (RstGraphicsLib) Display Graphics Library (DisplGraphicsLib) Lock Library (LockLib) Rowio Library (RowioLib) Segment Library (SegLib) Vask Library (VaskLib) Projection and Datum Libraries (ProjLib) Grid3D Library (Grd3dLib) Date/time Library (DateTimeLib) OpenGL gsurf Library (OpenglLib) Numerical math lib with optional interface to LAPACK/BLAS (GmathLib) Bitmap Library (BitmapLib)
A few things of which I think are not libraries but could be are the following:
GUI programming Digitizer/Mouse/Trackball
The main idea is to provide good documentation for functions in each library in the next paragraphs, and add usable examples of how to use and implement functions. Basically there's a need for a good overview, creating and/or explaining API's etc.
GIS library
Info and functions explained will come here: GisLib
Header files
System header files List (At least for Grass57)
bdlg_bm.h bitmap.h blas.h btree.h CC.h codes.h colors.h config.h.in datetime.h dbmi.h devlib.h D.h dig_atts.h display.h dlg_bm.h dlg.h edit.h form.h G3d.h geo.h geom.h gisdefs.h Gis.H glocale.h gmath.h gproj_api.h gprojects.h graph.h help.h ibtree.h icon_bm.h icon.h imagedefs.h imagery.h label_bm.h la.h lapack.h libtrans.h linkm.h lock.h monitors.h ortholib.h Paintlib.h patterns.h pbmplus.h P_datetime.h proto_dbmi.h P_site.h raster.h readsites.h region_bm.h rowio.h search.h segment.h shhopt.h site.h sitelib.h sqlp.h std_incs.h symbol.h transform.h vask.h vbuildlib.h VecT.h VERSION version.h.in V_.h winname.h.in
Directory Make
Dir.make Grass.make.in Lib.make Module.make Platform.make.in Rules.make Script.make Shlib.make Stlib.make
Directory include/vect
These header files are used internally only.
dig_defines.h dig_externs.h dig_globs.h digit.h dig_macros.h dig_nodes.h Dig_Structs.h
Raster library
Used for display in the GRASS monitor (the raster map programming is within libgis!).
Vector library
Grass53 Info and functions explained will come here (outdated)
Grass6 Info and functions explained see next
VectLib
Vector library consists in fact of 5 libraries. Here's a list:
DgLib the Directed Graph Library for network analysis DigLib lower level functions Rtree spatial index for lower level functions TransForm transformation Vlib higher level functions, intended for programming modules
Imagery library
Info and functions explained will come here: ImgLib
Raster Graphics Library
Info and functions explained will come here: RstGraphicsLib
Display Graphics Library
Info and functions explained will come here: DisplGraphicsLib
Lock Library
Info and functions explained will come here: LockLib
Rowio Library
Info and functions explained will come here: RowioLib
Segment Library
Info and functions explained will come here: SegLib
Vask Library
Info and functions explained will come here: VaskLib
Projection and Datum Libraries
Info and functions explained will come here: ProjLib
Grid3D Library
Info and functions explained will come here: Grd3dLib
DateTime Library
Info and functions explained will come here: DateTimeLib
OpenGL gsurf Library
Info and functions explained will come here: OpenglLib
GMath numerical Library
Info and functions explained will come here: GMathLib
Bitmap Library
Info and functions explained will come here: BitmapLib
Database Structure
Everything related to internal Grass Database structure
Raster Modules
Everything related to Raster Modules (make list in progress)
GRASS53 and GRASS6
Example of raster module source code and Makefile: http://freegis.org/cgi-bin/viewcvs.cgi/grass6/doc/raster/raster.tar.gz?tarball=1
Raster Module Programming
After running G_parse(), input and output file names should be checked to make sure that the input file exists, the output file is not the input file, and the output filename is legal. The third argument can be one of the following:
- GR_FATAL_EXIT - Prints error and exits program
- GR_FATAL_PRINT - Prints warning and returns error code
- GR_FATAL_RETURN - Returns error code
G_check_input_output_name(input_name, output_name, GR_FATAL_EXIT);
Vector Modules
Everything related to vector modules again separated for Grass53 and Grass6
Beginning to list the modules and exploring the differences between Grass53 and Grass6. Most modules either are replaced or updated by new or updated version in Grass6 or need transformation to Grass6. It would be nice to make a sort of template structure, using the same names for arguments as (partially?) has been done by Radim in Grass6. For this look at my first try at explaining the module v.buffer in Grass6. I think describing a module will also provide necessary info on lib functions, so these will be linked.
Example of vector module source code and Makefile: http://freegis.org/cgi-bin/viewcvs.cgi/grass6/doc/vector/vector.tar.gz?tarball=1
GRASS6
Vector model, API introduction and example module source code is in the Grass6 source tree: grass6x/doc/vector/
Example of vector module source code:
grass64/doc/vector/v.example/
Grass6 vector module list
Vector modules are usually in $grass_source_home/vector/
Vector Module Programming
v.buffer: An Example
Well first try on description of a module in a almost dummy way. As I'm new to programming please try to help me out, I'm bound to make mistakes. I'll assume that you have a basic understanding of c and are familiar with pointers, arrays and structs. As I do often not have a nice picture in my head containing a array of structs, I often do to some extend know what's happening. Maybe in time we could visualize some of this in a nice way. This module has 1 file main.c and a makefile. I'll begin with main.c and the function main as this seems the most logical place to start. In Grass modules almost behave like a standalone program as they have standard arguments found in c programs in general. The main thing where there is a difference to a standalone program is the Grass parser. Let's take a look at standard command line arguments.
v.buffer help
Prints help.
v.buffer input=map output=buffer type=line buffer=100 tolerance=1 debug=clean
Arguments debug and tolerance are a bit unclear but the rest is self explanatory, this draws a 100 map units buffer around lines in "map", output it to "buffer" with a tolerance of 1 (default).
Now lets take a look at function main() in main.c which has to receive this and do something with it.
int Main(int argc, char *argv[]); /* with argc number of arguments including
program name and argv a pointer to a
arguments array, normal C form */
Variables depicting arguments for commands are always structured like this (I commented on this):
{
struct Map_info In, Out; /* info on map format */
struct line_pnts *Points, *BPoints; /* Structures to hold line *Points (old map)
*Bpoints (buffer map) see Dig_Structs.h */
struct line_cats *Cats; /* structure to hold Cats see Dig_Structs.h*/
char *mapset; /* pointer to mapset */
struct GModule *module; /* Simple description of module see GiS.h */
struct Option *in_opt, *out_opt, *type_opt, *buffer_opt, *tolerance_opt, *debug_opt;
/* pointers to program arguments */
double buffer, tolerance; /* declaration to hold answer you give to
option arguments */
int type, debug; /* declaration to hold answer you give to
option arguments */
int ret, nareas, area, nlines, line; /* declaration to hold answers to function
calls */
char *Areas, *Lines; /* declaration to hold answer you give to
option arguments */
There seems two are always be defined, being struct GModule * and struct Option *. GModule is for outputting a description of the module and what it does. Next function therefore always assigned to pointer
- module:
module=G_define_module() /* in GisLib ParseR.c */
module->description = "description"
Option refers usually to input and/or output map and other arguments. To have a better understanding see GiS.h as it also provides a way to define flags. So here also assigning to pointers, although some standard options have been constructed (vector map in/out etc).
in_opt = G_define_standard_option(G_OPT_V_INPUT);
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
snip see source code
debug_opt = G_define_option(); /* assign new (not standard) option see GiS.h
and GisLib ParseR.c */
debug_opt->key = "debug";
debug_opt->type = TYPE_STRING;
debug_opt->required = NO;
debug_opt->options = "buffer,clean";
debug_opt->description = "Stop the process in certain stage.";
This last option shows the main items of a Option structure. Key is argument given on the command line, type indicates what kind of answer is to be expected (string, int etc), required tells is an answer must be given in order to process command, options is to compare the answer given where multiple answers are possible, description is self explanatory.
Also Map_info could be considered as standard, see dig_structs.h as it gives info on all supported vector formats.
After declaring all variables to arguments the next function is:
G_gisinit(argv[0]) /* I'm not sure, seems some initialization with
the argument being program name see GisLib
GisInit.c */
And then most important the call to parser:
if (G_parser (argc, argv)) /* parsing everything and also provides a call
to go interactively in shell or Gui IF no
arguments are given and if its possible for
given module in GiSlib ParseR.c, test if
argc is same as argv*/
exit(EXIT_FAILURE); /* if not error code */
Input and output file names should be checked to make sure that the input file exists, the output file is not the input file, and the output filename is legal. The third argument can be one of the following:
- GR_FATAL_EXIT - Prints error and exits program
- GR_FATAL_PRINT - Prints warning and returns error code
- GR_FATAL_RETURN - Returns error code
Vect_check_input_output_name(input_name, output_name, GR_FATAL_EXIT);
Now passing answers to to options before doing the work:
type = Vect_option_to_types ( type_opt ); /* get type see VliB TypE.c here type refers
to point, line, area etc */
buffer = abs ( atof( buffer_opt->answer ) ); /* get buffer answer and convert to float */
tolerance = atof( tolerance_opt->answer ); /* get tolerance answer and convert to float */
debug = DEBUG_NONE;
if ( debug_opt->answer ) {
if ( debug_opt->answer[0] == 'b' )
debug = DEBUG_BUFFER;
else if ( debug_opt->answer[0] == 'c' )
debug = DEBUG_CLEAN;
} /* decide if option debug and what kind of
debug */
Points = Vect_new_line_struct (); /* assign to pointer (link to VliB line.c) */
BPoints = Vect_new_line_struct (); /* assign to pointer (link to VliB line.c) */
Cats = Vect_new_cats_struct (); /* assign to pointer (link to VliB cats.c) */
This all being done open the vector map
/* open input vector */
if ((mapset = G_find_vector2 (in_opt->answer, "")) == NULL) {
G_fatal_error ( "Could not find input map <%s>\n", in_opt->answer);
} /* see GisLib Find_Vect.c */
Vect_set_open_level (2); /* determine level 1 or 2 see VliB OpeN.c */
Vect_open_old (&In, in_opt->answer, mapset); /* Open map, take a good look at &In dealing
with reference? see VliB OpeN.c */
Vect_set_fatal_error (GV_FATAL_PRINT); /* error code see VliB ErroR.c */
if (0 > Vect_open_new (&Out, out_opt->answer, 0) ) {
Vect_close (&In);
exit (EXIT_FAILURE);
} /* open new map for buffer if error close all
see VliB OpeN.c and VliB ClosE.c */
Some housekeeping: copying of header and history data
Vect_copy_head_data (&In, &Out); /* see VliB Init_Head.c */
Vect_hist_copy (&In, &Out); /* copy see VliB HisT.c */
Vect_hist_command ( &Out ); /* actual write see VliB HisT.c */
I'm wondering, haven't figured it out yet, if some of these aspects can be hidden better. Take the history and header part, can we not provide for one command which does it all? Also defining options, flags etc seems unnecessary difficult, or am I thinking to OO (object orientated)? Maybe a good template with strict names for all elementary components should be constructed, look at difference for type, it refers to
- type of variable, eg string int etc
- or type of feature, eg point, line, string.
Could be changed to typeVar and typeFea? Comments please.
Now finally we get to the actual work, making the buffer.
/* Create buffers' boundaries */
/* Lines */
if ( (type & GV_POINTS) || (type & GV_LINES) ) {
int nlines, line, ltype; /* declaring number of, line itself and
line type. somehow for points or lines? */
nlines = Vect_get_num_lines ( &In ); /* see VliB Level_Two.c */
fprintf ( stderr, "Lines buffers ... ");
for ( line = 1; line <= nlines; line++ ) {
G_debug ( 3, "line = %d", line ); /* precious system debug see GisLib DebuG.c */
G_percent ( line, nlines, 2 ); /* give percentage see GisLib PercenT.c */
ltype = Vect_read_line (&In, Points, NULL, line); /* reading line see VliB ReaD.c
if ( !(ltype & type ) ) continue; /* comparing type instructed (option argument)
and actual line type */
Vect_line_buffer ( Points, buffer, tolerance, BPoints ); /* perform the math see VliB
BuffeR.c */
Vect_write_line ( &Out, GV_BOUNDARY, BPoints, Cats ); /* write the line to output
map see Vlib WritE.c */
}
fprintf ( stderr, "\n"); /* print newline */
}
Continuing with areas, same procedure in basis, just other variables and functions
if ( type & GV_AREA ) {
int i, nareas, area, centroid, nisles, isle; /* declare counter i, number of areas, area
itself, centroid, number of isle, isle
itself. this is strictly areas */
nareas = Vect_get_num_areas ( &In ); /* see VliB Level_Two.c */
fprintf ( stderr, "Areas buffers ... ");
for ( area = 1; area <= nareas; area++ ) {
G_percent ( area, nareas, 2 ); /* give percentage see GisLib PercenT.c */
centroid = Vect_get_area_centroid ( &In, area ); */ get centroid for given area see
VliB Level_Two.c */
if ( centroid == 0 ) continue; /* first found? */
/* outer ring */
Vect_get_area_points ( &In, area, Points ); /* see Vlib area.c */
Vect_line_buffer ( Points, buffer, tolerance, BPoints ); /* perform the math see VliB
BuffeR.c */
Vect_write_line ( &Out, GV_BOUNDARY, BPoints, Cats ); /* write the line to output
map see VliB WritE.c */
/* islands */
nisles = Vect_get_area_num_isles (&In, area); /* see VliB AreA.c */
for ( i = 0; i < nisles; i++ ) {
isle = Vect_get_area_isle (&In, area, i); /* see VliB AreA.c */
Vect_get_isle_points ( &In, isle, Points ); /* get the points see VliB AreA.c */
Vect_line_buffer ( Points, buffer, tolerance, BPoints ); */ perform the math see
VliB BuffeR.c */
Vect_write_line ( &Out, GV_BOUNDARY, BPoints, Cats ); /* write the line to
output map see Vlib
WritE.c */
}
}
fprintf ( stderr, "\n");
}
if ( debug == DEBUG_BUFFER ) {
stop ( &In, &Out );
exit (EXIT_SUCCESS);
} /* if defined debug stop */
Now some after work must be done before writing and building the map. Radim has put some thoughts here, so see if they make sense. It's mainly to do with building a new topology, a mystery to most. In general in grass vector handling there are two levels of interaction. One that deals with simple geometry (level 1) and one with geometry and topology (level 2). This not being documented well and discussed lately on grass list (Jan,Feb 2004), is still confusing especially when dealing with code. Clearly new common boundaries have to be defined and duplicates deleted etc, and some issues with snapping. The point being that all buffer areas are now new areas, so the new map has old areas and new buffer areas.
/* Create areas */
/* Break lines */
fprintf (stderr, "Building parts of topology ...\n" );
Vect_build_partial ( &Out, GV_BUILD_BASE, stderr ); /* Partial Building the map for main area
not isles?, I'm still confused, see
VliB BuilD.c */
/* Warning: snapping must be done, otherwise collinear boundaries are not broken and
* topology cannot be built (the same angle). But snapping distance must be
very, very
* small, otherwise counterclockwise boundaries can appear in areas outside
the buffer.
* I have done some tests on real data (projected) and threshold 1e-8 was not enough,
* Snapping threshold 1e-7 seems to work. Don't increase until we find example
* where it is not sufficient. */
/* TODO: look at snapping threshold better, calculate some theoretical value to avoid
* the same angles of lines at nodes, don't forget about LongLat data, probably
* calculate different threshold for each map, depending on map's bounding box */
fprintf ( stderr, "Snapping boundaries ...\n" );
Vect_snap_lines ( &Out, GV_BOUNDARY, 1e-7, NULL, stderr ); /* see VliB SnaP.c */
fprintf ( stderr, "Breaking boundaries ...\n" );
Vect_break_lines ( &Out, GV_BOUNDARY, NULL, stderr ); /* see VliB Break_Lines.c */
fprintf ( stderr, "Removing duplicates ...\n" );
Vect_remove_duplicates ( &Out, GV_BOUNDARY, NULL, stderr ); /* defining new common boundaries
and removing duplicates?, see
VliB Remove_Duplicates.c */
More thoughts and dealing with Isles
/* Dangles and bridges don't seem to be necessary if snapping is small enough. */
fprintf ( stderr, "Removing dangles ...\n" );
Vect_remove_dangles ( &Out, GV_BOUNDARY, -1, NULL, stderr ); /* see VliB DangleS.c */
fprintf ( stderr, "Removing bridges ...\n" );
Vect_remove_bridges ( &Out, NULL, stderr ); /* see VliB BridgeS.c */
fprintf ( stderr, "Attaching islands ...\n" );
Vect_build_partial ( &Out, GV_BUILD_ATTACH_ISLES, stderr ); /* again build partial now for
Isle see VliB BuilD.c */
/* Calculate new centroids for all areas */
nareas = Vect_get_num_areas ( &Out ); /* see VliB AreA.c */
Areas = (char *) G_calloc ( nareas+1, sizeof(char) ); /* dealing with memory for first time
why necessary here? see GisLib
alloc.c */
for ( area = 1; area <= nareas; area++ ) {
G_debug ( 3, "area = %d", area ); /* see GisLib DebuG.c */
ret = area_in_buffer ( &In, &Out, area, type, buffer, tolerance ); /* Function in main.c
of v.buffer, is not
yet explained but
checking if area in
output map is in
buffer */
if ( ret ) {
G_debug ( 3, " -> in buffer"); /* see GisLib DebuG.c */
Areas[area] = 1;
} /* if in buffer (valid) assign 1 */
/* Write out centroid (before check if correct, so that it is visible for debug option
argument) */
if ( debug == DEBUG_CLEAN ) {
double x, y;
ret = Vect_get_point_in_area ( &Out, area, &x, &y ); /* new centroid see VliB PolY.c */
if ( ret < 0 ) {
G_warning ("Cannot calculate area centroid" );
continue;
} /* Checking again */
Vect_reset_cats ( Cats ); /* see VliB CatS.c */
if ( Areas[area] )
Vect_cat_set (Cats, 1, 1 ); /* see VliB CatS.c */
Vect_reset_line ( Points ); /* some reset? see VliB LinE.c */
Vect_append_point ( Points, x, y, 0 ); /* append points for writing to
new map see VliB LinE.c */
Vect_write_line ( &Out, GV_CENTROID, Points, Cats ); /* write centroids to file see
VliB WritE.c */
}
}
if ( debug == DEBUG_CLEAN ) {
stop ( &In, &Out );
exit (EXIT_SUCCESS);
} /* again something for Option debug argument */
/* Make a list of boundaries to be deleted (both sides inside) */
nlines = Vect_get_num_lines ( &Out ); /* see Vlib Level_Two.c */
G_debug ( 3, "nlines = %d", nlines ); /* see GisLib DebuG.c */
Lines = (char *) G_calloc ( nlines+1, sizeof(char) ); /* again allocating memory is beyond me
now */
for ( line = 1; line <= nlines; line++ ) {
int j, side[2], areas[2];
G_debug ( 3, "line = %d", line ); /* see GisLib DebuG.c */
if ( !Vect_line_alive ( &Out, line ) ) continue; /* check for validness, see VliB ReaD.c */
Vect_get_line_areas ( &Out, line, &side[0], &side[1] ); /* checking for number of areas or
isle on right or left side of
lines see VliB Level_Two */
for ( j = 0; j < 2; j++ ) {
if ( side[j] == 0 ) { /* area/isle not build */
areas[j] = 0;
} else if ( side[j] > 0 ) { /* area */
areas[j] = side[j];
} else { /* < 0 -> island */
areas[j] = Vect_get_isle_area ( &Out, abs ( side[j] ) ); */ not null, not area,
must be isle see VliB
AreA.c */
}
}
G_debug ( 3, " areas = %d , %d -> Areas = %d, %d", areas[0], areas[1], Areas[areas[0]], Areas[areas[1]] );
/* see GisLib DebuG.c */
if ( Areas[areas[0]] && Areas[areas[1]] )
Lines[line] = 1; /* assign lines for deletion */
}
G_free( Areas ); /* release memory see GisLib AlloC.c */
/* Delete boundaries */
for ( line = 1; line <= nlines; line++ ) {
if ( Lines[line] ) {
G_debug ( 3, " delete line %d", line ); /* see GisLib DebuG.c */
Vect_delete_line ( &Out, line ); /* see VliB WritE.c */
}
}
G_free( Lines ); /* release memory see GisLib AlloC.c */
/* Create new centroids */
Vect_reset_cats ( Cats ); /* some resetting see VliB CatS.c */
Vect_cat_set (Cats, 1, 1 ); /* writing cats see VliB CatS.c */
nareas = Vect_get_num_areas ( &Out ); /* See VliB AreA.c */
for ( area = 1; area <= nareas; area++ ) {
double x, y;
G_debug ( 3, "area = %d", area ); /* see GisLib DebuG.c */
if ( !Vect_area_alive ( &Out, area ) ) continue; /* see VliB ReaD.c */
ret = area_in_buffer ( &In, &Out, area, type, buffer, tolerance ); /* check areas in
output map is in
buffer see main.c in
v.buffer */
if ( ret ) {
ret = Vect_get_point_in_area ( &Out, area, &x, &y ); /* get centroid in area not in
isle see Vlib poly.c */
if ( ret < 0 ) {
G_warning ("Cannot calculate area centroid." ); /* see LibGis ErroR.c
continue;
}
Vect_reset_line ( Points ); /* see Vlib LinE.c */
Vect_append_point ( Points, x, y, 0 ); /* see Vlib LinE.c */
Vect_write_line ( &Out, GV_CENTROID, Points, Cats ); /* write the centroids see Vlib
WritE.c
}
fprintf ( stderr, "Attaching centroids ...\n" );
Vect_build_partial ( &Out, GV_BUILD_CENTROIDS, stderr ); /* again build partial topology but
now for centroids, still
confusing, is everything build
now? */
stop ( &In, &Out ); /* defined in main.c of v.buffer */
exit(EXIT_SUCCESS); /* exit program */
}
The workings of the main function should at least be somewhat clear now. Most functions are in the gis or vector library, however there remain the new functions of the module itself being area_in_buffer() to check whether area in output map is in buffer. Basically compare area in input map with output map. This function will also call another function input_distance() being defined in main.c of v.buffer and has something to do with defining the distance from a point of an area in the output map to the nearest feature in the input map. I'm myself not quite sure I understand what is happening, if you do add you're insights. If I understand it's defining the distance to the nearest feature in the buffer and if none are found outside the buffer. If I try to rationalize this, it makes a little more sense. Let's take a look at this:
+ - + + - + + - +
+ - + A + - + B + - +
+ - + + - + + - +
+ = buffer
- = borderline
area = centroid
A centroid may fall into buffer but the area as a whole doesn't have to be. So distance to centroid can be smaller than buffer or bigger. That's I believe the main rationale for these functions. Again we'll walk through the code starting with area_in_buffer()
/*
* Test if area in Out is in buffer.
* Return: 1 in buffer
* 0 outside buffer
*/
int area_in_buffer ( struct Map_info *In, struct Map_info *Out, int area, int type, double
buffer, double tolerance ) /* taking arguments which should by now mean
something */
{
double dist; /* think declaration for something returned by
function input_distance */
int ret, i; /* seen before, declaration for return value
and counter */
static struct ilist *List = NULL; /* declaration of list of integers see
Dig_Struct.h */
static struct line_pnts *Points = NULL; /* setting point structure to NULL see
Dig_Struct.h */
double x, y; /* declaring as double x and y coordinates? */
G_debug ( 3, " area_in_buffer(): area = %d", area ); /* see GisLib DebuG.c */
if ( List == NULL )
List = Vect_new_list (); /* make new list see VliB LisT.c */
if ( Points == NULL )
Points = Vect_new_line_struct (); /* new line structure see Vlib LinE.c */
/* Warning: because of tolerance, centroid in area calculated by Vect_get_point_in_area()
* may fall into buffer (arc interpolated by polygon) even if area is outside!!! */
/* Because of finite number of decimal places in double representation, RET is used. */
/* Test1: Centroid (the only reliable, I think).
* There are 3 possible cases for the distance of centroid to nearest input feature:
* 1) < (buffer - tolerance) => area inside the buffer
* 2) > buffer => area outside the buffer
* 3) > (buffer - tolerance) and < buffer (that means within the tolerance) => uncertain */
ret = Vect_get_point_in_area ( Out, area, &x, &y ); /* get coordinates for a point in an area
in output map see VliB PolY.c */
if ( ret == 0 ) { /* centroid OK */
/* test the distance to nearest input feature */
dist = input_distance ( In, type, buffer, x, y); /* this function is in main.c of
v.buffer explained later */
G_debug ( 3, " centroid %f, %f dist = %f", x, y, dist ); /* see GisLib DebuG.c */
if ( dist < (buffer - tolerance - RET) ) {
G_debug ( 3, " centroid in buffer -> area in buffer"); /* see GisLib DebuG.c */
return 1;
} else if ( dist > (buffer + RET) ) {
G_debug ( 3, "centroid outside buffer -> area outside buffer"); /* see GisLib DebuG.c */
return 0;
} /* testing for case 1 as described above by
Radim */
/* Test2: counterclockwise (ccw) boundary
* Boundaries around features are written in ccw order, that means area inside buffer is on the
* left side. Area boundaries are listed in cw order (ccw boundaries as negative).
* So if any boundary in area list is negative, i.e. in cww order, i.e. buffer inside area,
* whole area _must_ be in buffer. But, also areas without such boundary may be in buffer,
* see below. */
/* TODO: The problem!!!: unfortunately it is not true. After snap, break and remove duplicate
* it may happen that ccw line appears in the area outside the buffer!!! */
Vect_get_area_boundaries ( Out, area, List ); /* get boundaries of areas of in output map and
put in list plus assign an integer id
number for looping through them see VliB
AreA.c */
for ( i = 0; i < List->n_values; i++ ) {
if ( List->value[i] < 0 ) {
G_debug ( 3, " negative boundary -> area in buffer"); /* see GisLib DebuG.c */
return 1;
}
} /* looping through the list to find negative
areas */
/* Test3: Input feature within area.
* I believe, that if no negative boundary was found and area is in buffer,
* the distance of the nearest input feature for at least one area vertex must be < (buffer-
tolerance)
* This test is left as last because it is slow (check each vertex). */
Vect_get_area_points ( Out, area, Points ); /* get the points=vertex of line see VliB
AreA.c */
for ( i = 0; i < Points->n_points - 1; i++ ) {
dist = input_distance ( In, type, buffer, Points->x[i], Points->y[i]);
G_debug ( 3, " vertex dist = %f", dist ); /* see GisLib DebuG.c */
if ( dist < (buffer - tolerance - RET) ) {
G_debug ( 3, " vertex in buffer -> area in buffer"); /* see GisLib DebuG.c */
return 1;
}
} /* looping though the list and doing the math */
/* Test4?: It may be that that Test3 does not cover all possible cases. If so, somebody must
* find such example and send it to me */
G_debug ( 3, " -> area outside buffer"); /* see GisLib DebuG.c */
return 0;
}
void stop ( struct Map_info *In, struct Map_info *Out ) /* function stop called in main() */
{
Vect_close (In); /* see VliB ClosE.c */
fprintf ( stderr, "Rebuilding topology ...\n" );
Vect_build_partial ( Out, GV_BUILD_NONE, NULL ); /* see VliB BuilD.c */
Vect_build (Out, stderr); /* see VliB BuilD.c */
Vect_close (Out); /* see VliB ClosE.c */
} /* closing the input map and building and closing the output map */
Ok now for the last function which support the last one
described. Though I do not fully understand, it has to with topology
and the distance of the a point to a other feature in the map.
double input_distance ( struct Map_info *In, int type, double buffer, double x, double y )
{
int i; /* for looping */
static struct ilist *List = NULL; /* setting the list to null */
static struct line_pnts *Points = NULL; /* setting the list of line points to null */
double current_distance; /* declaration of */
BOUND_BOX box; /* see Dig_Structs.h assigning name box */
if ( List == NULL )
List = Vect_new_list (); /* See Vlib LisT.c */
else
Vect_reset_list ( List );
if ( Points == NULL )
Points = Vect_new_line_struct (); /* See Vlib LinE.c */
/* Inside area ? */
if ( type & GV_AREA) {
int area, centroid;
/* inside */
area = Vect_find_area ( In, x, y ); /* find a point in an area see Vlib find.c */
centroid = 0;
if ( area )
centroid = Vect_get_area_centroid ( In, area ); /* get the centroid of that area see VLiB
AreA.c */
G_debug ( 3, "area = %d centroid = %d", area, centroid ); /* see GisLib DebuG.c */
if ( centroid )
return 0.0;
} /* find the area and the centroid, if none return 0.0 */
/* outside area, within buffer? */
/* The centroid is in buffer if at least one point/line/boundary is in buffer distance */
box.E = x + buffer; box.W = x - buffer;
box.N = y + buffer; box.S = y - buffer;
box.T = PORT_DOUBLE_MAX; box.B = -PORT_DOUBLE_MAX; /* defining the bounding box see
Dig_Structs.h */
Vect_select_lines_by_box ( In, &box, GV_POINTS | GV_LINES, List ); /* see Vlib SelecT.c */
G_debug ( 3, " %d lines selected by box", List->n_values ); /* see GisLib DebuG.c */
current_distance = 2*buffer; /* defining the the distance as 2* buffer */
for ( i = 0; i < List->n_values; i++) {
int line, ltype; /* declare line and line type */
double dist; /* declare distance */
line = List->value[i]; /* make line addressable in list by a number */
G_debug ( 3, " line[%d] = %d", i, line ); /* see GisLib DebuG.c */
ltype = Vect_read_line (In, Points, NULL, line); /* get line types in input map see VliB
ReaD.c */
Vect_line_distance ( Points, x, y, 0, 0, NULL, NULL, NULL, &dist, NULL, NULL);
/* distance of a line? see VliB line.c */
G_debug ( 3, " dist = %f", dist ); /* see GisLib DebuG.c */
if ( dist > buffer ) continue;
/* lines */
if ( type & ltype ) {
if ( dist < current_distance )
current_distance = dist;
} /* test for initial distance if greater than buffer or smaller */
/* Areas */
if ( (type & GV_AREA) && ltype == GV_BOUNDARY ) {
int j, side[2], centr[2], area_in;
Vect_get_line_areas ( In, line, &side[0], &side[1] ); /* getting the lines for an area see
VliB Level_Two.c */
for ( j = 0; j < 2; j++ ) {
centr[j] = 0;
if ( side[j] > 0 )
area_in = side[j];
else /* island */
area_in = Vect_get_isle_area ( In, abs ( side[j] ) ); /* testing for isles see VliB
AreA.c */
if ( area_in > 0 )
centr[j] = Vect_get_area_centroid ( In, area_in ); /* get centroid see VliB AreA.c */
}
if ( centr[0] || centr[1] ) {
if ( dist < current_distance )
current_distance = dist;
}
}
}
return current_distance;
}
OK we should be able now to at least have a little understanding of what is happening now. Only thing we did not take a look at are the defines:
#include <stdlib.h> /* standard lib unixsystem */
#include <string.h> /* standard string lib unixsystem */
#include "gis.h" /* include gis.h standard Grass system for
all the structs we've seen */
#include "Vect.h" /*vect.h standard Grass system for all the
structs we've seen */
#define DEBUG_NONE 0 /* for the option debug */
#define DEBUG_BUFFER 1 /* for the option debug */
#define DEBUG_CLEAN 2 /* for the option debug */
/* TODO: look at RET value and use, is it OK? */
#define RET 0.000000001 /* Representation error tolerance */
/* for the math if the area centroid is in or outside buffer */
Well that's it for now. I'll hope you get the grasp and help out a little because it's by no means perfect now. Also a better system of referral must be thought of. I'm thinking of numbering everything like VliB 1.1 (1.= e.g AreA.c & .1 = first function declared) etc.
Sites Modules
Sites are obsolete in GRASS6. Use the Vector Library.
Image Data Modules
Everything related to Image Data Modules
Region and Mask Modules
Everything related to Region and Mask Modules
Environment Variables Modules
Everything related to Environment Variables Modules
Paint Modules
These are mostly retired in GRASS 6. Functionality is covered by either ps.map or display commands followed by screen output or d.out.png. For labeling vector fetaures see v.label and d.labels. See also d.graph.
Display Symbols
For user contributed symbols see the IconSymbols wiki page.
Symbol drawing commands are stored in a text file (one file for each symbol)
which are located in $GISBASE/etc/symbol/group/name or mapset/symbol/group/name
S_read() searches first in mapsets and then in $GISBASE (symbol names in
$GISBASE may be overidden by symbols of the same group/name in mapset).
File format example:
VERSION 1.0
BOX -10 -10 10 10
STRING
LINE
0 -10
0 10
END
END
A symbol file is composed of objects, each object starts with a keyword and
ends with a new line or 'END'. Some objects may contain other objects.
- Drawing is composed of POLYGONS and STRINGS.
- POLYGONS are composed of RINGS. The first ring is outer (clockwise for now),
following rings are inner (counter clockwise currently).
- STRINGS and RINGS are composed of LINES and ARCS.
POLYGONS and STRINGS are drawn using the default color and fill color
(specified in module options such as 'd.vect color='). This color may be
overwritten by using the COLOR and FCOLOR keywords.
---- DESCRIPTION OF OBJECTS ----
In the following descriptions these variables are used:
N - integer number
R G B - RGB color as 3 numbers ranging 0-255
X,Y - coordinates (double precision floating point)
Keywords:
VERSION N.N
Version number
BOX X Y X Y
Bounding box. Lower left and upper right corner. Larger side of box is taken
as 1 for symbol scaling. That means that if size of box is 10 x 5, and scale
in application is 2, all geometry will be multiplied by: 2 / 10.
Origin of symbol is 0,0.
COLOR NONE | R G B
Color used for string or outline of polygon. If NONE is used, then the polygon
outline is not drawn.
FCOLOR NONE | R G B
Fill color used for polygon. If NONE is used, then the polygon is not filled.
LINE
X Y
X Y
[X Y]
END
Line constructed from coordinate pairs.
ARC X Y R A1 A2 [C]
Arc specified by center (X Y), radius (R), start/end angle in degrees (A1/A2)
and optionaly direction (C for clockwise, default counter clockwise)
Note degrees are polar theta not navigational, i.e. 0 is east and
values increase counter-clockwise from there. Reducing the angle
from a full circle does not make it draw like a pacman but rather
with a straight chord-line from the start/end points on the circle
as defined by the two angles.
STRING
[COLOR R G B]
LINE | ARC
[LINE | ARC]
END
String composed of any number of lines and/or arcs.
RING
LINE | ARC
[LINE | ARC]
END
Ring composed of any number of lines and/or arcs. A ring is automatically
closed with a line from last point to first point.
POLYGON
[COLOR R G B]
[FCOLOR R G B]
RING
[RING]
END
A polygon is composed of any number of rings. The first ring is outer
(clockwise) and others are inner (counter-clockwise).
------------------------------------
Note: restriction for clockwise/counter clockwise rings should be removed
later and replaced by reordering in symbol library. Order is required by
ps.map because PS uses even-odd rule for polygon filling.
Internal structure:
A symbol read into memory is stored in structure SYMBOL containing an array
of parts (SYMBPART). Each part may be of type S_STRING or S_POLYGON, and
contain an array of chains (SYMBCHAIN) where both STRING or POLYGONS are
stored. One chain contains an array of elements (SYMBEL).
Compiling and Installing
Everything related to Compiling and Installing. Makefiles, configure files etc.
How to structure
Numerical math interface
GUI programming
Digitizer/Mouse/Trackball