I am looking to list with Geo commands all the cities such that
Point is 333 miles from Houston
Point is 755 miles from San Antonio
Point is 460 miles from El Paso.
How do I write up the code that gives a list of cities that meet three criteria?
Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this communityI am looking to list with Geo commands all the cities such that
Point is 333 miles from Houston
Point is 755 miles from San Antonio
Point is 460 miles from El Paso.
How do I write up the code that gives a list of cities that meet three criteria?
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]}]
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.