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 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.
>>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.
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.
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
- GRASS network not connected
- GRASS network connected
- GRASS network allocation
- Network Split with 1,2,3km from any school
- GRASS isonlines specific school
Pingback: NOW HIRING! | Information Technology Job
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!
What version of v.db.addcol are you using? The first step would be to check some examples from the docs:
http://grass.osgeo.org/grass65/manuals/v.db.addcol.html
Pingback: Back To GIS | Castells
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.
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/
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)
Sorry 😀 I don’t remember exactly
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.
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!
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.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!
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:
Olá, não achei essa opção no Qgis (v.in.ogr )
First things first… Portuguese is not a language that I know. Luckily for you this sentence is short and the structure is similar to spanish, so I gather that you did not find the option in Qgis.
Sort answer : Qgis != GRASS
Long Answer : You may be interested in the Qgis GRASS plugin (https://docs.qgis.org/2.6/en/docs/user_manual/grass_integration/grass_integration.html) and work from there.
Pingback: Barcelona trees | Castells
Pingback: Network Analysis QGIS. Barcelona distances | Castells
Thanks grreat blog post