Skip to main content

geography - How to test if a GeoPosition is inside a geographic Polygon?


I have tried the following:


usentity = Entity["Country", "UnitedStates"];
uspoly = usentity["Polygon"];
stdpoly = uspoly /. GeoPosition -> Identity;
RegionQ[stdpoly]

[1] True



RegionMember[stdpoly, {35, -90}]

returns unevaluated.


Is there a better approach?



Answer



If you look at the structure of the Polygon, it seems there is an extra level of {} that needs to be removed


Short@stdpoly
(* Polygon[{{{49.0021,-122.758},<<2620>>}}] *)

So you can use FlattenAt to remove that layer of brackets,



stdpoly2=FlattenAt[stdpoly,{1,1}];
Short@stdpoly
(* Polygon[{{{49.0021,-122.758},<<2620>>}}] *)

The polygons appear the same


RegionPlot /@ {stdpoly, stdpoly2}

Mathematica graphics


and now you can use RegionMember all you like


RegionMember[stdpoly2, {35, -90}]

(* True *)

There is also the fact that for some reason, geographic polygons have their coordinates reversed, as I learned here. This is what I did there,


stdpoly = Polygon @@ First@First@uspoly // Map[Reverse, #, {-2}] &;
Graphics@stdpoly
RegionMember[stdpoly, {-90, 35}]

Mathematica graphics


(* True  *)


As rcollyer points out, you can use the built-in GeoWithinQ to test this as well,


GeoWithinQ[Entity["Country", "UnitedStates"], GeoPosition[{35, -90}]]
(* True *)

Comments