# GeoArea of an Intersection

I am looking to calculate the geographic area of an intersection of a GeoDisk with a geographic entity such as a county.

I was hoping for a simple construct like:

state = Entity["Country", "PuertoRico"];
disk = GeoDisk[{18, -67}, Quantity[50, "Kilometers"]];
Print[ GeoGraphics[{disk, state}, GeoScaleBar -> True] ];
Print[" The overlap area is: ",
GeoArea[RegionIntersection[ Polygon[ state ], disk] ] ]


Unfortunately, the RegionIntersection function doesn't handle GeoGraphics structures, and I can find no such thing as a "GeoRegionIntersection" native function. Is there such a function hidden, but in trial, perhaps?

Alternatively, is there a reasonably straightforward way to convert the GeoDisk to a Polygon and do the calculation with something like:

Area[RegionIntersection[ Polygon[state], function[disk] ] ]?


I hope to apply this to several hundred such disks and am also interested in doing this for about a hundred counties for which I will take a weighted sum of population per county.

Thank you in advance, Steve

• Dec 5 '17 at 4:39

## 1 Answer

Here's an adaption of chuy's function from Jason B's link in the comments:

geoIntersectionPolygon[p1_, p2_] :=
Module[
{polys, regs, coords},
polys =
FirstCase[
InputForm[
GeoGraphics[{GeoStyling[Opacity[1]], #},
GeoBackground -> None]], p_Polygon :> p,
Point[{0, 0}],
\[Infinity]
] & /@ {p1, p2};
regs =
BoundaryDiscretizeGraphics /@ polys;
regs = Apply[RegionIntersection, regs];
If[Head[regs] === EmptyRegion,
None,
coords = MeshCoordinates[regs] /. {a_, b_} :> {b, a};
Polygon@
GeoPosition[coords[[FindShortestTour[coords][[2]]]]]
]
]


Testing:

geoIntersectionPolygon[
Entity["Country", "PuertoRico"]["Polygon"],
disk
] // GeoGraphics


Note that the shape is right, but the coordinates are off.

We can define a companion union function:

geoUnionPolygon[p1_, p2_] :=
Module[
{polys, coords},
polys =
Cases[
InputForm[
GeoGraphics[{GeoStyling[Opacity[1]], p1, p2},
GeoBackground -> None]], p_Polygon :> p,
\[Infinity]
];
coords =
MeshCoordinates[
BoundaryDiscretizeGraphics[polys]] /. {a_, b_} :> {b, a};
Polygon@
GeoPosition[coords[[FindShortestTour[coords][[2]]]]]
];

geoUnionPolygon[
Entity["Country", "PuertoRico"]["Polygon"],
disk
] // GeoGraphics


And we see that it seems to be primarily an issue of scaling. How to appropriately handle that is left to the reader. We do, however, recover the general shape, which should be enough for computing the area up to a constant scaling factor.

One further issue here is that the way this works doesn't handle disjoint regions appropriately. E.g. the islands off of the East coast of Puerto Rico shouldn't be connected to the region like that.

That's a kink I don't want to deal with though.

• A way to get the correct coordinates would be to add the option GeoProjection->"Equirectangular" in the GeoGraphics call inside geoIntersectionPolygon.
– jose
Dec 6 '17 at 3:06
• Actually, instead of the FirstCase[...]& you can use the function GeoGraphicsGeoEvaluate. That is, use polys = GeoGraphicsGeoEvaluate /@ {p1, p2}; . There is some info on how this function works in the GeoGraphics tutorial in the documentation.
– jose
Dec 6 '17 at 3:18
• @jose I think you may know more about this system than do I. Want to provide an answer instead? I would rather have your good answer in place than my bad one. Dec 6 '17 at 3:20
• @jose, wow, I'm impressed with the compactness of the GeographicsGeoEvaluate function, but I've never seen the apostrophe used that way before. I read the documentation and have the near right answer shown by @b3m2a1 with polys = GeoGraphicsGeoEvaluate /@ {Polygon[ state ], disk}, but no luck with: polys = GeoGraphics`GeoEvaluate[#, GeoProjection -> "Equirectangular"] /@ {Polygon[ state ], disk}. I will keep exploring now that you've both been very helpful. Thank you! Dec 7 '17 at 3:36