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["https://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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
10
$\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}

Improve this answer
$\endgroup$
3
  • $\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, 2017 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, 2017 at 17:26
  • $\begingroup$ I hadn't thought of just using the coordinates themselves. A very elegant solution! $\endgroup$
    – Shane
    Jul 10, 2017 at 17:43
3
$\begingroup$

You could try with GeoWithinQ:

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

(* Out: True *)
Improve this answer
$\endgroup$
1
  • $\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, 2017 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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