[GRASS5] [5.7] remove islands and U-dangles?

GRASS GRASS developer
1

Hi,

more questions about 5.7 vectors....

I've got a vector area which contains a bunch of islands and polygons
attached to the inside edge of the outer polygon. What commands would
one use to remove and 'flatten' the overlapping areas?

eg

/-------+---+----------\
| | X | |
| | X | <- area |
| |---| |
| |
| |=| <- island |
| |
\----------------------/

In the end I want a convex hull of the outer polygon as a single entity,
without losing resolution from a points->v.hull step.

I've tried a few things with v.overlay, but they just produced empty
maps (No tables or databases, probably bad topology in there.. d.vect
does draw the input as an area though, and the output coor file is huge
while other vector files are just a couple of hundred bytes).

'v.overlay op=or' might do the trick if it made 1 polygon out of two
input polygons instead of three and I could use the same map for ainput=
and binput=.

thanks for any ideas,
Hamish

ps [grass5 list]:
  I fixed the hidden dangling node bug in 5.7's v.digit.

2

On Thu, 12 Feb 2004 10:06:09 +0100
Radim Blazek <blazek@itc.it> wrote:

On Thursday 12 February 2004 06:42, you wrote:
> > > any chance I could try the buggy version of extract.c? I think
> > > '-d' will do exactly what I need, and I don't care about
> > > centroids.
> >
> > Go to v.extract/main.c and enable d_flag (commented),
> > then in v.extract/extract.c in block starting at the row 69 add
> > something which checks if area on both sides are exracted.
>
> I'll try, but I'm not too sure how to do "check if area on both
> sides are exracted" yet..

You have to check if both catr and catl are in num_array (number of
items num_index) so you need a loop as num_array is not sorted. It is
awful v.extract should be rewritten a bit I think, maybe
Vect_str_to_cat_list(), Vect_array_to_cat_list(),
Vect_cat_in_cat_list().

see also http://article.gmane.org/gmane.comp.gis.grass.user/3870

The following patch makes the -d (dissolve) function of v.extract work.
(this is against today's CVS)

It works, but probably needs some cleanup before putting in CVS..
(i.e. I've taken it as far as I can, someone else please do it)

It's faster than I thought it would be.. (for 15,000 cats)

Usage: v.extract -d in=vect1 out=vect2 type=area,boundary list=1-50000

Also, 'v.extract type=area list=1,2,3' keeps areas with missing
centroids; how do you automatically add missing centroids so you can
then get rid of the zombie areas with another 'v.extract list=1,2,3'?

These are islands within a bigger area, in the following I want to get
rid of "A2". The outer boundary is made up of multiple lines. A1 has a
centroid.

+------------+
| |
| +--+ |
| |A2| |
| +--+ + |
| A1 |
+------------+

Hamish

(attachments)

v.extract_dissolve.diff (2.59 KB)

3

On Tue, 11 May 2004 20:59:25 +1200
Hamish <hamish_nospam@yahoo.com> wrote:

On Thu, 12 Feb 2004 10:06:09 +0100
Radim Blazek <blazek@itc.it> wrote:

> On Thursday 12 February 2004 06:42, you wrote:
> > > > any chance I could try the buggy version of extract.c? I think
> > > > '-d' will do exactly what I need, and I don't care about
> > > > centroids.
> > >
> > > Go to v.extract/main.c and enable d_flag (commented),
> > > then in v.extract/extract.c in block starting at the row 69 add
> > > something which checks if area on both sides are exracted.
> >
> > I'll try, but I'm not too sure how to do "check if area on both
> > sides are exracted" yet..
>
> You have to check if both catr and catl are in num_array (number of
> items num_index) so you need a loop as num_array is not sorted. It is
> awful v.extract should be rewritten a bit I think, maybe
> Vect_str_to_cat_list(), Vect_array_to_cat_list(),
> Vect_cat_in_cat_list().

see also http://article.gmane.org/gmane.comp.gis.grass.user/3870

The following patch makes the -d (dissolve) function of v.extract work.
(this is against today's CVS)

It works, but probably needs some cleanup before putting in CVS..
(i.e. I've taken it as far as I can, someone else please do it)

It's faster than I thought it would be.. (for 15,000 cats)

Usage: v.extract -d in=vect1 out=vect2 type=area,boundary list=1-50000

Also, 'v.extract type=area list=1,2,3' keeps areas with missing
centroids; how do you automatically add missing centroids so you can
then get rid of the zombie areas with another 'v.extract list=1,2,3'?

These are islands within a bigger area, in the following I want to get
rid of "A2". The outer boundary is made up of multiple lines. A1 has a
centroid.

+------------+
| |
| +--+ |
| |A2| |
| +--+ + |
| A1 |
+------------+

Thank you very much for working on this problem just a few days after I sent my 'Wish' to the bugtracker!
I have tested the module patch with part of the spearfish location and made some screeshots. Maybe I can help to make it clearer.

1) http://www.gdf-hannover.de/dassau/muell/d.vect.png

Visualize spearfish soils data:
d.vect soils where='cat<20'

2) http://www.gdf-hannover.de/dassau/muell/v.extract.png

Using v.extract without dissolve option
v.extract in=soils out=soils2 type=area,centroid where='cat<20'

-> extracted areas within areas don't have an centroid anymore but are still visualized, although they should be deleted, as with d.vect. Topologically, this might be correct, because an area is surrounded by boundaries, but it should not be displayed, I think.

3) http://www.gdf-hannover.de/dassau/muell/v.extract_d.png

Using v.extract with dissolve option
v.extract -d in=soils out=soils3 type=area,centroid where='cat<20'

-> areas with cat<20 are dissolved (very nice), but now areas with dublicated centroids exist. And areas cat='21 and 43' still exist
but without centroid, as in example 2.

  Otto

4

> The following patch makes the -d (dissolve) function of v.extract
> work.(this is against today's CVS)
>
> It works, but probably needs some cleanup before putting in CVS..

I'll try and do this today..

> Also, 'v.extract type=area list=1,2,3' keeps areas with missing
> centroids; how do you automatically add missing centroids so you can
> then get rid of the zombie areas with another [xxxxxxxx]?

=> 'v.extract -d list=1-50000'

I thought 'v.clean tool=bpol' would do this, but it doesn't.
Should it? Or should/is this better in another tool?

If the number of left-over broken areas is small, then for now you can
remove them with v.digit.. :confused:

I have tested the module patch with part of the spearfish location and
made some screeshots. Maybe I can help to make it clearer.

1) http://www.gdf-hannover.de/dassau/muell/d.vect.png

Visualize spearfish soils data:
d.vect soils where='cat<20'

2) http://www.gdf-hannover.de/dassau/muell/v.extract.png

Using v.extract without dissolve option
v.extract in=soils out=soils2 type=area,centroid where='cat<20'

-> extracted areas within areas don't have an centroid anymore but are
still visualized, although they should be deleted, as with d.vect.
Topologically, this might be correct, because an area is surrounded by
boundaries, but it should not be displayed, I think.

See above. I don't know the solution to this.
Is this a display bug or a topology bug?

** note ps.map->vector fill has this problem too **

3) http://www.gdf-hannover.de/dassau/muell/v.extract_d.png

Using v.extract with dissolve option
v.extract -d in=soils out=soils3 type=area,centroid where='cat<20'

-> areas with cat<20 are dissolved (very nice), but now areas with
dublicated centroids exist. And areas cat='21 and 43' still exist but
without centroid, as in example 2.

v.clean tool=rmdac # remove duplicate area centroids

I'm not sure how it decides which one to keep, but I'd suspect the
area's cat might not be correct after this step.

Probably extract.c's dissolve needs to work on centroids somehow (right
now it's only boundaries).

Hamish

5

On Wednesday 12 May 2004 01:42, Hamish wrote:

> > The following patch makes the -d (dissolve) function of v.extract
> > work.(this is against today's CVS)
> >
> > It works, but probably needs some cleanup before putting in CVS..

I'll try and do this today..

But the problem is not to dissolve boundaries, the problem is to choose
only one centroid. OK, put it to CVS and I'll look at the centroids,
I think that I'll generate new centroids for merged areas, and query
the area in the input map. It would be also possible (faster) to simply
remove all duplicate centroids (negative LINE->left).

> > Also, 'v.extract type=area list=1,2,3' keeps areas with missing
> > centroids; how do you automatically add missing centroids so you can
> > then get rid of the zombie areas with another [xxxxxxxx]?

=> 'v.extract -d list=1-50000'

I thought 'v.clean tool=bpol' would do this, but it doesn't.
Should it? Or should/is this better in another tool?

So again, any space closed by boundary, even if it does not have any centroid
is area. It is area without centroid, without category.
We can discuss, if d.vect should (and when) display those areas, se more below.
It is not a bug! It is not a bug! It is not a bug!

v.clean tool=bpol breaks closed polygons, like those imported from Shapefile,
the result is similar as tool=break, but polygons are broken only on
common vertices. 'bpol' is much faster than 'break', but usually
(because simple feature data are not clean) 'break' must be used as well.

If the number of left-over broken areas is small, then for now you can
remove them with v.digit.. :confused:

No! The areas are enclosed by other areas which cannot be removed!

> I have tested the module patch with part of the spearfish location and
> made some screeshots. Maybe I can help to make it clearer.
>
> 1) http://www.gdf-hannover.de/dassau/muell/d.vect.png
>
> Visualize spearfish soils data:
> d.vect soils where='cat<20'
>
> 2) http://www.gdf-hannover.de/dassau/muell/v.extract.png
>
> Using v.extract without dissolve option
> v.extract in=soils out=soils2 type=area,centroid where='cat<20'
>
> -> extracted areas within areas don't have an centroid anymore but are
> still visualized, although they should be deleted, as with d.vect.
> Topologically, this might be correct, because an area is surrounded by
> boundaries, but it should not be displayed, I think.

See above. I don't know the solution to this.
Is this a display bug or a topology bug?

** note ps.map->vector fill has this problem too **

It is not a bug! Neither in display nor in topology.
We can change how d.vect works if you are sure that it
should be changed. I am not sure, it will solve your present
problem but creates another.
Say that d.vect displays only areas with category of given field
(that is what you want, I think?) then people will complain
that they cannot see areas without centroids, e.g. imported
from CAD. Here we must decide if it is better to see sometimes
more than user wants, or less than he expects. I am sure that
it is better to see more, you can think about waht you see and why.
To import data and see nothing is depressing.
You can use d.vect map=soils2 cats=1-10000 to see only areas with category.

In QGIS, one field means layer. If an area in GRASS has a category of
field 1, it is represented as feature in the layer polygon_1.
But it is different from GRASS, because user can see in the dialog,
that there is also layer polygon_0 (islands without category)
and he can add it to the map as another layer.

http://mpa.itc.it/radim/isle3.png
http://mpa.itc.it/radim/isle2.png
http://mpa.itc.it/radim/isle1.png

I don't know if it is good or not, we have to test it,
but it seems that it is more suitable in most cases and for
most users.

My suggestion is to leave d.vect and d.what.vect as it is
with all the complexity for experts and follow layer approach
(as described above) in user-friendly viewers.
I would appreciate, if you could try QGIS and tell me if it
does what you expect. If we find the QGIS approach good, the
same will be used for GRASS driver in OGR.

> 3) http://www.gdf-hannover.de/dassau/muell/v.extract_d.png
>
> Using v.extract with dissolve option
> v.extract -d in=soils out=soils3 type=area,centroid where='cat<20'
>
> -> areas with cat<20 are dissolved (very nice), but now areas with
> dublicated centroids exist. And areas cat='21 and 43' still exist but
> without centroid, as in example 2.

v.clean tool=rmdac # remove duplicate area centroids

This is right.

I'm not sure how it decides which one to keep, but I'd suspect the
area's cat might not be correct after this step.

The first found centroid is keept and all others are duplicate (removed).
The category will be correct, because boundaries are dissolved only
between areas with the same category.

Probably extract.c's dissolve needs to work on centroids somehow (right
now it's only boundaries).

Yes, as I said above. Let me know once you commit dissolve boundary patch
and I'll look at centroids. Or you can do it yourself,
Vect_build()
for each centroid
  if area < 0 // duplicate
     Vect_delete_line

see also grass51/vector/v.clean/rmdac.c

The same must be done in v.reclass

Radim

6

So it's clear: here's a picture to describe what I am trying to do.
(d.vect map=land cat=1-100000)

The holes/islands/areas-without-centroids are (real) islands in rivers.
I removed the rivers with 'v.extract -d list=1-100000', but the inner
areas remain as artifacts. (as from original shape file)

I am trying to extract a vector to tag as the coastline at 0 elevation
for a DEM interpolation. These left over islands can be very high up in
the mountains and make giant holes to sea level in the DEM.

'd.vect cat=1-100000' works great, I don't mean to change that.

What I need to do is remove those areas from the file altogether.

guesses at solutions:
. If there (is?) was a tool to add centroids to all areas/voids, I could
run 'v.extract -d' on it and be done?
. Or convert inner boundaries to lines..?
. Maybe dissolve needs to be smarter?

The QGIS polygon_0 layer should work well to show me where they are;
some of the areas are tiny and very hard to find when the rest of the
vector is shown. I would still have to use v.digit to remove them by
hand though if no other tool was available. (doing this can cause
v.digit to SegFault, BTW)

For more headaches, there is at least one spot with an island in a lake
on an larger island in a river.

thanks,
Hamish

(attachments)

extract_poly.png

7

On Wednesday 12 May 2004 10:54, Hamish wrote:

> > Probably extract.c's dissolve needs to work on centroids somehow
> > (right now it's only boundaries).
>
> Yes, as I said above. Let me know once you commit dissolve boundary
> patch and I'll look at centroids.

Ok, done.

I have done the centroids.

Radim

8

> > > Also, 'v.extract type=area list=1,2,3' keeps areas with missing
> > > centroids; how do you automatically add missing centroids so you
> > > can then get rid of the zombie areas with another [xxxxxxxx]?
>
> => 'v.extract -d list=1-50000'
>
> I thought 'v.clean tool=bpol' would do this, but it doesn't.
> Should it? Or should/is this better in another tool?

So again, any space closed by boundary, even if it does not have any
centroid is area. It is area without centroid, without category.
We can discuss, if d.vect should (and when) display those areas, se
more below. It is not a bug! It is not a bug! It is not a bug!

Can those lines/boundaries be picked by using 'v.select list=0' ?
Would/could that work?

In QGIS, one field means layer. If an area in GRASS has a category of
field 1, it is represented as feature in the layer polygon_1.
But it is different from GRASS, because user can see in the dialog,
that there is also layer polygon_0 (islands without category)
and he can add it to the map as another layer.

I don't know if it is good or not, we have to test it,
but it seems that it is more suitable in most cases and for
most users.

My suggestion is to leave d.vect and d.what.vect as it is
with all the complexity for experts and follow layer approach
(as described above) in user-friendly viewers.

yes.

I would appreciate, if you could try QGIS and tell me if it
does what you expect. If we find the QGIS approach good, the
same will be used for GRASS driver in OGR.

Ok, got it to compile after forcing the configure script to look in
/usr/include/qt3 [needed the libqt3-mt-dev debian package].

After installing, and starting qgis from within GRASS 5.7, I could plot
up the polygon_0 layer to highlight where the areas without centroids
that I needed to get rid of were (line legend -> red,50 wide). Then,
knowing where they were, I could go back to v.digit and remove them.
Success!

Two QGIS glitches (latest qgis CVS):
a) when zooming in a moderate amount (with the polygon_1 enabled) QGIS
would freeze my X-windows for 2-3 minutes at a time. Mouse worked,
everything else was frozen. After the wait it would redraw. Had to kill -9
it remotely a couple of times..

b) Some of the polygon_0 areas wouldn't show up. With a huge line width,
some would show up as a vertical bar, others blank. By playing around
with clicking on values in the attribute table to highlight the
individual polygons and changing the search radius for querying
by mouse [blind], I was finally able to track down the last few by
brute force. Yellow on white highlight for the attr table is really hard
to see..

Who is the polygon_2 layer and what does it do?

cheers,
Hamish

9

On Wednesday 12 May 2004 11:51, Hamish wrote:

The holes/islands/areas-without-centroids are (real) islands in rivers.
I removed the rivers with 'v.extract -d list=1-100000', but the inner
areas remain as artifacts. (as from original shape file)

What I need to do is remove those areas from the file altogether.

guesses at solutions:
. If there (is?) was a tool to add centroids to all areas/voids, I could
run 'v.extract -d' on it and be done?

Yes, islands are probably without category, assign a category (new centroid)
also to islands before v.extract: v.category type=area

For more headaches, there is at least one spot with an island in a lake
on an larger island in a river.

Should not change anything.

Radim

10

On Tuesday 18 May 2004 13:48, Hamish wrote:

> So again, any space closed by boundary, even if it does not have any
> centroid is area. It is area without centroid, without category.
> We can discuss, if d.vect should (and when) display those areas, se
> more below. It is not a bug! It is not a bug! It is not a bug!

Can those lines/boundaries be picked by using 'v.select list=0' ?
Would/could that work?

No.

Two QGIS glitches (latest qgis CVS):
a) when zooming in a moderate amount (with the polygon_1 enabled) QGIS
would freeze my X-windows for 2-3 minutes at a time. Mouse worked,
everything else was frozen. After the wait it would redraw. Had to kill -9
it remotely a couple of times..

Do you think that it is because drawing is slow (many lines) or
it is realy frozen? top?

Who is the polygon_2 layer and what does it do?

v.in.ogr writes the number of overlapping input areas as category of field 2,
usually errors in input data.

Radim

11

On Wednesday 19 May 2004 02:15, you wrote:

> > Two QGIS glitches (latest qgis CVS):
> > a) when zooming in a moderate amount (with the polygon_1 enabled)
> > QGIS would freeze my X-windows for 2-3 minutes at a time. Mouse
> > worked, everything else was frozen. After the wait it would redraw.
> > Had to kill -9 it remotely a couple of times..
>
> Do you think that it is because drawing is slow (many lines) or
> it is realy frozen? top?

X is really frozen for the time. Including xclock, gkrellm, NumLock
LED, keyboard.. Not sure if that is a bug in X or QGIS.
(Debian's X is pretty good though)

Running top (remotely) shows that XFree86 is using 99% of the CPU, and
QGIS is using 0% (QGIS is running tunneled over ssh, and the processor
of the machine it is on isn't under another load, and the network
traffic is low)

Even though there are complicated boundaries, it seems fast enough when
zoomed out.

I just tried with the huge 20m contour vector file for all of New
Zealand (lines only). It's slower (but X is fine) when zoomed out, but
fast when zoomed in. The other file I was trying (coastline + islands as
closed polygons) containing areas is fast when zoomed out and slow
(Xfree86) when zoomed in..

Also note in the attached screenshots that the entire screen was flooded
with the fill color when it finally rendered after a minute or so
(zoomed in). It's fine when zoomed out, but still 10sec of waiting at
that scale.

I think I know why it is slow, I could reproduce it.
When coordinates of the polygon are recalculated from map units to display
(canvas) units it may happen (if you zoom out) that some values overflow
integer maximum value. Result in the display is (known) mess:
http://mpa.itc.it/radim/qgis/poly-zoom-out.png (should be only filled by blue)

I guess that this also makes rendering (on X side) slow, because X has to render
very complex (wrong, more or less random) polygons.

Just now, I don't know how to solve this, because any pre-clipping would be slow?
Anybody has idea?

Radim

12

Hi

I think I know why it is slow, I could reproduce it.
When coordinates of the polygon are recalculated from map units to display
(canvas) units it may happen (if you zoom out) that some values overflow
integer maximum value. Result in the display is (known) mess:
http://mpa.itc.it/radim/qgis/poly-zoom-out.png (should be only filled by blue)

I guess that this also makes rendering (on X side) slow, because X has to render
very complex (wrong, more or less random) polygons.

Just now, I don't know how to solve this, because any pre-clipping would be slow?
Anybody has idea?

I hope it doesn't sound dumb but why not use doubles instead of integer to circumvent this?
It would need much more memory but we are on the save side then.

Ciao,
  Jens

13

On Wednesday 19 May 2004 12:46, Jens Oberender wrote:

Hi

> I think I know why it is slow, I could reproduce it.
> When coordinates of the polygon are recalculated from map units to
> display (canvas) units it may happen (if you zoom out) that some values
> overflow integer maximum value. Result in the display is (known) mess:
> http://mpa.itc.it/radim/qgis/poly-zoom-out.png (should be only filled by
> blue)
>
> I guess that this also makes rendering (on X side) slow, because X has to
> render very complex (wrong, more or less random) polygons.
>
> Just now, I don't know how to solve this, because any pre-clipping would
> be slow? Anybody has idea?

I hope it doesn't sound dumb but why not use doubles instead of integer to
circumvent this? It would need much more memory but we are on the save side
then.

XFillPolygon(display, d, gc, points, npoints, shape, mode)
             Display *display;
             Drawable d;
             GC gc;
             XPoint *points;
             int npoints;
             int shape;
             int mode;

typedef struct {
     short x, y;
} XPoint;

So not integer, only short.

Radim

14

On Monday 09 February 2004 16:37, Hamish wrote:

> > I've got a vector area which contains a bunch of islands and
> > polygons attached to the inside edge of the outer polygon. What
> > commands would one use to remove and 'flatten' the overlapping
> > areas?
> >
> >
> > eg
> >
> >
> > /-------+---+----------\
> >
> > | | X | |
> > | | X | <- area |
> > | |---| |
> > | |
> > | |=| <- island |
> >
> > \----------------------/
> >
> >
> > In the end I want a convex hull of the outer polygon as a single
> > entity, without losing resolution from a points->v.hull step.
>
> Sorry, I don't quite understand what you need, could you send me any
> more intuitive images/description?

my ascii art is poor..

see attached PNGs. (1 is from d.vect, 2 & 3 are from v.digit)

I have a vector file with an island as a big polygon, it has both the
coastline + river boundaries. The river polygons share a common boundary
with the coastline where they meet the ocean. Inland the river polygons
break off into topologic "islands", maybe because of a lake or a marsh
breaks the river.

I want a vector of only the outside coastline, with the rivers and
topo-islands removed.

By using v.extract I can get rid of the inner topo-islands by copying
out only the outer big islands' category number.

But I still can't get rid of the connected river polygons.

If I understand well, v.extract type=area new=1 list=1-10000 (categories of all areas)
with -d (Dissolve common boundaries) flag shoud do that, unfortunately
-d flag is disabled because there is a bug in the code which removes centroids.

Radim