Seeking tool/algorithm for assigning code to Enumeration Areas (polygons) using serpentine numbering method?

One of the projects I am working on has asked me to assign a code (starting from 1 and counting up) to polygons that represent enumeration areas so the statistician can sample the areas. I have been asked that the assigning of the code uses a serpentine numbering method. The image below shows an example of what I am trying to achieve. You start in the North West part of the area and snake bank and forth assigning numbers.

Now I've searched various sites to see if anyone had developed a tool or even an algorithm that could be implemented that would achieve this. All I have found are numerous powerpoint slides and pdf files that simply say they have used serpentine numbering. I can only assume they have done this manually.

Does anyone know of an existing tool that would take a shapefile and number as described?

It seems to me that someone may have attempted this before but I have not been successful in finding a suitable reference as I don't usually get involved in census type data processing.

Answer

Assume 5 rows with five gridded polygons for each row. Further assume each row of height 1 and each polygon of width 1, with centroids at the physical centers, sequentially named "A" through "Y" in random order:

The query to grid the cells in this table (ArcGIS 10.2 geodatabase in PostgreSQL with SDE.ST_GEOMETRY) would be:

SELECT floor(sde.ST_MaxY(shape)) as ycell,

       floor(sde.ST_MinX(shape) + sde.ST_MaxX(shape) / 2.0) as xcell,
       objectid
FROM   polytab

And the code to label top-to-bottom, left to right would be:

UPDATE polytab SET seq_order = vtab.rownum
FROM 
(   SELECT row_number()
        over (order by ycell desc, xcell) as rownum, fkey
    FROM 

    (   SELECT floor(sde.ST_MaxY(shape)) as ycell,
               floor(sde.ST_MinX(shape) + sde.ST_MaxX(shape) / 2.0) as xcell,
               objectid as fkey
        FROM   polytab
    ) AS vpoly
) AS vtab
WHERE polytab.objectid = vtab.fkey;

Which would result in:

So far, this is the same solution as in the previous question, but now we break new ground: In order to "reverse the polarity" on even scanlines, a modulus function can be used to alter the "xcell" values, based on whether "ycell modulus 2" is one or zero:

UPDATE polytab SET serp_order = vtab.rownum
FROM 
(   SELECT row_number()
        over (order by ycell desc, xcell) as rownum, fkey
    FROM 
    (   SELECT floor(sde.ST_MaxY(shape)) as ycell,
                (floor(sde.ST_MinX(shape) + sde.ST_MaxX(shape) / 2.0) *
            (cast(floor(sde.ST_MaxY(shape)) as integer) % 2) +
                (5 - floor(sde.ST_MinX(shape) + sde.ST_MaxX(shape) / 2.0)) *

            (cast((1 + floor(sde.ST_MaxY(shape))) as integer) % 2)) as xcell,
                objectid as fkey
        FROM   polytab
    ) AS vpoly
) AS vtab
WHERE polytab.objectid = vtab.fkey;

(As you can see, modulus operations must be integer, so the terms start to get unwieldy as a CAST is necessary; also note that the width is hard-coded as "5").

Now the labels are:

Porting this to an application which works with a generic file geodatabase is left as an exercise.