Thank you!

One final question: for the parks on the outskirts of town should I just follow the roads? Swimming across the longest river in North America and running across an airport are not very practical. However, we will be on foot so some of the parks in town might be able to get to by cutting across some blocks/empty lots.

Thanks again for all the help.

Patrick

Hi,

I have just had another idea which is certainly not doable with WolframAlpha alone. I found this website, which lists city parks in Great Falls. We can import the website

url = Import["http://www.greatfallsmt.net/recreation/city-parks", "Data"];

and then extract a list of all parks.

parks = url[[1, 2, 1, 1, 2]][[All, 1]]

This gives a rather long list:

{"Anaconda Hills Golf Course", "BelView", "Bloomingdale", "Boston 
Heights", "Broadwater Bay", "Carter", "Charles Russell", "Chowen 
Springs", "Clara", "Community Hall", "Dudley Anderson", "Eagle Falls 
Golf Club", "Eagles Crossing", "Elks Riverside", "Fox Hollow", 
"Garden Home", "Gibson", "Grande Vista", "Graybill", "Heren", 
"Highland", "Horizon", "Hylande Heights", "Jaycee (Gallatin)", 
"Kiwanis - North", "Kiwanis - West", "Kranz", "Lion's", "Madison", 
"Margaret", "Meadowlark", "Memorial", "Mitchell Pool", "Montana", 
"Morony", "Multi Sports Complex", "Noah's Ark", "Oddfellow", 
"Overlook", "Pinski", "Rhodes", "Riverview", "Roosevelt", 
"Sacajawea", "Sandhills", "Sight & Sound", "Skyline", "Skyline 
Optimist", "Sunnyside", "Sunrise", "Tourist", "Valleyview", "Verde", 
"Vereran's Memorial", "Wadsworth", "Warden", "West Bank", "West 
Viaduct", "Westwood", "Whittier"}

In order to get the geo positions of these parks we use Google Maps. First we need to generate a query string; I'll show it for the first park:

StringJoin["http://maps.google.com/maps/api/geocode/json?address=", StringReplace[parks[[1]], " " -> "+"], "+Great+Falls", ",+Montana&sensor=false"]

This gives:

"http://maps.google.com/maps/api/geocode/json?address=Anaconda+Hills+Golf+Course+Great+Falls,+Montana&sensor=false"

We then need to extract the coordinates from the website. This would have been easy had I used:

Import["http://maps.google.com/maps/api/geocode/json?address=Anaconda+\Hills+Golf+Course,+Montana&sensor=false", "JSON"]

But when I first tried it I did not use the "JSON" option. I extracted the coordinates for all parks using:

coordsallparks = 
 DeleteDuplicates[
  Quiet[GeoPosition[{ToExpression@
        StringTake[
         StringSplit[
           StringSplit[StringSplit[#, "location"][[2]], "lat"][[2]], 
           ","][[1]], -10], 
       ToExpression@
        StringTake[
         StringSplit[
           StringSplit[StringSplit[#, "location"][[2]], "lat"][[2]], 
           ","][[2]], -27 ;; -15]}] & /@ (Import[#] & /@ (StringJoin[
          "http://maps.google.com/maps/api/geocode/json?address=", 
          StringReplace[#, " " -> "+"], "+Great+Falls", 
          ",+Montana&sensor=false"] & /@ parks))]]

instead. I then deleted the failed requests and one "outlier".

parkpos = Select[coordsallparks, (NumberQ[#[[1, 1]]] && NumberQ[#[[1, 2]]] && #[[1, 1]] < 50) &]

This gives:

{GeoPosition[{47.4942, -111.283}], GeoPosition[{47.4731, -111.341}], 
 GeoPosition[{47.5158, -111.34}], GeoPosition[{47.5059, -111.249}], 
 GeoPosition[{47.4955, -111.309}], GeoPosition[{47.7854, -111.334}], 
 GeoPosition[{47.5102, -111.286}], GeoPosition[{47.4981, -111.279}], 
 GeoPosition[{47.4951, -111.231}], GeoPosition[{47.4963, -111.248}], 
 GeoPosition[{47.4763, -111.324}], GeoPosition[{47.4653, -111.235}], 
 GeoPosition[{47.4735, -111.333}], GeoPosition[{47.5101, -111.23}], 
 GeoPosition[{47.4611, -111.256}], GeoPosition[{47.5325, -111.314}], 
 GeoPosition[{47.4865, -111.303}], GeoPosition[{47.5113, -111.267}], 
 GeoPosition[{47.5097, -111.325}], GeoPosition[{47.5014, -111.279}], 
 GeoPosition[{47.4955, -111.259}], GeoPosition[{47.5494, -111.499}], 
 GeoPosition[{47.582, -111.063}], GeoPosition[{47.4926, -111.307}], 
 GeoPosition[{47.492, -111.306}], GeoPosition[{47.4987, -111.226}], 
 GeoPosition[{47.5073, -111.326}], GeoPosition[{47.5304, -111.307}], 
 GeoPosition[{47.5483, -111.505}], GeoPosition[{47.5333, -111.295}], 
 GeoPosition[{47.4889, -111.267}], GeoPosition[{47.5341, -111.301}], 
 GeoPosition[{47.5373, -111.302}], GeoPosition[{47.4904, -111.282}], 
 GeoPosition[{47.502, -111.222}], GeoPosition[{47.5216, -111.333}], 
 GeoPosition[{47.4918, -111.305}], GeoPosition[{47.5129, -111.312}], 
 GeoPosition[{47.5065, -111.341}], GeoPosition[{47.5088, -111.293}]}

If we assume that these are the coordinates of the parks, we can find the order of park visits for our tour:

order = Last[
  FindShortestTour[
   GeoPosition /@ (Select[
      coordsallparks, (NumberQ[#[[1, 1]]] && 
         NumberQ[#[[1, 2]]] && #[[1, 1]] < 50) &])]]

Last but not least we can plot that:

GeoGraphics[GeoPath[parkpos[[order]][[All, 1]]]]

Hope this helps,

Marco

PS: Looking at the graphic I don't think that all of the parks qualify as city parks....

Using Marco's approach in a bit different way:

montParks = GeoEntities[Entity["AdministrativeDivision", {"Montana", "UnitedStates"}], "Park"]

order = Last[FindShortestTour[GeoPosition /@ montParks]]
Out[] = {1, 4, 5, 6, 3, 2, 1}

ImageMultiply[

  GeoListPlot[montParks, GeoZoomLevel -> 7, 
   GeoBackground -> "StreetMap", ImageSize -> 900],

  GeoListPlot[montParks[[order]], PlotMarkers -> GeoMarker, 
   GeoBackground -> "ContourMap", Joined -> True, ImageSize -> 900]

  ] // Sharpen