[GRASS5] Nested vector islands, and incorrect left/right poly id's for island edges...

[David D Gray]

Eric G. Miller wrote:

> It appears that GRASS doesn't write correct information for
> left/right adjaceny for edges of areas nested more than
> one polygon deep.  That is: If polygon 1, contains polygon
> 2, and polygon 2 contains polygon 3, the lines comprising
> polygon 3 will have polygon 1 written as the adjacent area
> for either left or right, when I would think it should be
> polygon 2.
> 


Yes I would have thought so.


> I was trying to work up a way to do the following:
> 
>   Given an area with N islands/holes.  Create a set of
>   interior boundaries such that islands/holes that
>   share one or more edges are combined into a single
>   boundary (with interior edges being dropped).
> 


I don't quite understand how islands can share edges. I thought in the 
earlier posting you meant that an edge of an island can be the same as 
that of an (outer) ring. If two 'islands' share a boundary, then they 
are actually the constituent cells of an island, not the island itself. 
If the build routine is creating islands within a polygon that abut, 
then that is a bug.

I have recently been developing a terminology that helps to create a 
more formal link between the intuitive notions and the topological and 
geometric concepts in common circulation.

Four concepts:

Holes <-------------------------> Hulls

   |                                  |

  |                                  |
  |                                  |
  |                                  |
  |                                  |

Islands <-------------------------> Cells



Holes and hulls are defined in relation to an individual polygon. In 
practical terms, a polygon is the interior of a fully connected open 
subspace of the whole, union'd with its closure. Without going into that 
too deeply, it can be thought of as the boundary, but purely in terms of 
a limiting bound to the interior, with no conception of edges, nodes 
etc. Exterior boundaries (only one) consitutes the 'hull' and interior 
boundaries the 'holes'. You can tell if its a hole or a hull by tracknig 
from an arbitrary interior point to a boundary point on the subject 
boundary. Then track right, and follow the boundary round to its 
starting point  (we've dropped to geometrical operations here). If the 
boundary is traversed clockwise it is a hull; if anti-clockwise it is a 
hole.

Islands and cells are defined in relation to the boundary system of a 
coverage, which really is a network with some constraints. You start at 
a node or track from an arbitrary point to a node, and continue 
tracking, always following the right-most link at each node, until you 
get back to the starting point. If you circulate clockwise, it is a 
cell, and if anti-clockwise, it's an island. An 'island' is the boundary 
  set of a subspace whose interior is an 'aggregate' of cells, union'd 
with that interior. A 'cell' is a subspace whose interior union'd with 
the cell boundary co-incides with a polygon. If an island has only a 
single cell it is a 'simple island' (though it can have islands of its 
own), and is also a cell in its own right.

The first set of definitions is similar in structure to the kind of 
polygon coverage that is found in shapefiles for examples, while the 
latter is like that found in GRASS, Arc/Info etc.

This is quite detailed, as it ties in to formal topological theory, but 
it's useful for us as it highlights some important areas - how the 
different kinds of coverage that are used relate to different aspects of 
the topology of the coverage, polygons relate to interior sets, 
'topological' coverages to boundary sets. There's also the practical 
trick of tracking right or left to distinguish inner and outer components.


> I tried to write a little program that iterated through
> areas, and their islands, reporting the lines for all
> and the reported line left/right indices.  That's when
> I discovered, I was getting islands where none of the
> edges reported being adjacent to their immediate parent.
> 
> I guess, next I'll try using nodes instead of edges.  If
> anyone familiar with the guts of the vector libraries
> (Radim, David ?) wants to point me in the right direction,
> I'd appreciate it...
> 


I haven't tried this yet - I'll look now and see if I can confirm this 
bug. We have to sort that out priority A1 if there is, and then I think 
that should solve the problem.

David


> Attached:  area_report.c & Gmakefile
> 
> 
> 
> ------------------------------------------------------------------------
> 
> #include <stdio.h>
> #include <stdlib.h>
> #include "gis.h"
> #include "Vect.h"
> 
> int main (int argc, char **argv)
> {
> 	struct Option *opt_map;
> 	struct Map_info map;
> 	char *mapset;
> 	int i, j, k, att;
> 	int isle_area;
> 	P_AREA *area;
> 	P_LINE *line;
> 	P_ISLE *isle;
> 	
> 	opt_map = G_define_option();
> 	opt_map->key = "map";
> 	opt_map->required = YES;
> 	opt_map->gisprompt = "dig,old,vector";
> 
> 	G_gisinit(argv[0]);
> 
> 	if (G_parser (argc, argv))
> 		exit (EXIT_FAILURE);
> 
> 	if (NULL == (mapset = G_find_file2 ("dig", opt_map->answer, "")))
> 		G_fatal_error ("Vector map [%s] not found!", opt_map->answer);
> 
> 	if (V2_open_old (&map, opt_map->answer, mapset))
> 		G_fatal_error ("Couldn't open vector map at level 2!");
> 
> 	for (i = 1; i <= map.n_areas; i++)
> 	{
> 		fprintf (stdout, "Area %d {\n", i);
> 		att = V2_area_att (&map, i);
> 		fprintf (stdout, "  dig_att = %d\n", att);
> 		
> 		V2_get_area (&map, i, &area);
> 
> 		fprintf (stdout, "  lines {\n");
> 		for (j = 0; j < area->n_lines; j++) 
> 		{
> 			line = &(map.Line[abs(area->lines[j])]);
> 			fprintf (stdout, "    line %d, left %d, right %d\n",
> 					area->lines[j],
> 					line->left, line->right);
> 		}
> 		fprintf (stdout, "  }\n");
> 
> 		for (j = 0; j < area->n_isles; j++)
> 		{
> 			fprintf (stdout, "  Island %d {\n", area->isles[j]);
> 			isle = &(map.Isle[area->isles[j]]);
> 			fprintf (stdout, "    Lines {\n");
> 			for (k = 0; k < isle->n_lines; k++)
> 			{
> 				line = &(map.Line[abs(isle->lines[k])]);
> 				fprintf (stdout, "      Line %d, LeftArea %d, RightArea %d\n",
> 						isle->lines[k],
> 						line->left,
> 						line->right);
> 			}
> 			if (i == abs(line->left)) 
> 				isle_area = abs(line->right);
> 			else if (i == abs(line->right))
> 				isle_area = abs(line->left);
> 			else 
> 				isle_area = -1;
> 			
> 			fprintf (stdout, "    } = Area %d\n  }\n", isle_area );
> 		}
> 		
> 		fprintf (stdout, "}\n");
> 	}
> 			
> 			
> 			
> 	return 	0;
> }
> 
>