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}
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}]
(* 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
Post a Comment