plot_knuth_miles.py 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142
  1. """
  2. ===========
  3. Knuth Miles
  4. ===========
  5. `miles_graph()` returns an undirected graph over 128 US cities. The
  6. cities each have location and population data. The edges are labeled with the
  7. distance between the two cities.
  8. This example is described in Section 1.1 of
  9. Donald E. Knuth, "The Stanford GraphBase: A Platform for Combinatorial
  10. Computing", ACM Press, New York, 1993.
  11. http://www-cs-faculty.stanford.edu/~knuth/sgb.html
  12. The data file can be found at:
  13. - https://github.com/networkx/networkx/blob/main/examples/drawing/knuth_miles.txt.gz
  14. """
  15. import gzip
  16. import re
  17. # Ignore any warnings related to downloading shpfiles with cartopy
  18. import warnings
  19. warnings.simplefilter("ignore")
  20. import numpy as np
  21. import matplotlib.pyplot as plt
  22. import networkx as nx
  23. def miles_graph():
  24. """Return the cites example graph in miles_dat.txt
  25. from the Stanford GraphBase.
  26. """
  27. # open file miles_dat.txt.gz (or miles_dat.txt)
  28. fh = gzip.open("knuth_miles.txt.gz", "r")
  29. G = nx.Graph()
  30. G.position = {}
  31. G.population = {}
  32. cities = []
  33. for line in fh.readlines():
  34. line = line.decode()
  35. if line.startswith("*"): # skip comments
  36. continue
  37. numfind = re.compile(r"^\d+")
  38. if numfind.match(line): # this line is distances
  39. dist = line.split()
  40. for d in dist:
  41. G.add_edge(city, cities[i], weight=int(d))
  42. i = i + 1
  43. else: # this line is a city, position, population
  44. i = 1
  45. (city, coordpop) = line.split("[")
  46. cities.insert(0, city)
  47. (coord, pop) = coordpop.split("]")
  48. (y, x) = coord.split(",")
  49. G.add_node(city)
  50. # assign position - Convert string to lat/long
  51. G.position[city] = (-float(x) / 100, float(y) / 100)
  52. G.population[city] = float(pop) / 1000
  53. return G
  54. G = miles_graph()
  55. print("Loaded miles_dat.txt containing 128 cities.")
  56. print(G)
  57. # make new graph of cites, edge if less then 300 miles between them
  58. H = nx.Graph()
  59. for v in G:
  60. H.add_node(v)
  61. for u, v, d in G.edges(data=True):
  62. if d["weight"] < 300:
  63. H.add_edge(u, v)
  64. # draw with matplotlib/pylab
  65. fig = plt.figure(figsize=(8, 6))
  66. # nodes colored by degree sized by population
  67. node_color = [float(H.degree(v)) for v in H]
  68. # Use cartopy to provide a backdrop for the visualization
  69. try:
  70. import cartopy.crs as ccrs
  71. import cartopy.io.shapereader as shpreader
  72. ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal(), frameon=False)
  73. ax.set_extent([-125, -66.5, 20, 50], ccrs.Geodetic())
  74. # Add map of countries & US states as a backdrop
  75. for shapename in ("admin_1_states_provinces_lakes_shp", "admin_0_countries"):
  76. shp = shpreader.natural_earth(
  77. resolution="110m", category="cultural", name=shapename
  78. )
  79. ax.add_geometries(
  80. shpreader.Reader(shp).geometries(),
  81. ccrs.PlateCarree(),
  82. facecolor="none",
  83. edgecolor="k",
  84. )
  85. # NOTE: When using cartopy, use matplotlib directly rather than nx.draw
  86. # to take advantage of the cartopy transforms
  87. ax.scatter(
  88. *np.array(list(G.position.values())).T,
  89. s=[G.population[v] for v in H],
  90. c=node_color,
  91. transform=ccrs.PlateCarree(),
  92. zorder=100, # Ensure nodes lie on top of edges/state lines
  93. )
  94. # Plot edges between the cities
  95. for edge in H.edges():
  96. edge_coords = np.array([G.position[v] for v in edge])
  97. ax.plot(
  98. edge_coords[:, 0],
  99. edge_coords[:, 1],
  100. transform=ccrs.PlateCarree(),
  101. linewidth=0.75,
  102. color="k",
  103. )
  104. except ImportError:
  105. # If cartopy is unavailable, the backdrop for the plot will be blank;
  106. # though you should still be able to discern the general shape of the US
  107. # from graph nodes and edges!
  108. nx.draw(
  109. H,
  110. G.position,
  111. node_size=[G.population[v] for v in H],
  112. node_color=node_color,
  113. with_labels=False,
  114. )
  115. plt.show()