Basic Network Analysis With GRASS

I performed a little network analysis with ArcGIS on the Schools in Enschede. Given a shapefile with roads and the school points I obtained the service areas at 1km and 2km from the schools. Can I obtain the same with GRASS?

First Notes

For those who never heard of GRASS (If you know, you can skip to the fun part).
GRASS is an impressive GIS package with a lot of functionalities. It is designed as an environment where different tools execute GIS computations.
By default all the tools can be executed directly from GRASS command line:
tool op1=val1 op2=val2
or by a GUI:
tool will launch the graphical interface

In this post I always use the command line tool with all the needed options set, but you can perform the same using the GUI.

Data Data!

To import the data:
v.in.ogr -o dsn=NWB_Enschede.shp output=roads
v.in.ogr -o dsn=Schools_corrected.shp output=schools

Notice the -o flag. This is needed because the projection description did not match exactly, nevertheless it was the same. When GRASS finds mismatching projection it shows the PROJ_INFO of both projections (file and project), a little comparison was enough to determine that the projection was apparently the same, so in this case it’s safe to override the check.

While inspecting the data I encountered the first “surprise”. In ArcGIS the road layer has an attribute called “Shape Length” that I use for the network analysis. In this case, there’s no such attribute (automagic by ArcGIS). In GRASS I need to add and populate a LENGTH attribute.

Add a new attribute of type Double:
v.db.addcol map=roads columns=LENGTH Double

Populate the attribute with the length of the road:
v.to.db map=roads option=length units=meters columns=LENGTH

Creating the Network

Everything needed is under v.net

First connect the school points to the road network (this adds new lines):
v.net input=roads points=schools output=network operation=connect thresh=500

GRASS network not connected


GRASS network connected

GRASS network connected

GRASS created lines to fully connect the schools to the road network. The schools and roads are joined in a layer called network:
v.category network op=report

Layer/table: 1/network
type       count        min        max
point          0          0          0
line        7132          1       7004
boundary       0          0          0
centroid       0          0          0
area           0          0          0
all         7132          1       7004
Layer: 2
type       count        min        max
point         64          1         64
line           0          0          0
boundary       0          0          0
centroid       0          0          0
area           0          0          0
all           64          1         64

On this result, and with a little bit of inspection, is easy to realize that layer 2 is not assigned to any database (with various information of the schools). But this information is already on our workspace so I’ll connect it. Connecting the layer to the information can be useful in the future.
v.db.connect map=network table=schools layer=2

Now we have the network prepared, let’s get the party started. The network utilities provided by GRASS are:

  • Network Maintenance (v.net)
  • Allocate Subnets (v.net.alloc)
  • Split net (v.net.iso)
  • Shortest Path (v.net.path)
  • Visibility Network (v.net.visibility)
  • Steiner Tree (v.net.steiner)
  • Traveling salesman analysis (v.net.salesman)

For this case I’m only interested in v.net.alloc and v.net.iso.
The first one creates sub networks according to nearest centre. The second one can split the network in various distances.

>>Network Allocation

First, let’s see how to allocate the network:
v.net.alloc --overwrite --verbose input=network output=network_alloc ccats=1-63
The --overwrite, input and ouptut are self explanatory. The ccats sets which schools are used in the allocation, the range 1-63 contains all the schools.

Using v.category again grass shows the different categories on the newly created layer.
The result looks boring in the map view. That’s because there is no color assigned to each class:
d.vect -c map=network_alloc
-c stands for “ramdomize randomize colors”.

Now I can see that the allocation worked, the different allocated zones have different colors.

GRASS network allocation in Rainbow mode 😛

>>Split net by cost

I’m interested to see who is more lucky. And in this situation I define luckiness as a function of closeness to school.

  • Less than 1km: You are very lucky! Sleeep!
  • 1km~2km: Not the best, not the worst
  • 2km~3km: Not the worst, not the best
  • More than 3km: I’m sorry, wake up at 6

Execute the v.net.iso module with 1000m,2000m and 3000m thresholds:
v.net.iso --overwrite input=network output=network_isolines ccats=1-64 costs=1000,2000,3000 afcolumn=LENGTH abcolumn=LENGTH

The ccats (again) represents the schools to use in the analysis.
costs is the isolines to separate (1km,2km,3km).
afcolumn and abcolumn is the dataset column for forward and backward cost.

GRASS v.net.iso Isolines

Network Split with 1,2,3km from any school

But this is a general view, I want to know my exact situation. I’m only interested in my school:
v.net.iso --overwrite input=network output=network_isolines_specific ccats=5 costs=1000,2000,3000 afcolumn=LENGTH abcolumn=LENGTH

d.vect -c map=network_isolines_specific

d.vect map=schools type=point where=cat=5 fcolor=9:157:9 icon=basic/circle size=10

First, specific analysis. Second, colorize the network. Third, show only the specific school on the map.

GRASS isonlines specific school

GRASS isonlines specific school

I know I know. Green usually means better, but in this case means less lucky!. That’s the result of the random colors.

Conclusion

I can obtain the same basic result in ArcGIS and GRASS.
It’s just a basic example, in other exercises with ArcGIS I performed a Location-Allocation analysis with Demand points and Facilities to properly decide which new facilities are more suitable, I’m not really sure if that can be done directly with GRASS, but I’m planning to check it out 🙂


References

GIS Pages, the University of Trento ⇒GO
v.net.alloc source ⇒GO
v.net.iso source ⇒GO
v.net.salesman source ⇒GO


Advertisements

15 Comments

Filed under gis

15 responses to “Basic Network Analysis With GRASS

  1. Pingback: NOW HIRING! | Information Technology Job

  2. Robert Capps

    Hey, thanks for the post!

    When I run this line:

    v.db.addcol map=roads columns=LENGTH Double

    I get an INVALID SYNTAX error at ‘columns.’ Any idea why?

    Thanks!

  3. Pingback: Back To GIS | Castells

  4. D G Rossiter

    So where are the enschede data files? My children went to some of those schools… I am trying to reproduce your analysis in QGIS with the GRASS algorithms.

    • kxtells

      HI 😀

      Well, that was quite some time ago and it was an exercise from the ITC faculty, so they provided the data to me.

      But I remember Enschede having a good background with an amazing group of people working for data openness.

      If I were you, my first stop would be:
      http://opendata.enschede.nl/

  5. D G Rossiter

    Thanks, I will check with ITC faculty — do you remember which instructor? There are no streets or schools in opendata.enschede.nl (yes I do speak Dutch)

  6. D G Rossiter

    OK I will follow up with my comments and if I work this example into a more detailed tutorial I will check back with you.

  7. anthom

    Thanks very much for the post.
    I performed exact the same by using GUI (as I use windows and there is no command line). After connecting the road network (based on OSM Berlin) with 45 points (geocoded service stations), I performed the allocation function. However, some points were allocated to the networks of others (every point has its own cat) – maybe similar to your green networks in the centre of your map – and some roads were “lost”. When I performed the iso function, for some of the 45 points a network split had been done and for others not. What am I doing wrong?
    Thanks!

    • kxtells

      Hi! I am glad to see that this post helped so many people :-D.

      I’ll try to guess (because that’s the only thing I can do, make an educated guess..)

      To the point. The -c option (randomize colors) is documented, and it tells us that it will generate colors randomly.

      It was okay for my example (that was not the main point of the post). But when they say random, you should always ask yourself… how random is it? How many colors will be generated?. What I’m pointig at is: The categories may be well allocated, but the colors are far too few to represent a difference.

      When I check the source code of lines.c, I can see a palette array with the colors: light gray, dark gray, bright red, dark red, bright green, dark green… and so on up to 16 colors , it is suspiciously similar to the colors in the map 😀

      I had 65 schools, and you have 45 service stations. Some colors will be repeated.

      So, I say, first double check the categories information on your generated dataset. Because they may be set correctly and it is simply a visualization problem.

      If that’s the case, you can generate your own palette and feed it to d.vec with the -a option. Or to make things easier, export the layer and generate the coloring with Qgis that I think is better suited, or at least more user friendly, for that task.

  8. anthom

    Thanks for your help!
    After checking the data, I suppose, it’s not a visualization problem, but a double-/triple-cat-problem. The new “road-servicestation-net” contains around 114,000 objects, but listed 113,900 different cats (cat-numbers) only. The created lines to the service stations have the same cat-number as the next road. When performing the alloc function, some lines (both new lines and “old” ones) are allocated to their own cats which differ from the remaining road network around the service stations. I don’t know how to manage this problem. Do you have any idea?
    Thanks!

    • kxtells

      Wow, sorry. Right now from what you say and without the data in front of me is difficult to say.
      But it may be quite an interesting problem.

      Keep me posted with:

  9. Olá, não achei essa opção no Qgis (v.in.ogr )

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s