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

Question

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?

Answer 1

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"];

disks = MapThread[GeoDisk, {coords /@ threecities, radii}];

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]}]

Answer 2

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.