Digitizing in GRASS

GRASS GRASS user
1

Hello All.

There is an option to break a line in the digitizer. But is there an
option to eliminate a node ?
Or, is there a function to eliminate end point nodes that are snapped?

Here is the problem in detail:

I am digitizing at a high zoom with the mouse and a backdrop image. So I
have to stop digitizing while I pan to the section.
After I finish I run Edit mode and snap all nodes, and label the areas.
When I run v.support and tell it to snap nodes to a given threshold, it
goes mad with the vector file. Because of this I would like to know if
there is manual method of eliminating unneeded nodes that are on the
same line.

Any suggestions?

Thanks,
Alejandro Imass

2

Alejandro Imass wrote:

Hello All.

There is an option to break a line in the digitizer. But is there an
option to eliminate a node ?
Or, is there a function to eliminate end point nodes that are snapped?

Here is the problem in detail:

I am digitizing at a high zoom with the mouse and a backdrop image. So I
have to stop digitizing while I pan to the section.
After I finish I run Edit mode and snap all nodes, and label the areas.
When I run v.support and tell it to snap nodes to a given threshold, it
goes mad with the vector file. Because of this I would like to know if
there is manual method of eliminating unneeded nodes that are on the
same line.

Any suggestions?

Thanks,
Alejandro Imass

Hi Alejandro

v.trim does this - as long as you set the threshold option to 0, then
it won't make any other changes. But you must copy the dig_att file from
the old map to the new, as v.trim doesn't (it seems) copy attributes.
It just builds an empty attribute file. Then once you've copied it,
run v.support. Voila - all excess nodes are gone!

There's also a recent module called v.build.polylines which also
is supposed to do this but I haven't been able to check it yet.
But there's no manual way of elliminating just one node while
in v.digit.

David

3

David D Gray wrote:

Hi Alejandro

v.trim does this - as long as you set the threshold option to 0, then
it won't make any other changes. But you must copy the dig_att file from
the old map to the new, as v.trim doesn't (it seems) copy attributes.
It just builds an empty attribute file. Then once you've copied it,
run v.support. Voila - all excess nodes are gone!

There's also a recent module called v.build.polylines which also
is supposed to do this but I haven't been able to check it yet.
But there's no manual way of elliminating just one node while
in v.digit.

I thought there was some secret keys (like line un-break) because I have
noticed some keys I press accidentally which are not in the menus do some
things, but I haven't looked closely to see what they do.

David

Thank you David.
Do you have any clue as to what could be happening in my second mail.

Sorry that my character sketch got all messed up. It's just two big polygons
sharing a middle line, and the middle line is cut by a smaller polygon.

Thanks a million.

4

Alejandro Imass wrote:

David D Gray wrote:

>
> Hi Alejandro
>
> v.trim does this - as long as you set the threshold option to 0, then
> it won't make any other changes. But you must copy the dig_att file from
> the old map to the new, as v.trim doesn't (it seems) copy attributes.
> It just builds an empty attribute file. Then once you've copied it,
> run v.support. Voila - all excess nodes are gone!
>
> There's also a recent module called v.build.polylines which also
> is supposed to do this but I haven't been able to check it yet.
> But there's no manual way of elliminating just one node while
> in v.digit.
>

I thought there was some secret keys (like line un-break) because I have
noticed some keys I press accidentally which are not in the menus do some
things, but I haven't looked closely to see what they do.

You may well be right, but I don't know if there are any.

Do you have any clue as to what could be happening in my second mail.

Sorry that my character sketch got all messed up. It's just two big polygons
sharing a middle line, and the middle line is cut by a smaller polygon.

Actually, I was unable to verify this, as I tried a map with the
topology you described, and it behaved just as expected. The
listings produced by the export modules are arcs not rings. The six
listings produced therefore form only part of a polygon each, the
so-called `arc-node' format, unless it is a simple island, of
course, but there is none in this topology. If v.out.ascii or
v.out.arc is producing actual closed rings there must be
something wrong with the original map.

David

5

David D Gray wrote:

so-called `arc-node' format, unless it is a simple island, of
course, but there is none in this topology. If v.out.ascii or
v.out.arc is producing actual closed rings there must be
something wrong with the original map.

Thanks David. This is correct and I verified this, it is indeed producing a
correct arc-node format.
(I suppose then that area labeling creates the centroids?)
I have the same problem as Justin Hickey describes, but I don't know how to deal
w/ it. I too was hacking my way last night at the ascii files, but there should
be an explanation, and it's probably a simple one.

I noticed that gen2shp does not use the lab file (I really don't know if it
supposed to), so I guess that the centroids are ignored. This is perhaps the
reason that the shape viewer I'm using interprets 6 polygons instead of 3 and
closes each arc separately by throwing a straight line between the nodes.

So now some more questions come to mind:

1) Is there a utility such as Justin Hickey's hack that figures out a closed
polygon taking the centroids from the area labels and generating the appropriate
ascii file?
2) Is there really a problem w/ gen2shp and it should take the lab file into
account (or perhaps the problem is mine and I did not understand the source
code) ? Maybe gen2shp was designed with (1) in mind ?!
3) I noticed a centroid calculator in the contrib dir of shapelib. Maybe gen2shp
expects shapelib to figure this out?
4) Is there a problem with the shape file viewer I'm using (MapServer) and IT
should figure out what the polygons should like? (But without a centroid I
consider this a tough job, unless it would use something like mentioned in (3)!)
5) Am I completely off course here ?

If (2) is true I would be more than happy to throw in my 2 cents and modify
gen2shp "...with a little help from my friends..." of course!
As I said before I'm so happy w/ GRASS that I would really like to contribute to
this project, in any way that I can.

Thanks to all.

Alejandro Imass

6

Alejandro Imass wrote:

If (2) is true I would be more than happy to throw in my 2 cents and modify
gen2shp "...with a little help from my friends..." of course!
As I said before I'm so happy w/ GRASS that I would really like to contribute to
this project, in any way that I can.

Well I found out w/ some friends that shape does not support topology that well and
that polygons are defined just a set of points.
So my best conclusion at the moment would be to suspect that gen2shp or shapelib
should resolve this (i.e. they should be able to read the lab file and close the
polygons correctly), because it seems that GRASS and the visualizer I'm using
(MapServer) are behaving as expected.

So maybe there is probably need to create a new version of gen2shp perhaps
gen2shpcent that would actually resolve the (or my) polygon dilema.

What do y'all think?

Best Regards,
Alejandro

7

Hello

I don't know if this will help you or not, but about a year ago I wrote a
module v.out.poly. If I remember correctly, this module's output is the same as
v.out.ascii except that it guarantees that each list of points is a complete
polygon instead of possibly being only an edge of a polygon, and the
coordinates are reversed.

I wrote this when we discovered that both GRASS and arc-info ouput area edges
and not complete polygons. We had a large polygon that was broken up into 3
edges but we needed to have a complete polygon as output. If the polygons are
broken up into edges, it is practically impossible to know which edges belong
to which polygon based on the output of v.out.ascii.

I never got around to submitting it to the source tree because it was a bit of
a hack to get our job done and I simply did not have a lot of time to clean it
up. Besides, I'm not a GIS person and I do not know if I did the proper thing
or not. I simply looked at GRASS function calls and data structures, found that
there was an area id and used it to output the points of the polygon. I copied
the v.out.ascii code and went from there. I don't even know if it is compatible
with version 5.

If you want to give it a try, I can send the code to you but I would strongly
recommend that you look through the code to see if it does what you expect.
It's not too long and is contained in a single file. I also put some comments
in to describe what the code does.

Again, I'm not sure if this will help you or not, but you are welcome to try it
if you want.

--
Sincerely,

Jazzman (a.k.a. Justin Hickey) e-mail: jhickey@hpcc.nectec.or.th
High Performance Computing Center
National Electronics and Computer Technology Center (NECTEC)
Bangkok, Thailand

People who think they know everything are very irritating to those
of us who do. ---Anonymous

Jazz and Trek Rule!!!

8

Justin Hickey wrote:

Hello

I don't know if this will help you or not, but about a year ago I wrote a
module v.out.poly. If I remember correctly, this module's output is the same as
v.out.ascii except that it guarantees that each list of points is a complete
polygon instead of possibly being only an edge of a polygon, and the
coordinates are reversed.

Mostly I have been plagued by the opposite problem where you
are given complete rings, and have to split them into arcs and
nodes, but it can be a problem this way round as well.

If you want to give it a try, I can send the code to you but I would strongly
recommend that you look through the code to see if it does what you expect.
It's not too long and is contained in a single file. I also put some comments
in to describe what the code does.

Yes. I would like to have a look at that

David

9

Justin Hickey wrote:

If you want to give it a try, I can send the code to you but I would strongly
recommend that you look through the code to see if it does what you expect.
It's not too long and is contained in a single file. I also put some comments
in to describe what the code does.

Again, I'm not sure if this will help you or not, but you are welcome to try it
if you want.

Yes. I would like to see this very much.

Thank You,
Alejandro

BTW. Take a look at my response to David and tell me if this is on the same lines
as what you were experiencing. It seems very much so to me, but I could be wrong.

10

David,

todays I tried your updated v.in.shape (CVS version
from yesterday) on some SHAPE-files. These shape files
are in quite bad condition and your new
v.in.shape creates a proper topology from it!

  CONGRATULATIONS!

In my original dataset there have been up to
11 lines parallel (or cutting each other sometimes),
your v.in.shape removed the rubbish perfectly.
So it was a hard test for your module :slight_smile:

Keep up the good work!

Kind regards
  
   Markus

--
Dipl.-Geogr. Markus Neteler * University of Hannover
Institute of Physical Geography and Landscape Ecology
Schneiderberg 50 * D-30167 Hannover * Germany
Tel: ++49-(0)511-762-4494 Fax: -3984

11

Hello to everybody
I just compiled GRASS5b7 as I was very keen on the newly now implemented G3D
modules. I tried s.vol.idw (as s.vol.rst is not yet implemented in the
official release) in order to run an interpolation on a 3 dimensional site
file (x|y|z|%w). However I ran into some troubles:
I couldn't import (s.in.ascii) my sites. Although my ascii file with the
site information looks like the new specified site format
(http://www.geog.uni-hannover.de/grass/gdp/html_grass5/programming/sites-api
/index.html), but GRASS couldn't read the file properly (warnings were
printed) and only the first line of the file was read. Do I need header
information in the ascii file, or does s.in.ascii not properly import ascii
site files?
And for the G3D part: will s.vol.idw produce a raster file that can be read
by the standard 2d GRASS modules? (such as d.rast)

Thanks for your support
Bernhard Sturm
---
University of Berne
Centre for Development and Environment
Hallerstrasse 12
3012 Bern

12

Hi Bernhard,

On Fri, May 19, 2000 at 10:37:16PM +0200, Bernhard Sturm wrote:

Hello to everybody
I just compiled GRASS5b7 as I was very keen on the newly now implemented G3D
modules. I tried s.vol.idw (as s.vol.rst is not yet implemented in the
official release) in order to run an interpolation on a 3 dimensional site
file (x|y|z|%w).

Note: s.vol.rst is not released yet by its authors.

However I ran into some troubles:
I couldn't import (s.in.ascii) my sites. Although my ascii file with the
site information looks like the new specified site format
(http://www.geog.uni-hannover.de/grass/gdp/html_grass5/programming/sites-api
/index.html), but GRASS couldn't read the file properly (warnings were
printed) and only the first line of the file was read. Do I need header
information in the ascii file, or does s.in.ascii not properly import ascii
site files?

Which warning did you get? s.in.ascii should accept ASCII files with a
structure as follows:

x1 y1 z1 attribute1
x2 y2 z2 attribute2
x3 y3 z3 attribute3
[...]

And for the G3D part: will s.vol.idw produce a raster file that can be read
by the standard 2d GRASS modules? (such as d.rast)

r3.to.rast was planned but never written (any volunteers?).

Hopefully forthcoming is r3.to.sites (under development), then you
could use s.to.rast. Note: s.to.rast should be updated with an
option which attribute to select (could be taken from s.surf.rst,
searching for volunteers as well...)

To display in 3D use r3.mkdspf to create isosurfaces and then r3.showdspf
to display these isosurfaces in 3D.

Jaro Hofierka wrote s.to.rast3 last week, find it in CVS (and forthcoming
GRASS 5 release). Use this to convert your 3D sites into grid3D. You
can grab the code from CVS web interface.

Best wishes

Markus

13

Hi Markus

Which warning did you get? s.in.ascii should accept ASCII files with a
structure as follows:

x1 y1 z1 attribute1
x2 y2 z2 attribute2
x3 y3 z3 attribute3
[...]

This works now fine, the problem was, that I wanted to import too many
attributes (I had also string attributes that I wanted to import, but I
dropped them, and now it's working alright)

But whenever I try to run s.vol.idw I get an error message telling me, that
the module checks for a missing PERMANENT/DEFAULT_WIND3 and then
PERMANENT/WIND3 file. I assume I have to generate these files when building
a new mapset database, but I couldn't figure out how and when those files
are created...

Any idea?

thanks for the help and (I am enjoing GRASS more and more, since I started
writting my own scripts :wink:

Bernhard Sturm
--
Department of Geography
Centre for Development and Environment
University of Berne
Hallerstrasse 12
3012 Bern - SWITZERLAND
www.cde.unibe.ch

14

On Sun, May 21, 2000 at 12:38:18PM +0200, Bernhard Sturm wrote:

But whenever I try to run s.vol.idw I get an error message telling me, that
the module checks for a missing PERMANENT/DEFAULT_WIND3 and then
PERMANENT/WIND3 file. I assume I have to generate these files when building
a new mapset database, but I couldn't figure out how and when those files
are created...

Any idea?

Yes, seems there's nothing building these if they don't already exist.
Buried in the documentation/APIs is the format (which is slightly
different than WIND/DEFAULT_WIND:

Most important is to have all the field names exactly right:

Example DEFAULT_WIND3:

Proj: 99
Zone: 0
North: 4452506.013229
South: 3392725.581724
East: 542351.332382
West: -374215.527298
Top: 10000.000
Bottom: -100.00
nofRows: 2120
nofCols: 1833
nofDepths: 100.00
e-w resol: 500.03647555
n-s resol: 499.89642996
t-b resol: 100.0

Example WIND3:

Proj: 99
Zone: 0
North: 10.00000000000000000000000000000000000000000000000000
South: -10.00000000000000000000000000000000000000000000000000
East: 10.00000000000000000000000000000000000000000000000000
West: -10.00000000000000000000000000000000000000000000000000
Top: 10.00000000000000000000000000000000000000000000000000
Bottom: 0.00000000000000000000000000000000000000000000000000
nofRows: 20
nofCols: 20
nofDepths: 10
e-w resol: 1.00000000000000000000000000000000000000000000000000
n-s resol: 1.00000000000000000000000000000000000000000000000000
t-b resol: 1.00000000000000000000000000000000000000000000000000

This stuff's still pretty experimental though, I guess.

--
¶ One·should·only·use·the·ASCII·character­set·when·compos­

» ing·email·messages.


15

On Sun, May 21, 2000 at 12:38:18PM +0200, Bernhard Sturm wrote:

Hi Markus

> Which warning did you get? s.in.ascii should accept ASCII files with a
> structure as follows:
>
> x1 y1 z1 attribute1
> x2 y2 z2 attribute2
> x3 y3 z3 attribute3
> [...]

This works now fine, the problem was, that I wanted to import too many
attributes (I had also string attributes that I wanted to import, but I
dropped them, and now it's working alright)

But whenever I try to run s.vol.idw I get an error message telling me, that
the module checks for a missing PERMANENT/DEFAULT_WIND3 and then
PERMANENT/WIND3 file. I assume I have to generate these files when building
a new mapset database, but I couldn't figure out how and when those files
are created...
Any idea?

Of course :slight_smile:

You need to write above files. Perhaps I write a script to
create them automatically from WIND and DEFAULT_WIND (or you
write it?). Structure is:

PERMANENT/WIND/ and PERMANENT/DEFAULT_WIND

---------cut here----------
Proj: 0
Zone: 0
North: 10.00000000000000000000000000000000000000000000000000
South: 0.00000000000000000000000000000000000000000000000000
East: 10.00000000000000000000000000000000000000000000000000
West: 0.00000000000000000000000000000000000000000000000000
Top: 10.00000000000000000000000000000000000000000000000000
Bottom: 0.00000000000000000000000000000000000000000000000000
nofRows: 10
nofCols: 10
nofDepths: 10
e-w resol: 1.00000000000000000000000000000000000000000000000000
n-s resol: 1.00000000000000000000000000000000000000000000000000
t-b resol: 1.00000000000000000000000000000000000000000000000000
------ cut again ----------------

Edit and change above file to your projection and coordinates.
Take care: The capital letters are required.

Then run g3.region to re-check it.

thanks for the help and (I am enjoing GRASS more and more, since I started
writting my own scripts :wink:

In case of interest for others, I will be happy to include your
script(s). Perhaps we start an own

"GRASS Scripts page"

on the server? Then scripts could be shared with others.

Best wishes

Markus

--
Dipl.-Geogr. Markus Neteler * University of Hannover
Institute of Physical Geography and Landscape Ecology
Schneiderberg 50 * D-30167 Hannover * Germany
Tel: ++49-(0)511-762-4494 Fax: -3984

16

Alejandro Imass wrote:

Justin Hickey wrote:

>
>
> Hope this helps.
>
> --
> Sincerely,
>

Thank you Justin.
Even though I've a lot of stuff on my mind I've spent the last few days trying to
determine an algorithim that figures out if a given centroid is in the area
(using the formula for centroid) of which adjacent arcs (arcs that are snapped at
some nodes) and then build arc/view type polygons (arcs are redundant in the area
edges) by concatenating these arcs. For example, a six arc ungenerate file from
GRASS would probably turn out to be a 3 polygon file by adding the adjacenta arcs
to each other arc.
My idea is(was?) to process the arc/info genrate output from GRASS, which is in
arc-node format and make the arc/view type ascii file from this one and then
feeding this new ascii file into gen2shp. Perhaps the centroids are not needed,
and I guess, from what you described, that your code does this in a simple form.
I sure hope so, because not bieng at all from this area, I have had to un-dust my
Calculus notes and books, even though being the GIS polygons defined by straight
edges (instead of curves) these areas can be quite "easily" ;-] (yeah right)
computed using triangles ! I guess.

Well anyway, I hope your code is much simpler than my way of thinking !

Thanks again and best regards,
Alejandro Imass

Hi Alejandro

Do you need a real centroid, ie. a mass point of the area, or would a
centrally located interior point do? If you just want a point for
labels, attribute assignments, perhaps the second option is what
you need. But there might be calculations you need to do that
require the actual centroid... The problem with using a real
centroid for labelling, attributes is that you can't be sure it
will not be inside an island, even if you take interior holes
into account - so you would have to move it anyway.

It's easier to calculate an any_old_point than a true centroid. You
just have to circulate the rings and find places where they intersect
the polygon's equator. A point can go in int1<->int2, int3<->int4 etc.
regions because these are `inside'. Of course it's more complicated
if there are holes, but you just have to circulate the holes' rings
as well to find further regions where points are forbidden. inside a
hole = outside the main polygon.

David

17

Justin Hickey wrote:

Alejandro

On May 22, 7:12am, Alejandro Imass wrote:
> Well anyway, I hope your code is much simpler than my way of thinking !

As you will see, the code actually uses the "area" data structure in GRASS. It
contains all the points for a particular area. At least it worked that way for
us. The polygon that was in three pieces using v.out.ascii was in one piece
using v.out.poly. I did not do any calculations such as calculating centroids
to combine arcs.

Hi

Isn't there a function that can calculate an area point _from_ GRASS if
the topo
information has been built? I'm sure I saw something but I can't
quite recall where.

Regards

David

18

David D Gray wrote:

Do you need a real centroid, ie. a mass point of the area, or would a
centrally located interior point do? If you just want a point for
labels, attribute assignments, perhaps the second option is what
you need. But there might be calculations you need to do that

Sorry for the misunderstanding, I was referring to the cuasi-centroid from the
labelling process (the second option you mentioned above).
It's just that I was thinking about using areas to figure out which arcs surrounded
this interior point, I got this idea from a centroid calculator I found in the contrib
dir of shapelib. It is probably not a good idea since as you can tell, I'm by no means
a GIS person.

OK, perhaps if I define the problem in more detail.... I will do the best I can.

So suppose that you only had three well defined areas. Two big areas, and a smaller
one in the middle, something like Ucrania, Rumania, and Moldavia. We have digitized
and labelled these areas, trimmed the file and come up with a total of 6 arcs, 4
nodes, and 3 center-or-inner points that define which country is which.

Now when you export this to any of two (GRASS and ARC/INFO) ascci formats, you will
end up with an arc-node file and a labels file which contains the centerpoints.
This is correct and perfect, but ARC/View, or rather MapServer which reads ARC/View
data, is kind'a dumb and what it interprets are 6 distinct polygons that will be
closed arbitarirly at at all end-nodes. This is because (and was explained to my by an
ARC/View user, so it could be wrong) ARC/View does not really support topology, so in
order to display this data correctly, it would need only three polygons defined as a
set of points, and it will close it at the end points if they are not snapped. What
you would call an Island in a package like GRASS.

So basically what I'm trying to do is create a little program that, known the center
points, can figure out which arcs surround it and can build these Islands which are
snapped to one another like jigsaw pieces. The ascii file from the GRASS export would
be fed to this progam and the resulting file would be in turn fed to gen2shp. Another,
probably better idea would be to add a parameter to gen2shp and have it read the .lab
file. The work of figuring out which arcs surround which centerpoint has to be done
either way.
It's just that building this into gen2shp, would probably be better becuase I guess a
lot of people already use this program.

BTW I think that Justin Hickey's program does exaclty this, in this case no further
work needs to be done, since it uses the area definitions in GRASS. I will test it
today and comment on this later.

What I cannot understand, and has me very curious, is why this has not been widely
needed before. It is hard for me to believe that only Justin and I have this
requirement !

Thanks, and hoping to hear your comments.

Best Regards,
Alejandro