4

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?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Alan
  • 43
  • 3

1 Answers1

4

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"] *)
chuy
  • 11,205
  • 28
  • 48
  • Nice work with the clever approach. It's too bad Mathematica can't do this directly with geographic regions like it does with geometric ones. Thank you. – Alan May 23 '17 at 00:33