19

I need to apply the four colors theorem in a polygonal shape in a way that I do not need to choose manually each color to put in each region. I wish to knows if there is any extension, plug-in, script or data base that may be used with ArcGIS and ArcToolBox to do it mathematically or programatically, so I could use it for now on with every map I come to create.

enter image description here enter image description here enter image description here

  • 1
    I also would like to know if there is this functionality in other systems other than ArcGIS, like QuantumGIS... – Please_Dont_Bully_Me_SO_Lords Mar 10 '13 at 17:45
  • 2
    I posted a suboptimal solution on GIS (with working R code) and an optimal solution (which will use three or even two colors if they can be found to work) on Mathematica. That solution is recursive; the reply to my post gives a linear programming solution. Manifold GIS has long had a five-color algorithm built in. (Four-coloring is hard to do; five-coloring is relatively simple to achieve.) – whuber Mar 11 '13 at 17:02
  • If you have no "code so far" my ArcGIS for Desktop recommendation would be to start with the Polygon Neighbors tool to get a table listing all neighbours of each polygon. – PolyGeo Jan 8 '14 at 9:25
  • @PolyGeo: thanks for the tools (I didn't know it) but I could not use it to solve my problem – radouxju Jan 8 '14 at 9:45
12

First of all, thanks for all answers and comments. Unfortunetaly, the existing tools were not fully compatible with the latest versions of QGIS and ArcGIS. Therefore I made my own solution using the tool indicated by @polygeo, the QGIS plugin from @Alexandre and the name of the algorithm (four color map) from @Jens.

Here is my code for those interested (for ArcGIS but the second part could be used in QGIS as well).

arcpy.MakeFeatureLayer_management(fc, fc[:-4]+ "_lyr" )
try:
    arcpy.AddField_management(fc[:-4] + "_lyr", "color", "SHORT")
except:
    print "field alread exists"   
arcpy.CalculateField_management(fc[:-4] + "_lyr", "color",  "10" , "PYTHON")

arcpy.PolygonNeighbors_analysis(fc[:-4] + "_lyr", fc[:-4] + "_tb.dbf" )
graph = []
cursor=arcpy.da.SearchCursor( fc[:-4] + "_tb.dbf" , ("src_FID","nbr_FID") )
for row in cursor:
    graph.append(row)


pols = arcpy.da.UpdateCursor(fc[:-4] + "_lyr", ("OID@","color"))
colored = []
for pol in pols:
    nbrs = [ second for first, second in graph if first == pol[0]]
    usedcolors = []
    for nbr in nbrs:
        usedcolors += [second for first, second in colored if first == nbr]
    pol[1]=[color for color in range(10) if color not in usedcolors][0]
    colored.append(pol)
    pols.updateRow(pol)

Note that the algorithm does not guarantee that only 4 colors are used: though it has been proven that the solution exists, the "brute force" is necessary to achieve it. In my case, I got 7 colors which is small enough. The script could have an additional loop until the solution is found, but I need to do it for hundreds of maps and 7 colors is OK.

  • 2
    This is brilliant - thanks so much for sharing it. I noticed at ArcGIS 10.2 the field names on the PolygonNeighbors output table have changed slightly - the fields are now called src_OBJECT and nbr_OBJECT – Stephen Lead May 15 '14 at 1:01
  • Is this script optimal i.e. does it ensure that a minimum of colors would be used? – Below the Radar Feb 8 '17 at 14:50
  • 1
    As far as I understood, raw force is needed. As mentioned in my post, you must run it several times to have a chance to reach the 4 colours. – radouxju Feb 8 '17 at 15:17
  • Still works great! Maybe the src_* and nbr_* field names depend on input type. I ran it now with a geodatabase fc input and Desktop 10.5 and they were named src_OBJECTID and nbr_OBJECTID. The script could be adjusted to list the fields that start with src and nbr so input type (or version of ArcGIS) does not matter. – BERA Jan 31 '18 at 10:35
4

There is a VB6 developer sample and an ArcGIS 9.x geoprocessing tool but from the comments on this ArcGIS Idea they don't work at 10.0+.

Perhaps someone would be interested in porting it.

A QGIS solution called TopoColour is given in the comments of this related question: Color polygons so each is distinct from its neighbors

3

If you are using QGIS, I believe that what you need is the Coloring a map plugin.

Unfortunatly, the plugin is only available for QGIS 1.8 version, but you can always download and see how the code works!

3

This is an adaptation of @radouxju's answer into a function. It will add a color field to the input feature layer and calculate. It should work regardless of the field name endings of PolygonNeighbors (they seem to be different for different users/inputs/arcgis versions (?))

def color_me(feature_layer):
    import arcpy
    try:
        arcpy.AddField_management(feature_layer, 'color', 'SHORT')
    except:
        print 'field alread exists'   

    arcpy.CalculateField_management(feature_layer, 'color',  '10' , 'PYTHON')

    arcpy.PolygonNeighbors_analysis(feature_layer, r'in_memory\neighbor_table' )
    graph = []
    neighbor_fields = [f.name for f in arcpy.ListFields(r'in_memory\neighbor_table') if f.name.startswith(('src', 'nbr'))]
    cursor=arcpy.da.SearchCursor(r'in_memory\neighbor_table' , neighbor_fields)
    for row in cursor:
        graph.append(row)

    pols = arcpy.da.UpdateCursor(feature_layer, ('OID@','color'))
    colored = []

    for pol in pols:
        nbrs = [ second for first, second in graph if first == pol[0]]
        usedcolors = []
        for nbr in nbrs:
            usedcolors += [second for first, second in colored if first == nbr]
        pol[1]=[color for color in range(10) if color not in usedcolors][0]
        colored.append(pol)
        pols.updateRow(pol)
    arcpy.Delete_management(r'in_memory\neighbor_table')

enter image description here

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.