For the sake of loading speed and to eliminate dependencies on outside files not allowed in Demonstrations, I have initialized this Demonstration with precomputed station locations and the distances between every pair of stations for each subway line.
I used borough outlines from [1] and subway station locations from [2]. I found three potentially tricky elements to the calculation — first, the city provides locations for subway entrances but not subway platforms, and the entrances are not necessarily centered around the platform; I chose a single entrance arbitrarily for each station and if I chose an outlying one, it will introduce a small error in the result. Second, the subway entrance locations are in {latitude, longitude} format (actually, {latitude, longitude} × 1,000,000), and this is the order expected by
Mathematica's
GeoDistance function, but the basemap is formatted in {longitude, latitude} format, so we have to reverse the order of the elements in each location when we switch between mapping and determining distances. Third, the city's database of subway entrances identifies an entrance as pertaining to a given line even if it is connected to that line only by an underground pedestrian tunnel. Since for the purposes of my distance calculations, I cared about what the trains are doing only, these extra stations had to be screened out before that computation was run.
Rather than checking only adjacent stations on each line for the closest stations, and only the ends of each line for the most distant points visited by a given train, I had
Mathematica check all possible permutations of stations on each line. As it turned out, the nearest stations were always adjacent, but the farthest ones were not necessarily line-ends. See the M train, for example.