2
$\begingroup$

The goal of this code is to identify and graph a "circle of cities" about a center city.

cityCenter = Interpreter["City"]["St. Louis"]
stateCenter = Interpreter["USState"]["Missouri"]
nearbyStates = 
 Join[AdministrativeDivisionData[Interpreter["USState"][stateCenter], 
   "BorderingStates"], {stateCenter}]
distanceToCenter[s_] := 
 QuantityMagnitude[
  TravelDistance[Interpreter["City"][cityCenter], 
   Interpreter["City"][s]]]
cityList = 
  Sort[Flatten[
    Table[CityData[{All, nearbyStates[[k]][[2]][[1]], 
       "UnitedStates"}], {k, 1, Length[nearbyStates]}]]];
table1 = Table[{cityList[[k]], distanceToCenter[cityList[[k]]]}, {k, 
    Length[cityList]}];
radius = 50;
tolerance = 0.05;

table2 = Join[{Interpreter["City"][cityCenter]}, 
   Select[table1, 
    radius*(1 - tolerance) <= #[[2]] <= radius*(1 + tolerance) &]];

graph1 = GeoListPlot[(First /@ table2), PlotMarkers -> Point, 
   GeoLabels -> Automatic, GeoRangePadding -> Scaled[0.5], 
   ImageSize -> Medium, GeoStyling["StreetMapLabels"]];
Drop[(First /@ table2), 1]
Length@(First /@ table2)
graph2 = GeoGraphics[GeoMarker[Interpreter["City"][cityCenter]]];
Show[graph1, graph2]

This seems unnecessarily slow (i.e., there is something better). What's wrong? Is it the TravelDistance?

$\endgroup$
  • $\begingroup$ The poor performance is all the Interpreter[...] stuff and getting the CityData over the network. I even get a timeout error: EntityValue::timeout: A network operation for EntityValue timed out. Please try again later. The variable cityList is huge - I ended up with 6469 cities for Length[cityList] - most of the time spent in a runtime profile was just downloading the data. $\endgroup$ – flinty Jul 20 at 12:59
  • $\begingroup$ It runs as mine as-is. Odd! $\endgroup$ – Veritas Lux Jul 20 at 20:36
  • 1
    $\begingroup$ As a side note, it's possible to run your own instance of Graphhopper and query that without any restrictions or network issues (after a long setup time to build the graph). Graphhopper is almost certainly what WRI are using under the hood. It's a pity it seems impossible to point TravelDistance etc to a different instance. $\endgroup$ – Carl Lange Sep 20 at 8:35
2
$\begingroup$

As flinty said, most of the time is required to download data. There are a couple of changes that will help that problem: 1) get CityData only for states near cityCenter, and 2) use TravelDistanceList instead of TravelDistance.

If we test fewer cities, the process will be faster. We can limit the number of cities (and make fewer calls to CityData) by selecting the states with boundaries close to cityCenter. Here's how to to select the states.

nearbyStates = Select[
   Join[AdministrativeDivisionData[stateCenter,
     "BorderingStates"], {stateCenter}], 
   QuantityMagnitude@
      GeoDistance[cityCenter, #, DistanceFunction -> "Boundary"] < 
     radius*(1 + tolerance) &];

For St. Louis and radius = 50, we call CityData for two states instead of nine.

Next, we need to get distances more efficiently because network calls are the slowest part of the task. TravelDistance needs 1 network call for each distance it finds. For me, finding the distances ran for more than 15 minutes.

Instead, a call to TravelDistanceList returns the all the distances between every pair of locations in a list. We can get many distances with one network call. I found the task was complete in less than 4 minutes.

However, TravelDistanceList doesn't accept long lists of locations (250 seems to work), so the list of cities must be grouped into "edible" chunks. TravelDistanceList returns distances between pairs of locations (saved as distList), but we need only the odd-numbered results. Combine cityList with the odd-numbered distances, and group each city and its distance as table1.

Here's the code to limit nearbyStates and use TravelDistanceList. I've simplified calls to Interpreter[...], and changed the code for cityList and table1.

cityCenter = Interpreter["City"]["St. Louis"];
stateCenter = 
  cityCenter[EntityProperty["City", "AdministrativeDivision"]];
radius = 50;
tolerance = 0.05;
(*remove states if cityCenter is too far from a state's border*)
nearbyStates = Select[
   Join[AdministrativeDivisionData[stateCenter,
     "BorderingStates"], {stateCenter}], 
   QuantityMagnitude@
      GeoDistance[cityCenter, #, DistanceFunction -> "Boundary"] < 
     radius*(1 + tolerance) &];
(*distanceToCenter isn't needed, but it's useful for checking results*)
distanceToCenter[s_] := 
   QuantityMagnitude[TravelDistance[cityCenter, s]]
cityList = Sort[Flatten[
    CityData[{All, ##}] & @@@ 
     EntityValue[nearbyStates, "CanonicalName"]]];
distList = QuantityMagnitude[
   TravelDistanceList /@ Partition[
     Riffle[ConstantArray[cityCenter, Length[cityList]], cityList],
     UpTo[250]]
   ];
table1 = Partition[
   Riffle[cityList, 
    Flatten[#[[Range[1, Length[#], 2]]] & /@ distList]],
   2];
table2 = Join[{cityCenter}, 
   Select[table1, 
    radius*(1 - tolerance) <= #[[2]] <= radius*(1 + tolerance) &]];
graph1 = GeoListPlot[First /@ Rest[table2], PlotMarkers -> Point,
      GeoLabels -> Automatic, GeoRangePadding -> Scaled[0.5],
      ImageSize -> Medium, GeoBackground -> GeoStyling["StreetMap"]];
graph2 = GeoGraphics[GeoMarker[cityCenter]];
Show[graph1, graph2]
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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