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;
}