3
$\begingroup$

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

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
3
$\begingroup$

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

asd

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

union

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.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ A way to get the correct coordinates would be to add the option GeoProjection->"Equirectangular" in the GeoGraphics call inside geoIntersectionPolygon. $\endgroup$
    – jose
    Dec 6, 2017 at 3:06
  • 2
    $\begingroup$ Actually, instead of the FirstCase[...]& you can use the function GeoGraphics`GeoEvaluate. That is, use polys = GeoGraphics`GeoEvaluate /@ {p1, p2}; . There is some info on how this function works in the GeoGraphics tutorial in the documentation. $\endgroup$
    – jose
    Dec 6, 2017 at 3:18
  • 1
    $\begingroup$ @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. $\endgroup$
    – b3m2a1
    Dec 6, 2017 at 3:20
  • $\begingroup$ @jose, wow, I'm impressed with the compactness of the Geographics`GeoEvaluate 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 = GeoGraphics`GeoEvaluate /@ {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! $\endgroup$
    – Stephen Wilkus
    Dec 7, 2017 at 3:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.