Numbering polygons according to their spatial relationships?

Question

I have, for example, this layer in postgresql (polygons):

table : poly (objectid: integer, num: integer, shape: geometry)

I want to enumerate the polygons from left to right and top to bottom like this:

I am working with Vb.net, ArcEngine, PostgreSQL, and ArcSDE.

Solution that I proposed

I think to create grid lines (red layer; see picture) and iterate it, then verify intersection and finally enumerate the polygons. Will this work?

Answer 1

If you create an ancillary table with the attributes (fkey,xcell,ycell), populated from the objectid and a gridding of the maximum Y coordinate (or centroid coordinate or upper-center envelope coordinate), then you could use SQL to update the table, with something like:

    UPDATE poly SET num = vtab.rownum
    FROM 
    (   SELECT row_number() over (order by ycell desc, xcell) as rownum, fkey
        FROM ancillary_tab
    ) AS vtab
    WHERE poly.objectid = vtab.fkey;

(The UPDATE syntax is different from RDBMS to RDBMS, and PostgreSQL doesn't have the ROWNUM reserved word of ORACLE, but the moral equivalent would work in other databases.)

You didn't specify the cell size, but if it's just 1 degree, then a floor() math function would work. The more generic formula would be:

     cell = floor(coord / cellsize)

[Updated for gridding algorithm]

UPDATE: Thinking about it, you don't even need an ancillary table. I had to stand up a new database for this to prototype it, but after loading a similar dataset in PG 9.1.9 with ArcSDE 10.1sp1, this query:

UPDATE poly SET num = 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   poly
    ) AS vpoly
) AS vtab
WHERE poly.objectid = vtab.fkey;

Turned this:

Into this:

Answer 2

Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite, it's quite simple.

Suppose we want to order the labels from left to right. Firstly, we have to calculate the distance of each feature from a common reference on the left of all geometries (e.g. -180th meridian) and then order the geometries by distance.

ogr2ogr FRA_adm1_temp.shp FRA_adm1.shp -dialect sqlite -sql "SELECT *, ST_Distance(geometry, GeomFromText('LINESTRING(-180 90, -180 -90)')) AS DIST FROM FRA_adm1 ORDER BY DIST"

Finally, simply update the table adding a LABEL column with autoincremental numbers:

ogr2ogr FRA_adm1_label.shp FRA_adm1_temp.shp -sql "SELECT FID+1 AS LABEL, * FROM FRA_adm1_temp"

This approach can be used also in PostgreSQL.

Answer 3

Since it's unlikely that any two polygons have exactly the same northmost extent, you will need some compromise to order by both north-to-south and east-to-west.

Personally, I'd use either bounding boxes (looking at first maxY, then minX) or centroids (looking at first y, then x). I'd decide how much tolerance I'd allow and still consider y values to be approximately equal. Then I would collect sets from northmost to southmost (descending y ranges), arrange each set by ascending x values, and number across the sets. Ideal tolerance would probably need to be determined by trial and error.

This is very like a grid approach but probably simpler.

Note that I can't see latitude lines in your drawing but I assume not curved since your grid is not curved. I think both my approach and your grid approach would put 7 in the same set with 2-6, and 7 would be numbered 2.

A centroid approach might group 2,3, and 7 in one group, 4-6 and 9 in another. Bounding box might put 2-5 and 7 in one, 12, 9 and 6 in another.