Find the GeoArea of overlapping GeoDisks

Question

Given a set of GeoDisks, how can one determine the combined region area, presumably through GeoArea? GeoDisks may overlap or be disjoint. Disjoint disks should sum to the total area, while overlapping disks should only contribute their union area.

With regular Disk objects, one can use RegionUnion on the set of Disks and then use Area, and this works fine. However, there doesn't seem to be a similar "GeoUnion" for geographical objects. Any ideas?

Answer 1

Depending on your accuracy goals, one possible way is to generate GeoGraphics, extract out the Polygon(s), use BoundaryDiscretizeGraphics to convert to a MeshRegion, extract the coordinates from the mesh, reverse order so it's {lat, lon}, make it back into a single Polygon, and finally use GeoArea.

geoAreaDisks[d__GeoDisk] :=
 Module[{polys},
  polys = 
   Cases[InputForm[
     GeoGraphics[{GeoStyling[Opacity[1]], d}, 
      GeoBackground -> None]], p_Polygon :> p, \[Infinity]];
  UnitConvert[
   GeoArea[Polygon@
     GeoPosition[
      MeshCoordinates[
        BoundaryDiscretizeGraphics[polys]] /. {a_, b_} :> {b, a}]], 
   "Acres"]]

d1 = GeoDisk[Entity["City", {"Topeka", "Kansas", "UnitedStates"}], 
   Quantity[100, "Miles"]];
d2 = GeoDisk[Entity["City", {"Manhattan", "Kansas", "UnitedStates"}], 
   Quantity[75, "Miles"]];

geoAreaDisks[d1, d2]
(* Quantity[2.68948*10^7, "Acres"] *)