GeoArea of an Intersection

Question

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

Answer 1

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.