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 =
      Rectangle @@ (GeoBoundingBox[#, 
          geoPadding] /. {GeoPosition :> Identity})] &,
    entities, 1];
  (* Faster to split them *)
  intersections = ParallelMap[
     With[{boxToPoly = 
         RegionIntersection[boxIntersections[[index]], borderPoly]},
      If[Head[boxToPoly] === EmptyRegion,
        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 = 
   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 = 
     Entity["AdministrativeDivision", {EntityProperty[
         "AdministrativeDivision", "ParentRegion"] -> #}], 
     "Entities"] & /@ borderStates ;

(* Get counties bordering Gulf of Mexico *)

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

(* Show *)

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

enter image description here

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

  • $\begingroup$ 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/… $\endgroup$ – Brian G Sep 28 '15 at 4:03
  • $\begingroup$ @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. $\endgroup$ – Edmund Sep 28 '15 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).

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


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 with Count

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

  • $\begingroup$ 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. $\endgroup$ – Edmund Sep 29 '15 at 9:10
  • $\begingroup$ No noticeable difference in execution time. :-( $\endgroup$ – Edmund Sep 29 '15 at 10:11
  • $\begingroup$ @Edmund OK I made some changes to completely avoid polygon region intersection. $\endgroup$ – Peter Roberge Sep 29 '15 at 21:21
  • $\begingroup$ 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.). $\endgroup$ – Edmund Oct 1 '15 at 9:20
  • $\begingroup$ @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 $\endgroup$ – Peter Roberge Oct 1 '15 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.

  • 1
    $\begingroup$ When you say everything are you referring to all geo entities? $\endgroup$ – Edmund Sep 28 '15 at 18:34
  • $\begingroup$ Could you clarify what you mean in your answer? I don't think I fully understand your suggestion. $\endgroup$ – MarcoB Sep 29 '15 at 23:01
  • $\begingroup$ 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. $\endgroup$ – laudiacay Sep 30 '15 at 20:28
  • $\begingroup$ 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. $\endgroup$ – laudiacay Sep 30 '15 at 20:29
  • $\begingroup$ In this site we value posting working code as answers. $\endgroup$ – Dr. belisarius Sep 30 '15 at 20:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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