# Using Geo commands, how do I get a list of cities that all meet the following criteria

I am looking to list with Geo commands all the cities such that

1. Point is 333 miles from Houston

2. Point is 755 miles from San Antonio

3. Point is 460 miles from El Paso.

How do I write up the code that gives a list of cities that meet three criteria?

• The wording of the question makes it sound to me like all 3 criteria must be met simultaneously, but I don't believe it's possible to satisfy all 3 constraints at once. Is it safe to assume you are just looking for a list of cities that satisfy #1, a separate list of cities that satisfy #2, etc.? Also, is there some leeway in the exact distance? Like would a city that is 340 miles from Houston still count? Jun 9 '20 at 21:54
• Yes, some small additional (or less) mileage is ok. I agree with you that there are no cities such that all three criteria are met. This is actually good. Either the point is on a road or such a city does not exist. Jun 9 '20 at 22:13

GeoWithinQ is reasonably fast:

threecities = Entity["City", {#, "Texas", "UnitedStates"}] & /@
{"Houston", "SanAntonio", "ElPaso"};

allcities = CityData[{All, "Texas", "UnitedStates"}];

Length @ allcities

1512

coords = Association[# -> CityData[#, "Coordinates"] & /@ allcities];

radii = Quantity[{333, 755, 460}, "Miles"];

selector = And @@@ Transpose[GeoWithinQ[disks, Values[coords]]];

citiesinintersection = Pick[allcities, selector]


GeoGraphics[{disks, PointSize[Small],
Point[GeoPosition[coords@#] & /@ allcities],
Red, PointSize[Large], Point[GeoPosition[coords@#] & /@ threecities],
Blue, PointSize[Medium], Point[GeoPosition[coords@#] & /@ citiesinintersection]}]


– kglr
Jun 26 '20 at 2:23

You can see from the following code that there is no point where the three circles intersect. If you want cities within those circles, though, the San Antonio requirement is redundant: as you can see, the intersection between the Houston and El Paso disks is the region of interest:

GeoGraphics[{
GeoDisk[Entity["City", {"Houston", "Texas", "UnitedStates"}], Quantity[333, "Miles"]],
GeoDisk[Entity["City", {"SanAntonio", "Texas", "UnitedStates"}], Quantity[755, "Miles"]],
GeoDisk[Entity["City", {"ElPaso", "Texas", "UnitedStates"}], Quantity[460, "Miles"]]}
]


One would hope to calculate the intersection of those regions, then lookup cities in that region. Unfortunately, however, there is no GeoRegionIntersection function (see Find the GeoArea of overlapping GeoDisks and GeoArea of an Intersection).

A second approach would be to get the list of all cities in those two relevant regions, then take the intersection of those lists:

houston = GeoNearest[
"City",
Entity["City", {"Houston", "Texas", "UnitedStates"}],
{All, Quantity[333, "Miles"]}
]


... and a similar one for El Paso (let's call it elpaso). Unfortunately again, though, that command times out for me; it does work for a smaller radius (e.g. 20 miles or so), so it must be a humongous list and my system / connection chokes on it.

If you are able to get the lists above, then Intersection[houston, elpaso] would do the trick.

• Is there a way to make it think of cities with, say, over 20,000 citizens or whatever? That criteria might help the system not choke. Maybe another intersection? Jun 9 '20 at 22:18