4

I am attempting to assign one of 16 colors to map objects (e.g. counties, states, ZIP Codes, etc.) so that no two adjacent objects share the same color.

I'm seeing tons of articles online about "edge coloring" and the "four color theorem," but I cannot seem to find a way to apply these theorems as an algorithm in PostgreSQL/PostGIS (or any practical application, for that matter).

I'm 100% sure this is a solved problem, but since Google isn't revealing an answer, I'm guessing I'm too ignorant on this topic to ask the right question.

Can someone point me in the right direction?

4

I've had good luck using the 6-color algorithm described in the introduction to Two Linear-Time Algorithms for Five-Coloring a Planar Graph. Although it's certainly possible to color a graph using fewer colors, it may not look any better than using 5 or 6.

ALGORITHM 6-COLOR.

Given an n-vertex planar graph G in adjacency list form, this algorithm determines a 6-coloring of G.

Step 1. Establish degree lists.

For each j where 0 < j < n-1, form a doubly-linked list of all vertices of G of degree j.

Step 2. Label vertices smallest degree last.

For i=n, n-1, n-2, ... , 1, designate the first vertex of the non-vacuous j-degree list of smallest j as vertex vi.

Delete vi from the j degree list.

For each vertex v’ that was adjacent to vi in G and remains in some degree list, say j’, delete v’ from the j’ degree list and insert v’ in the j’ - 1 degree list.

Step 3. Color vertices.

For i =1,2,...,n, assign vertex vi the smallest color value (which must be some integer between one and six) not occurring on the vertices adjacent to vi that have already been colored.

|improve this answer|||||
3

I'm sure this function is incredibly inefficient (not to mention that it's coded very poorly -- it was written as a "once-off" quick and dirty implementation), but here's one whose results work quite well for me:

CREATE OR REPLACE FUNCTION public.calc_colors(
    tbl text,
    unique_field text,
    neighbour_style text,
    search_distance real DEFAULT 0)
  RETURNS SETOF record AS
$BODY$
DECLARE
 first_vertex text;
 r RECORD;
 next_color INTEGER;
 max_used_color INTEGER;
 available_colors INTEGER;
 next_vertex text;
 create_contig TEXT;
 color_string TEXT;
BEGIN

    --find neighbours
    IF neighbour_style = 'no_corners' THEN
        create_contig = 'CREATE TEMP TABLE contig_temp ON COMMIT DROP AS SELECT p1.' || quote_ident(unique_field) || '::text AS c1, p1.geom as geom, p2.' ||  quote_ident(unique_field) || '::text AS c2 FROM ' || tbl || ' p1 CROSS JOIN ' || tbl || ' p2 WHERE ( ST_Relate( p1.geom, p2.geom, ''2********'') OR ST_Relate( p1.geom, p2.geom, ''****1****'')) AND p1.' || quote_ident(unique_field) || ' != p2.' ||  quote_ident(unique_field);
    ELSIF neighbour_style = 'within' THEN
        create_contig = 'CREATE TEMP TABLE contig_temp ON COMMIT DROP AS SELECT p1.' || quote_ident(unique_field) || '::text AS c1, p1.geom as geom, p2.' ||  quote_ident(unique_field) || '::text AS c2 FROM ' || tbl || ' p1 CROSS JOIN ' || tbl || ' p2 WHERE (ST_Intersects( p1.geom, p2.geom) Or ST_DWithin(  st_envelope(p1.geom), st_envelope(p2.geom), ' || search_distance || ')) AND p1.' || quote_ident(unique_field) || ' != p2.' ||  quote_ident(unique_field);
    ELSE
        create_contig = 'CREATE TEMP TABLE contig_temp ON COMMIT DROP AS SELECT p1.' || quote_ident(unique_field) || '::text AS c1, p1.geom as geom, p2.' ||  quote_ident(unique_field) || '::text AS c2 FROM ' || tbl || ' p1 CROSS JOIN ' || tbl || ' p2 WHERE ST_Intersects( p1.geom, p2.geom) AND p1.' || quote_ident(unique_field) || ' != p2.' ||  quote_ident(unique_field);
    END IF;

    EXECUTE create_contig;

    CREATE INDEX sidx_contig_temp
          ON contig_temp
          USING gist
          (geom);


      CREATE TEMP TABLE vertex_degree_temp ON COMMIT DROP AS SELECT c1, count(*) as neighbour_count FROM contig_temp GROUP BY c1 ORDER BY count(*) DESC, c1;
      CREATE INDEX dv_tmp_idx ON vertex_degree_temp (c1);

      --color first vertex, which is vertex with highest number of neighbours
      SELECT c1 INTO first_vertex FROM vertex_degree_temp LIMIT 1;

      CREATE TEMP TABLE vertex_colors_temp ON COMMIT DROP AS SELECT first_vertex c1, 1::int color_idx;
      CREATE INDEX vc_tmp_idx ON vertex_degree_temp (c1);
      max_used_color = 1;

    LOOP

    --find next vertex, which is vertex with largest number of distinct colours in neighbours
    --for ties, choose vertex with highest number of neighbours
      SELECT d.c1 INTO next_vertex FROM (
    SELECT DISTINCT d.c1, d.neighbour_count, c.color_idx FROM vertex_degree_temp d, contig_temp n 
        LEFT JOIN vertex_colors_temp c ON (c.c1 = n.c2)
        WHERE d.c1 NOT IN (SELECT c1 FROM vertex_colors_temp) AND d.c1 = n.c1 ) d
        GROUP BY d.c1, d.neighbour_count
        ORDER BY count(d.color_idx) DESC, neighbour_count DESC
        LIMIT 1;

    EXIT WHEN next_vertex IS NULL;

    --find least used available color for vertex
    DROP TABLE IF EXISTS used_colors_temp;
    CREATE TEMP TABLE used_colors_temp AS SELECT * FROM (SELECT generate_series(1,max_used_color) a ) seq WHERE a NOT IN 
        (SELECT c.color_idx FROM contig_temp n, vertex_colors_temp c
        WHERE (c.c1 = n.c2) AND n.c1 = next_vertex ORDER BY color_idx);

    SELECT count(*) FROM used_colors_temp INTO available_colors;

    IF available_colors = 0 THEN
        next_color = max_used_color + 1;
        max_used_color = next_color;
    ELSE
        --choose colour with smallest count, but furthest distance from current geometry
        SELECT d.color_idx FROM
        (SELECT avail.color_idx, ST_Distance(st_envelope(c.geom), st_envelope(avail.geom)) dist FROM contig_temp c, 
        (SELECT u.c1, u.color_idx, c.geom FROM vertex_colors_temp u LEFT JOIN contig_temp c ON u.c1 = c.c1 WHERE color_idx IN (SELECT a FROM used_colors_temp) ) avail
        WHERE c.c1 = next_vertex) d GROUP BY d.color_idx ORDER BY min(dist) DESC LIMIT 1 INTO next_color;

    END IF;

    --insert
      INSERT INTO vertex_colors_temp (c1, color_idx) VALUES (next_vertex, next_color);

      max_used_color = GREATEST(max_used_color, next_color);

    END LOOP;

      RETURN QUERY SELECT * FROM vertex_colors_temp;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 1000;

It's used by joining it to the existing table, eg:

SELECT l.id,
    colors.color_idx,
    l.geom
   FROM public.localities l
     LEFT JOIN public.calc_colors('public.localities'::text, 'id'::text, 'within'::text, 2000::real) colors(id text, color_idx integer) ON colors.id::bigint = l.id

(Again, apologies for the horrid syntax... it was a throw-away implementation!)

This function will return a color number for each feature in the table specified as the first argument passed to the function. The second argument is a column which uniquely identifies each feature (used for joining the colors back to the original table). The magic is in the 3rd and 4th arguments - these control how the colors are assigned. The 3rd argument "neighbour_style" can be either:

  • blank, meaning that colours will be assigned in a way so that no intersecting features share the same colour number
  • or 'no_corners', meaning that no intersecting or touching features will share the same colour
  • or 'within', in which case the 4th parameter is used to specify a distance tolerance. No features within this distance of each other will share the same colour (can be useful if features almost touch each other and you don't want these to share the same colour).

The algorithm will attempt to assign colours so that the overall number of features with a given colour number is roughly equal, and so that the colours are roughly evenly spread over the geographic bounds of the input table.

Any suggestions for improvements to this script would be greatly appreciated!

|improve this answer|||||
2

Answering my own question here in the hope that this helps others (or my future self). In the intervening weeks since I posted the question, I was unable to find an elegant solution to this problem that could be done algorithmically in PostgreSQL. Instead, I broke the task down into constituent parts and more or less brute-forced it.

For the purposes of this walk-through, I'll use ZIP Codes as the example, but this approach worked just as well for counties, cities, states, and countries. And rather than posting very system-specific code that wouldn't be reusable, I'll just explain my approach.

1) I created an interim relational table that contained the IDs of each neighboring ZIP. (e.g. 12345 is adjacent to 12346, 12347, and 12348). To determine adjacency, I used ST_Intersects(). (See notes below about cities). To speed things up, I could have pre-filtered the query to limit the search for adjacent ZIPs within a certain distance, but I found I could process everything within a tolerable coffee break for this one-time run.

2) I then created a PHP script that processed every ZIP by a) picking a random color from a pool, and b) looking to see if any of the adjacent ZIPs via #1 above have already used that color. If the color had been used, then I looped through random colors from the pool until I found a fresh, unused color.

Notes:

  • For my color values, I simply used an int range from 0-15 in the first batch and 0-7 in a refined second attempt. This makes it much easier to pick a random value (e.g. rand(0,7)), and it provides more flexibility in styling down the line. In my case, I use that int value with conditional rules in CartoCSS to style the actual color values on the fly.
  • I found only one ZIP Code and three counties that had more than 7 neighbors (thus exhausting the available color pool in my 8 color set (these were extraordinarily tall shapes with lots of neighbors along each long side if you're wondering how this is possible). To prevent the script in #2 from looping infinitely, I simply slammed the last random value picked and lived with the duplicate color. This problem did not occur in the 16-color pool, obviously. I have read that the 50 US States can be done with as few as 4 colors, though I didn't try.
  • Cities are a little trickier since most don't butt right up against another. Thus, ST_Intersects() will not work. Instead, I defined adjacency as being within a certain short distance (too long, and you'll get too many neighbors). While this approach will not create a mathematically pure solution, it works in practice well enough to trick the eye when drawing the map.
  • The run-times to process the ZIPs and cities, respectively, with the script in #2 above were around an hour or so on my rig. Make sure your max_execution_time value in your php.ini is set to allow a really long job.

I hope this helps you. Feel free to ping me with any questions or clarifications.

|improve this answer|||||

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.