[GRASS-user] Intersections on World Map

GRASS GRASS user
1

Hi,

I still get stucked on an issue using v.overlay. I imported
a shapefile
with the world continents (the one from ESRI). Projection
is Lat-Long.
Then I draw a line using v.in.ascii representing a latitude
e.g.

VERTI
L 2 1
-180 -80
180 -80
1 1

Now I use the following command to merge the intersection
between the
latitude (80 South) and the areas in the world continents
vector.

v.overlay ainput=$vector atype=line alayer=1 binput=World
btype=area \
blayer=1 output=$out operator=and

Works pretty well so far, until I reach latitude South 75+.
Intersections doesn't
match continents anymore (please see attached PNG file).

I would appreciate any hint.

Regards,
Martin

--
Martin Bley

(attachments)

map.png

2

hmmm
really strange...
did you verify the topology of the world map?

v.build map=$map

ciao

Il giorno gio, 01/11/2007 alle 16.11 +0100, Martin Bley ha scritto:

Hi,

I still get stucked on an issue using v.overlay. I imported
a shapefile
with the world continents (the one from ESRI). Projection
is Lat-Long.
Then I draw a line using v.in.ascii representing a latitude
e.g.

VERTI
L 2 1
-180 -80
180 -80
1 1

Now I use the following command to merge the intersection
between the
latitude (80 South) and the areas in the world continents
vector.

v.overlay ainput=$vector atype=line alayer=1 binput=World
btype=area \
blayer=1 output=$out operator=and

Works pretty well so far, until I reach latitude South 75+.
Intersections doesn't
match continents anymore (please see attached PNG file).

I would appreciate any hint.

Regards,
Martin

_______________________________________________
grassuser mailing list
grassuser@grass.itc.it
http://grass.itc.it/mailman/listinfo/grassuser

--
Ti prego di cercare di non inviarmi files .doc, .xls, .ppt, .dwg.
Preferisco formati liberi.
Please try to avoid to send me .doc, .xls, .ppt, .dwg files.
I prefer free formats.
http://it.wikipedia.org/wiki/Formato_aperto
http://en.wikipedia.org/wiki/Open_format

Ivan Marchesini
Department of Civil and Environmental Engineering
University of Perugia
Via G. Duranti 93/a
06125
Perugia (Italy)
Socio fondatore GFOSS "Geospatial Free and Open Source Software" http://www.gfoss.it
e-mail: marchesini@unipg.it
        ivan.marchesini@gmail.com
tel: +39(0)755853760
fax (university): +39(0)755853756
fax (home): +39(0)5782830887
jabber: geoivan73@jabber.org

3

Hi Ivan,
thanks for your answer

On Fri, Nov 02, 2007 at 12:51:50PM +0100, ivan marchesini wrote:

hmmm
really strange...
did you verify the topology of the world map?

v.build map=$map

GRASS 6.0.2 (World):~ > v.build map=World
Building topology ...
4010 primitives registered
Building areas: 100%
1983 areas built
1953 isles built
Attaching islands: 100%
Attaching centroids: 100%
Topology was built.
Number of nodes : 3980
Number of primitives: 4010
Number of points : 0
Number of lines : 0
Number of boundaries: 2042
Number of centroids : 1968
Number of areas : 1983
Number of isles : 1953
Number of areas without centroid : 15
GRASS 6.0.2 (World):~ >

don't know what to do with this information. What about the 15 "areas
without centroid"?

I just started a new v.overlay operation after building ithe topology
using v.build - same effect.

Thanks,
Martin

--
Martin Bley

4

Hi Martin,
grass has a topological management of vector maps...
grass uses "lines" to define terrain elements like streams, roads, and
so on (linear features)
grass uses "boundaries" to define areas...
boundaries are special kinds of lines that should be close on themselfe
or on other boundaries...
if you have a map containing only a boundary closed on itself then you
probably have:
Number of nodes : 1
Number of primitives: 1
Number of points : 0
Number of lines : 0
Number of boundaries: 1
Number of centroids : 0
Number of areas : 1
Number of isles : 1
Number of areas without centroid : 1

instead, if you put a centroid inside this area:
Number of nodes : 2
Number of primitives: 2
Number of points : 0
Number of lines : 0
Number of boundaries: 1
Number of centroids : 1
Number of areas : 1
Number of isles : 1

the centroid gives the category and the attributes to the area defined
by the sourronding boundaries...

so, usually, you can have an area without centroid only if this area
represents an hole (for example a lake in a vegetational map, since in
the water I consider there aren't plants)

in your case you have:
4010 primitives (2042 boundaries + 1968 centroids)
the boundaries define 1983 areas
but we have 1968 centroids
so there are probably 15 holes...

the holes can exists.. in effect... (for example lakes or internal seas)

this to explain the results of v.build...

but this doesn't explain the reason of the strange overlay...

why did you try to do all the intersection line for line...
why didn't you create a single "lines map" that contain all the
latitudes and then perform the overlay?

May be I'm missing something but I can't explain this result...
sorry...
can someone else try to give an explanation?

Ivan

Il giorno ven, 02/11/2007 alle 17.33 +0100, Martin Bley ha scritto:

Hi Ivan,
thanks for your answer

On Fri, Nov 02, 2007 at 12:51:50PM +0100, ivan marchesini wrote:
> hmmm
> really strange...
> did you verify the topology of the world map?
>
> v.build map=$map
GRASS 6.0.2 (World):~ > v.build map=World
Building topology ...
4010 primitives registered
Building areas: 100%
1983 areas built
1953 isles built
Attaching islands: 100%
Attaching centroids: 100%
Topology was built.
Number of nodes : 3980
Number of primitives: 4010
Number of points : 0
Number of lines : 0
Number of boundaries: 2042
Number of centroids : 1968
Number of areas : 1983
Number of isles : 1953
Number of areas without centroid : 15
GRASS 6.0.2 (World):~ >

don't know what to do with this information. What about the 15 "areas
without centroid"?

I just started a new v.overlay operation after building ithe topology
using v.build - same effect.

Thanks,
Martin

--
Ti prego di cercare di non inviarmi files .doc, .xls, .ppt, .dwg.
Preferisco formati liberi.
Please try to avoid to send me .doc, .xls, .ppt, .dwg files.
I prefer free formats.
http://it.wikipedia.org/wiki/Formato_aperto
http://en.wikipedia.org/wiki/Open_format

Ivan Marchesini
Department of Civil and Environmental Engineering
University of Perugia
Via G. Duranti 93/a
06125
Perugia (Italy)
Socio fondatore GFOSS "Geospatial Free and Open Source Software" http://www.gfoss.it
e-mail: marchesini@unipg.it
        ivan.marchesini@gmail.com
tel: +39(0)755853760
fax (university): +39(0)755853756
fax (home): +39(0)5782830887
jabber: geoivan73@jabber.org

5

you might try v.patch + v.clean instead of "v.overlay op=and"

?

Hamish

__________________________________________________
Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around
http://mail.yahoo.com

6

Hi Ivan,
hi folks,

On Sun, Nov 04, 2007 at 03:25:17PM +0100, ivan marchesini wrote:

[...]
this to explain the results of v.build...

thanks a lot for your detailed explanation.

but this doesn't explain the reason of the strange overlay...

I suspect, v.overlay using operator "AND" is buggy. I found a bug report
(ID #5427) describing this behavior, which has still open status. Using
operator "NOT", v.overlay will work properly. In my case that solves the
issue.

why did you try to do all the intersection line for line...
why didn't you create a single "lines map" that contain all the
latitudes and then perform the overlay?

Well, I'm new to GRASS and I did it this way, because don't know to do
better this time - I'll work on it :wink:
My task is the following: I would like to figure out for every single
latitude, how much of it is covered by ocean surface and how much is
covered by a continent. So, I want a table like this

Lat|Ocean length (meters)|total length (meters)
[...]
10N|???|39470.171
[...]

Total length I calculated using a perl script. To get the "Ocean
Length", I export the vector maps with the intersections (using points),
afterwards I use another perl script, to calc length of the segments and
finally the sum of all segments.

That's why I didn't use a single "lines map" - just to keep things
clearly arranged.

My next excercise will be connection those segments lengths to the GRASS
database.

Regards,
Martin

7

Martin Bley wrote:

My task is the following: I would like to figure out for every single
latitude, how much of it is covered by ocean surface and how much is
covered by a continent. So, I want a table like this

Lat|Ocean length (meters)|total length (meters)
[...]
10N|???|39470.171
[...]

^^ I think you mean km

Total length I calculated using a perl script. To get the "Ocean
Length", I export the vector maps with the intersections (using points),
afterwards I use another perl script, to calc length of the segments and
finally the sum of all segments.

Attached is a shell script that calculates this from the ETOPO2 raster map.
Besides showing off that dataset it demonstrates the r.reclass, r.profile, and
d.graph modules.

See the GRASS Newsletter Vol 1, July 2004 for details on the ETOPO2 dataset.

It isn't nearly as efficient or accurate as doing it with vectors+perl, but
it's done by the time you get a coffee.

I thought it interesting to look at slicing by longitude, so it does that too.

enjoy,
Hamish

      ____________________________________________________________________________________
Be a better pen pal.
Text or chat with friends inside Yahoo! Mail. See how. http://overview.mail.yahoo.com/

(attachments)

ocean_perc_by_lat.sh (2.77 KB)