5
$\begingroup$

I'm trying to be able to calculate whether a coordinate is within a census tract. As an example, here's a particular census tract in Chicago defined as a set of polygon vertices:

coords = Uncompress[FromCharacterCode[
   Flatten[ImageData[Import["http://i.stack.imgur.com/fIe9F.png"],"Byte"]]]]

You can see the rough shape of the tract with Graphics[Polygon[GeoPosition[coords]]].

I was hoping one of the following would tell me whether a coordinate was in the tract:

RegionMember[GeoPosition[coords], {41.7284, -87.741}]
RegionMember[Polygon[GeoPosition[coords]], {41.7284, -87.741}]

But both complain of improperly specified regions. Any ideas?

FWIW: I hope to eventually be numerically integrating over the map. For each point in the numerical integration, I'll need to find the census tract to which it belongs. So this could get a bit dicey if I'm computationally inefficient.

$\endgroup$
9
$\begingroup$

You could use GeoWithinQ, as mentioned in JasonB's answer to How to test if a GeoPosition is inside a geographic Polygon?:

GeoWithinQ[GeoPosition[coords], GeoPosition[{41.7284, -87.741}]]

True

However, since you're interested in speed, you might consider RegionMember without GeoPosition:

GeoWithinQ[GeoPosition[coords],GeoPosition[{41.7284, -87.741}]] //AbsoluteTiming

rf = RegionMember[Polygon[coords]];
rf[{41.7284, -87.741}] //AbsoluteTiming

{0.001123, True}

{8.*10^-6, True}

$\endgroup$
  • $\begingroup$ This highlights how much can be gained by using the operator form of RegionMember. I had naively dismissed it because the 2-argument form seems to be so slow: RegionMember[Polygon@coords, {41.7284`, -87.741`}] // RepeatedTiming is slower than the GeoWithinQ so much that I was trying to roll my own winding-number-based compiled function. But if you build the operator first it's blazing fast. $\endgroup$ – Jason B. Jul 10 '17 at 17:13
  • $\begingroup$ For this particular application, making a region member function (variable rf above) for each census tract, then testing the points against each function will likely be the most efficient way to tackle this problem. $\endgroup$ – Peter Roberge Jul 10 '17 at 17:26
  • $\begingroup$ I hadn't thought of just using the coordinates themselves. A very elegant solution! $\endgroup$ – Shane Jul 10 '17 at 17:43
3
$\begingroup$

You could try with GeoWithinQ:

pol = Polygon[GeoPosition /@ coords];
GeoWithinQ[pol, GeoPosition[{41.7284`, -87.741`}]]

(* Out: True *)
$\endgroup$
  • $\begingroup$ Excellent. Thanks! I'll wait to see if there are other bright ideas with their own benefits but this definitely an acceptable answer. $\endgroup$ – Shane Jul 10 '17 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.