12
$\begingroup$

We can have road networks of the countries in GIS shape file.For exmple from here. Normally traffic flows are proportionate to populations of the cities. enter image description here

I wonder how we can have the same in Mathematica?Is it possible at all?

$\endgroup$
4
  • 2
    $\begingroup$ Alex, this is an interesting question but it would help if you posted what you tried and where you got stuck. $\endgroup$
    – Verbeia
    Sep 17, 2013 at 12:20
  • $\begingroup$ Are you sure that information about traffic intensity is in this data? If not, the answer is: yes but provide the data to work on. I'm asking because I'm not able to find it, only the description that the road is state/country etc. $\endgroup$
    – Kuba
    Sep 19, 2013 at 11:47
  • $\begingroup$ @Kuba Hi great to see you again.I haven't said that the traffic data are the same as populations,I said that it is proportionate to each other.I mean as far as we are getting closer to the cities the traffic flows are getting bigger.The bigger the cities are, the traffic around them are crowded.This Map has been generated by a commercial software based on the traffic flow.I want to actually have it based on the population intensity as I asked this question. $\endgroup$
    – Alex
    Sep 19, 2013 at 12:04
  • $\begingroup$ @Kuba First I put the cities on the Map of the roads and tried to have the intersection of the roads which are the same as the cities coordinates as nodes.And the roads connecting cities as edges.But I kind of stuck!! $\endgroup$
    – Alex
    Sep 19, 2013 at 12:09

1 Answer 1

9
$\begingroup$

I chose to use the Australian data for patriotic reasons.

Going to the web site you suggested, it is straightforward to get the data for the Australian roads network. On unpacking the resulting ZIP file, there are four files: AUS_roads.dbf, AUS_roads.prj, AUS_roads.shp and AUS_roads.shx. Mathematica supports .shp import, so it is a simple matter of doing the following:

SetDirectory["~/Downloads/AUS_rds"];
initialdata = Import["AUS_roads.shp"]

enter image description here
(That little dot in the bottom right corner must be an anomaly. I can't find a matching inhabited island on Google Maps.)

The underlying data does not seem to include traffic flow data so I will have to make up the next bit. You would need to translate the underlying Graphics output into something you can use as a graphics visualisation. This rather complex piece of code does that:

rules = #1 -> #2 & @@@ (Join @@ (Partition[#, 2, 1] & /@ 
  (initialdata[[1, 2, 1, 2]] /. Line -> Sequence)));

What I'm doing here is picking out the pieces of the graphic that represent the red lines’ coordinates; I looked at the InputForm of initialdata to get the right Part specification. I used the replacement rule to convert this into lists of points instead of lines, and then use Partition to convert this in to pairs of points that are connected. I can then translate that into rules suitable for use in Graph or GraphPlot with the pure function #1 -> #2 that is applied at the row level (@@@).

There are about a quarter of a million pairs of coordinates, so here is a small sample:

In[28]:= Take[rules, 10]

Out[28]= {{142.162, -10.1571} -> {142.149, -10.1496}, {142.149, \
-10.1496} -> {142.138, -10.1469}, {142.138, -10.1469} -> {142.125, \
-10.1378}, {142.125, -10.1378} -> {142.117, -10.1383}, {142.117, \
-10.1383} -> {142.111, -10.1291}, {142.111, -10.1291} -> {142.103, \
-10.1257}, {142.103, -10.1257} -> {142.101, -10.1212}, {142.101, \
-10.1212} -> {142.102, -10.114}, {142.102, -10.114} -> {142.099, \
-10.1094}, {142.33, -10.1855} -> {142.324, -10.1873}}

It should be possible to respecify this as a Graph using something like the following:

Graph[Thread[Range[Length[rules]] -> Range[2, Length[rules] + 1]], 
 VertexCoordinates -> (Join[rules[[All, 1]], {rules[[-1, 2]]]})]

Here is an example where I only showed the first 101 points. Note the use of the undirected edge symbol. You could alternatively just turn off arrow heads.

Graph[Thread[Range[100] \[UndirectedEdge] Range[2, 101]], 
 VertexCoordinates -> (Join[rules[[1 ;; 100, 1]], {rules[[100, 2]]}])]

enter image description here

To get some ideas on how to get varying edge thicknesses, have a look at the answers to this question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.