Geographic Information Systems Asked by FelixIP on January 16, 2021
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:
Rectangle size can change as soon as it stays within the limits, i.e. smaller is fine.
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.
Answered by FelixIP on January 16, 2021
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.
an implementation of this algorithm, only for circle, is here: http://jsfiddle.net/nwvao72r/3/
Answered by Farid Cheraghi on January 16, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP