How to make geometric intersection faster for geographic borders?

Question

I have developed a solution to identify geographic regions bordering a region in this post (94723).

Can this functionality be made to execute any faster? Find which entities in a set border other entities. For example, counties that border a neighbouring state (New York state to Pennsylvania counties, West Virginia state to Ohio counties, etc.), Gulf of Mexico to USA state counties, ect. The current execute takes over 40+ minutes for Gulf of Mexico to USA state counties. Can the run time be reduced by a factor of 10 or more.

geoBorderingEntities[region_, entities_List?VectorQ, 
  geoPadding_?QuantityQ] :=
 Module[{borderBox, borderPoly, boxIntersections, intersections, 
   positions, borderPoints},
  borderBox = 
   Rectangle @@ (GeoBoundingBox[region, 
       geoPadding] /. {GeoPosition :> Identity});
  borderPoly = region["Polygon"] /. {GeoPosition :> Identity};
  boxIntersections =
   ParallelMap[
    RegionIntersection[
      borderBox,
      Rectangle @@ (GeoBoundingBox[#, 
          geoPadding] /. {GeoPosition :> Identity})] &,
    entities, 1];
  (* Faster to split them *)
  intersections = ParallelMap[
    Function[{index},
     With[{boxToPoly = 
        Evaluate@
         RegionIntersection[boxIntersections[[index]], borderPoly]},
      If[Head[boxToPoly] === EmptyRegion,
       boxToPoly,
       RegionIntersection[
        boxToPoly, (entities[[index]][
           "Polygon"] /. {GeoPosition :> Identity})]
       ]
      ]
     ],
    Range[Length@entities], {1}];
  positions = 
   Position[intersections, x_ /; ! (Head@x === EmptyRegion), {1}, 
    Heads -> False];
  Extract[entities, positions]]

However, this takes too long (about 40+ minutes) to complete for larger problems. For example:

gulf = OceanData["GulfMexico"];
usStates = 
  EntityValue[
   Entity["AdministrativeDivision", {EntityProperty[
       "AdministrativeDivision", "ParentRegion"] -> 
      CountryData["UnitedStates"]}], "Entities"];

(* Get bordering US states *)

borderStates = geoBorderingEntities[gulf, usStates, Quantity[5, "Kilometers"]];

(* Get counties of bordering US states *)

borderStatesAdminDivs = 
  EntityValue[
     Entity["AdministrativeDivision", {EntityProperty[
         "AdministrativeDivision", "ParentRegion"] -> #}], 
     "Entities"] & /@ borderStates ;

(* Get counties bordering Gulf of Mexico *)

borderCounties = geoBorderingEntities[gulf, #, Quantity[5, "Kilometers"]] & /@ borderStatesAdminDivs ;

(* Show *)

GeoGraphics[{
   {EdgeForm[Black], White, Sequence @@ (Polygon[#] & /@ borderStatesAdminDivs)}, 
   {Red, gulf["Polygon"]}, 
   {EdgeForm[Gray], Blue, GeoStyling[Opacity[1]], Tooltip[Polygon[#], #["Name"]] & /@ Flatten[borderCounties]}}, 
 GeoBackground -> None]

How can I speed up the region intersections to find polygons in entities that are touching the region polygon?

Answer 1

An order of 10! I can't get you there, but I think I can add some average RT improvement.

You have some up-front time costs like fetching values from the research servers- we'll ignore these, as it can't be improved without caching / downloading.

Where I see enormous gains are using your simple operations (Fast) before you delve into the expensive operations (slow).

Your simple operations are your rectangular region intersections, whereas your expensive operations are your polygon region intersections.

If you were to process all US states in this example as rectangles, THEN process all the Administrative Divisions as rectangles you eliminate a bunch of your polygon region intersections:

We first chop off the second half of your function:

geoBorderingEntities[region_, entities_List?VectorQ, 
  geoPadding_?QuantityQ] := 
 Module[{borderBox, borderPoly, boxIntersections, intersections, 
   positions, borderPoints}, 
  borderBox = 
   Rectangle @@ (GeoBoundingBox[region, 
       geoPadding] /. {GeoPosition :> Identity});
  borderPoly = region["Polygon"] /. {GeoPosition :> Identity};
  boxIntersections = 
   ParallelMap[{#, 
      RegionIntersection[borderBox, 
       Rectangle @@ (GeoBoundingBox[#, 
           geoPadding] /. {GeoPosition :> Identity})]} &, entities, 1]

(I also changed your output just to make fiddling with this easier for me)

Now, assuming your previous initialization,

borderStates = geoBorderingEntities[gulf, usStates, Quantity[5, "Kilometers"]];
xStates = Select[borderStates, Head[#[[2]]] =!=  EmptyRegion &][[All, 1]]

xStates will now contain your states, and didn't do any expensive calculations.

Now, we get the sub-regions:

counties = 
  Flatten[EntityValue[
      Entity["AdministrativeDivision", {EntityProperty[
          "AdministrativeDivision", "ParentRegion"] -> #}], 
      "Entities"] & /@ xStates , 1] ;

Leaving us with 693 areas to check.

borderCount = geoBorderingEntities[gulf, counties, Quantity[5, "Kilometers"]];
borderCount = Select[borderCount, Head[#[[2]]] =!= EmptyRegion &][[All, 1]]

Leaving us with 180 areas to do expensive calculations on!

See EDIT from here on

From here you can call your original code to find the actual overlaps.

Now, my way wont always work, but if your shape is pretty box-shaped, it will speed it up. If you have an entity like a river with a diagonal shape, you will see basically no performance gains from the above.

EDIT

I decided to do away with the expensive calculations- we're only going to work with coordinates, because we don;t really need overlapping regions, we just need to know if a county is next to the region.

So I will think about this instead as is the Gulf inside the county

Obtain coordinates for the target Polygon, Obtain all bounding boxes from the original reduced set of counties (borderCount):

gulfCord = gulf["Polygon"] /. {Polygon -> Identity, GeoPosition -> Identity} ;
allBB = (GeoBoundingBox[#] /. {Polygon -> Identity, GeoPosition -> Identity} & /@ borderCount);

Make the bounding boxes into rectangles:

allRec = Apply[Rectangle, allBB, 1];

Now test all the gulf coordinates to see if they are inside of the counties:

pointTest = Boole[RegionMember[#, gulfCord[[1]]] & /@ allRec];

Grab only the indices of the regions detected as true:

indexOfIn = Position[pointTest, 1][[All, 1]] // DeleteDuplicates;
reducedSet = Flatten[borderCount[[#]] & /@ indexOfIn, 1];

Plot the results:

GeoGraphics[{EdgeForm[Gray], Blue, GeoStyling[Opacity[1]], 
  Tooltip[Polygon[#], #["Name"]] & /@ Flatten[reducedSet], Red, 
  gulf["Polygon"]}]

This gets you 10x faster

Edit 2

The above will not work when the polygon contains consecutive distant points.

So, the solution is to create more points across the distant points

gc = gulfCord[[1]];
midpoint[a_, b_] := ({(a[[1]] + b[[1]])/2, (a[[2]] + b[[2]])/2});
createPoints[x_] := 
  Flatten[Table[{x[[i]], midpoint[x[[i]], x[[i + 1]]]}, {i, 1, 
     Length@x - 1}], 1];
expandedPointSet = Nest[createPoints[#] &, gc, 10];

Now you can use this expanded point set to test if these points are in the counties

Note: the larger the factor in the Nest function, the more point samples

Answer 2

Pre-process everything into an enormous adjacency matrix data structure, or possibly a graph. It would take a very long time now, but be blazing fast later.