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: