I tried Sander's code for the initial post. (In Mathematica 12.) Everything worked fine up to the last step when it seemed to mess up plotting some of the routes. Maybe some change has broken the code. I found that it can be fixed by plotting the road segments individually. In the last line, replace

GeoPath[#2, "TravelPath"]

with

Map[GeoPath, #2]

I also found that the rounding strategy was not necessary. First, I experimented with just truncating the data to 5 decimal places, which seemed fine. Then when I produced a tidied up copy in a separate notebook, I forgot to include that step, and it was still fine. This was for the UK, with a grid of only 30x30. Whether it would be ok elsewhere, I don't know. This allowed simplification of the final steps, and I dispensed with the intermediate saving to disk. The grid for the UK came out not as regular as for France, with various empty spaces in the Scottish Highlands for example. My final code, for my town, was:

num = 30; country = "UnitedKingdom"; dest = 
 CityData[{"Bourne", "Lincolnshire", "UnitedKingdom"}, "Coordinates"];
a = 0.5; (*affects how the road width scales*)
cities = CityData[{All, country}];
citiesloc = CityData[#, "Coordinates"] & /@ cities;
citiesloc = DeleteMissing[citiesloc, 1, \[Infinity]];
reg = CountryData[country, "Polygon"];
reg = Identity @@ Identity @@ reg;
reg = First@MaximalBy[reg, Length];
region = MeshRegion[reg, Polygon[Range[Length[reg]]]];
rmf = RegionMember[region];
citiesloc = Select[citiesloc, rmf];
ccf = Nearest[citiesloc];
bounds = CoordinateBounds[reg];
n = {num, num};
pts = MapThread[Subdivide[#1[[1]], #1[[2]], #2] &, {bounds, n - 1}];
pts = Tuples[pts];
pts = DeleteDuplicates[First@*ccf /@ pts];
pts // Length;
Graphics[{{EdgeForm[Black], FaceForm[LightBlue], 
   Polygon[Reverse /@ reg]}, Point[Reverse /@ pts], 
  Style[Point[Reverse@dest], Red, PointSize[Large]]}]

Then...

routes = Map[{dest, #} &, pts];
td = Map[TravelDirections[GeoPosition /@ #] &, routes];
td = DeleteCases[td, $Failed];
tp = Map[#["TravelPath"][[1, 1]] &, td];
pieces = Map[Partition[DeleteDuplicates[Flatten[#, 1]], 2, 1] &, tp];
pieces = Flatten[pieces, 1];
tpieces = Tally[pieces];
tpieces = GatherBy[tpieces, Last];
tpieces = SortBy[tpieces, Part[#, 1, -1] &];
widths = tpieces[[All, 1, -1]];
tpieces = tpieces[[All, All, 1]];
GeoGraphics[
 MapThread[{Thickness[#1^a/1000], Map[GeoPath, #2]} &, {widths,
    tpieces}], GeoRange -> Entity["Country", country]]

And here is the result:

Realy cool! Beautiful post! I suggest you to use some compression in your file, they are too big. Maybe something like: Export["~/desktop/test.jpg",RandomImage[], "CompressionLevel"->0.2 ]

Also interesting is to look at the length of the route versus the direct distance:

copy=DeleteDuplicates[Flatten[#,1]]&/@alldata; (* merge data and delete duplicates*)
directlength=GeoDistance2@@@copy[[All,{1,-1}]]; (* calculate direct distance*)
routelength=BlockMap[GeoDistance2@@#&,#,2,1]&/@copy; (* calculate road distance *)
routelength=Total/@routelength; (* total it*)
ListPlot[{directlength,routelength}\[Transpose]/1000.0,Frame->True,AspectRatio->Automatic,Epilog->{Line[{{0,0},{1000,1000}}]},FrameLabel->(Style[#,14]&/@{"Direct distance [km]","Road distance [km]"})]
Histogram[routelength/directlength,Automatic,"PDF",Frame->True,FrameLabel->(Style[#,14]&/@{"\[CurlyPhi] = Road distance/Direct distance","PDF(\[CurlyPhi])"})]

where GeoDistance2 is a quick version of GeoDistance:

ClearAll[GeoDistance2]
GeoDistance2[{lat1_,lon1_},{lat2_,lon2_}]:=Module[{\[CurlyPhi]1,\[CurlyPhi]2,\[CapitalDelta]\[CurlyPhi],\[CapitalDelta]\[Lambda],a,c},
    {\[CurlyPhi]1,\[CurlyPhi]2}={lat1,lat2}Degree;
    \[CapitalDelta]\[CurlyPhi]=\[CurlyPhi]2-\[CurlyPhi]1;
    \[CapitalDelta]\[Lambda]=(lon2-lon1)Degree;
    a=Haversine[\[CapitalDelta]\[CurlyPhi]]+Cos[\[CurlyPhi]1]Cos[\[CurlyPhi]2]Haversine[\[CapitalDelta]\[Lambda]];
    c=2ArcTan[Sqrt[1-a],Sqrt[a]];
    6.3674446571 10^6c (* result in meters *)
]

I guess the average ratio describes how developed a country is/how mountainous it is/how many water bodies it has. A comparison for different countries would be very nice. But I don't want to kill the Wolfram servers...

Dear Sander,

what you (and Bernat) are doing is absolutely amazing! I will study your code in detail. Thanks a lot for sharing!

... you can just do this within 50 lines or so and it took me 30 minutes or so to write it! Imagine doing this in C, Python, Java, Swift, ....!

This is more than true!

Best wishes! -- Henrik

@Sander Huisman this is wonderful work! I liked your grid a lot too - good visual aid to understand the scales of computation. This, I bet, will get many readers and some might miss another nice visual explanation from the authors of original project. The explanation is a bit buried on their site. It also relates to what @Marco Thiel is suggesting, which is indeed an interesting direction of thinking. So here are a few images and quotes from the authors of the original project:

The idea behind the "Urban Mobility Fingerprint" is to calculate routes to a single destination (A) from many departure points (B). The points (B) are located on a radius around the destination. The size of the radius is defined by travel time with different transportation modes. All routes together then create the "Urban Mobility Fingerprint".

Once the routes are drawn, all starting points are rotated and align on to one ideal connection line. The direct line between these two dots is the ideal connection between the two locatoins. The length of the routes is scaled to fit the connecting ideal line. We call the graph that is created the "Street DNA". Since in reality street connections never follow a perfectly straight line deviations are unavoidable. The Street DNA graph shows exaclty these situations. The color-coding helps to identify in what compass direction routes are located on the map. The outline is created by connecting all locations (B) that are reached within the same travel time. These lines are called isochrones. The color of each route is defined by the compass direction it leads to, making it easier to relate routes within the fingerprint graph. We created the Urban Mobility Fingerprint and Street DNA graphs for a couple of cities around the world. The 13 examples we selected are supposed to resemble a lage diversity of culutral backrounds around the globe combined with objectivly compelling street networks. Next to different cultural backgrounds, we also looked for characteristic places and iconic street topography. Coastal locations like San Fransisco, Dubai's Palm Island or Manhanttan show the extreme deviations from the ideal connection we experience while travelling through our cities.

Creating these lines for the routes on a country-scale should be quite easy; just some rotation and scaling. I'll give it a try later. I chose a grid of cities rather than all cities so that it is more homogeneously distributed over the area, and thus not biasing certain areas.

Thanks for the nice comments!

Hi Marco,

Thanks for the kind words. Unfortunately I didn't save the times. But should be also possible. I might try tonight. Just to clarify: I'm also not using the distance now. I just color it by how often a piece of road is 'used' by all the routes. But it might be interesting to see it as a graph indeed. I will dig deeper later...

thanks!