Minimising number of dynamic pages to map scattered points using ArcGIS Desktop?

Question

From time to time I have to produce mapbook to show points of interest. First step to create pages, using regular mesh:

I don't like the solution because a) there are some pages with single points (e.g. page 25) sitting on the edge and b) too many pages.

First issue is easy to fix using code, - move rectangle of page extent to center of relevant points extent:

I still don't like it, it looks very crowded because number of pages remains the same. Remember, they all end up being actual A3 paper pages in multiple copies of report !.

So I've cooked a code that reduce number of pages. In this example from 45 to 34.

I am not sure if this the best result that can be achieved,

What is the best strategy (pseudo code, publication, Python library) to shuffle through points in order to minimise number of given size rectangles to capture all of the points? Surely, someone discovered it in game theory, military art or fishing industry

This is update to original question:

This shows real extent and page size required:

Closer zoom showing 10 out of 164 pages:

Sample Point Feature Class

Rectangle size can change as soon as it stays within the limits, i.e. smaller is fine.

Answer 1

This is not the answer, I just thought I post Python solution for those who interested:

# ---------------------------------------------------------------------------
# PAGE MAKER
# 
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy, traceback, os, sys
from arcpy import env

width=650
height=500

try:
    def showPyMessage():
            arcpy.AddMessage(str(time.ctime()) + " - " + message)
    mxd = arcpy.mapping.MapDocument("CURRENT")
    points = arcpy.mapping.ListLayers(mxd,"points")[0]
    pgons = arcpy.mapping.ListLayers(mxd,"pages")[0]

    g=arcpy.Geometry()
    geometryList=arcpy.CopyFeatures_management(points,g)
    geometryList=[p.firstPoint for p in geometryList]
    curT = arcpy.da.InsertCursor(pgons,"SHAPE@")
    while True:
        nPoints=len(geometryList)
        small=[geometryList.pop(0)]
        for p in geometryList:
            small.append(p)
            mPoint=arcpy.Multipoint(arcpy.Array(small))
            ext=mPoint.extent
            cHeight=ext.height
            cWidth=ext.width
            if cHeight>height or cWidth>width:
                small.remove(p)
        mPoint=arcpy.Multipoint(arcpy.Array(small))
        ext=mPoint.extent
        xC=(ext.XMin+ext.XMax)/2
        yC=(ext.YMin+ext.YMax)/2
        LL=arcpy.Point (xC-width/2,yC-height/2)
        UL=arcpy.Point (xC-width/2,yC+height/2)
        UR=arcpy.Point (xC+width/2,yC+height/2)
        LR=arcpy.Point (xC+width/2,yC-height/2)
        pgon=arcpy.Polygon(arcpy.Array([LL,UL,UR,LR]))
        curT.insertRow((pgon,))
        short=filter(lambda x: x not in small,geometryList)
        arcpy.AddMessage('Grabbed %i points, %i to go' %(len(small),len(short)))
        if len(short)==0: break
        geometryList=short[:]
    del mxd
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

applied it lately for survey planning:

UPDATE:

It seems that for some patterns dealing with ‘stray’ points first is the way to go. I’ve used ‘convex hull’ peel to identify them, idea of whuber, can’t find post, sorry.

Answer 2

This looks like the geometric version of the Maximum Coverage Problem which is closely related to the Set Cover Problem, and those two are NP-Complete.

So to solve it, one could use approximation. I would try the following algorithm and it seems to work perfectly. Although due to the complexity of problem, we cannot find the best answer.

1. Foreach point generate N=10 rectangles in random distances; just to make sure the rectangle covers the point (each rectangle has at least one point belongs to it and each point belongs to at least one rectangle)
2. Repeat until all points are covered: get rectangle covering max number of uncovered points. Mark points as covered.

an implementation of this algorithm, only for circle, is here: http://jsfiddle.net/nwvao72r/3/