# How to make geometric intersection faster for geographic borders?

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,
Module[{borderBox, borderPoly, boxIntersections, intersections,
positions, borderPoints},
borderBox =
Rectangle @@ (GeoBoundingBox[region,
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]},
boxToPoly,
RegionIntersection[
boxToPoly, (entities[[index]][
"Polygon"] /. {GeoPosition :> Identity})]
]
]
],
Range[Length@entities], {1}];
positions =
Position[intersections, x_ /; ! (Head@x === EmptyRegion), {1},
Extract[entities, positions]]


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

gulf = OceanData["GulfMexico"];
usStates =
EntityValue[
CountryData["UnitedStates"]}], "Entities"];

(* Get bordering US states *)

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

(* Get counties of bordering US states *)

EntityValue[
"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?

• You could try compiling some/all of your code. I've never actually done it so I can't say how much it would speed up computation, but I know compiled code should be faster than interpreted. This link might be helpful. mathematica.stackexchange.com/questions/1803/… Commented Sep 28, 2015 at 4:03
• @BrianG I'm having a lot of issues trying to compile this. The entity polygons can contain more than one polygon in some instances. I've extracted the points out Polygon with this rule //. {Polygon -> Identity, {{x__}} -> {x}} and set the compile parameter as {x, _Real, 2} but it still says that the values passed in are not of rank 2. Let us hope this catches the eye of someone more experienced with compile. Commented Sep 28, 2015 at 16:51

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).

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,
Module[{borderBox, borderPoly, boxIntersections, intersections,
positions, borderPoints},
borderBox =
Rectangle @@ (GeoBoundingBox[region,
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)

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[
"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.

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

• I believe I am doing this in the second ParallelMap in the function where I intersect the rectangle region with the border polygon and check for EmptyRegion before continuing with the entity polygon. Anything intersected with EmptyRegion immediately gives EmptyRegion so I'm not expecting a benefit of deconstructing the function like this. I'm currently running a test to see if wrapping the With in the second ParallelMap in an If that test Head[boxIntersections[[index]]] === EmptyRegion has any effect. Commented Sep 29, 2015 at 9:10
• No noticeable difference in execution time. :-( Commented Sep 29, 2015 at 10:11
• @Edmund OK I made some changes to completely avoid polygon region intersection. Commented Sep 29, 2015 at 21:21
• The updated solution does not work in cases were the polygons do not overlap. This is the case in most cases. I actually tried a point-wise solution first in the linked question. Your solution finds no bordering entities examples described in that link. It works for this particular gulf question but not for any other case mentioned (New York state to Pennsylvania counties, West Virginia state to Ohio counties, etc.). Commented Oct 1, 2015 at 9:20
• @Edmund Makes sense- so you need to increase the fidelity of the polygon. You can assume that all space between consecutive points is a line- so construct N midpoints for consecutive polygon coordinates and add to your test array. AND set the padding on the bounding boxes > 0 i.e. 1KM Commented Oct 1, 2015 at 16:14

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.

• When you say everything are you referring to all geo entities? Commented Sep 28, 2015 at 18:34
• Could you clarify what you mean in your answer? I don't think I fully understand your suggestion. Commented Sep 29, 2015 at 23:01
• Pre-process all the counties and bodies of water into an adjacency matrix, yes. Then you can easily run all sorts of pathfinding algorithms as well as easily seeing what borders what. Commented Sep 30, 2015 at 20:28
• math.uiuc.edu/~ash/Discrete/213Ch3.pdf has some interesting stuff about adjacency matrices. If you want to do stuff besides contiguity, you can also weight them for all sorts of fun pathfinding heuristics. Commented Sep 30, 2015 at 20:29
• In this site we value posting working code as answers. Commented Sep 30, 2015 at 20:51