Identify optimal subset of polygons based on cumulative overlap with a separate layer?

Question

Take the example below, where there are four shapes overlapped by a grid of points. Suppose that each shape has a "score" attribute linked to it, where score = number of sides of the shape (so circle = 1, triangle = 3, etc.). I want to find the optimal set (S) of k points (where k is a number decided beforehand), where S, when treated as a single unit, overlaps polygons with the highest possible summed score.

Provided k = 2 in the example below, the optimal set S would contain points 3 and 5, which collectively overlap shapes with a total score of 17 (hexagon [six] + rectangle [four] + square [four] + triangle [three]), higher than any other possible pair of points. Note that the answer is not points 2 and 3, because the score from the hexagon cannot be counted twice (so the total score of points 2 and 3 would be 13 [6+3+4] rather than 19 [6+3+6+4]). Thus, scores can't be calculated for each point individually using spatial join; they can only be calculated for specific pairs of points based on the unique shapes that the pair overlaps with.

Is there a way to have ArcMap determine the optimal set of points (S)? I also realize that this can be figured out manually for simple examples, but in my case there are thousands of points, which overlap range shapefiles for thousands of species, with each species having an associated score. My goal, specifically, is to identify the 20 points on the map that have the highest cumulative score-weighted species richness possible. I can't just choose the 20 points that have the highest scores individually, because they're all clustered next to each other on the map and thus harbor almost exactly the same species, which would limit the cumulative biodiversity captured by the top 20 points in total if using this method.

Answer 1

I hope I understand what you are after Find polygons intersecting each point by running:

GenerateNearTable("POINTS", "fish_net", ="C:/SCRATCH/SCRATCH.gdb/POINTS_neighbours", "0 Meters", closest="ALL", closest_count="0")

Transfer polygon weights into existing field "NEAR_RANK" in output table:

and run this script:

import arcpy,os, time
from collections import defaultdict
import itertools as itt

pgons = defaultdict(set)
ranks = dict()
points = set()
tbl= arcpy.da.TableToNumPyArray("neighbours",("IN_FID","NEAR_FID","NEAR_RANK"))
# populate dictionaries of intersections (pgons) and weight of indiv.polygon (ranks)
for point,pgon, rank in tbl:
    pgons[point].add(pgon)
    ranks[pgon]=rank
    points.add(point)
# iterate through combination by 2 and compute total score of unique polygons
list2sort=[]
for pA,pB in itt.combinations(points,2):
    comboSet = pgons[pA].union(pgons[pB])
    score = 0
    for item in comboSet:score+=ranks[item]
    list2sort.append((score,set([pA,pB])))
list2sort.sort(reverse = True)
first = list2sort.pop(0)
uniqueList = [first]
# weed out repeated points down sorted list
bigSet = first[1]
for item in list2sort:
    thisSet = item[1]
    unique = thisSet.difference(bigSet)
    if len (unique)<2: continue
    bigSet=bigSet.union(thisSet)
    uniqueList.append(item)
for item in uniqueList:      
    arcpy.AddMessage(item)

I ran it on:

These are output messages:

[17, set([3, 5])]
[15, set([8, 2])]
[10, set([1, 4])]
[6, set([9, 7])]

No wonder highest ranking polygon plays big role. With very little efforts you can force uniqueness of polygons as well, but this will result in just 1 pair being found, i.e. (3,5)