CGI: Export geometry area and length/perimeter in the GeoJSON properties.

This causes overhead on some geometries though.  Querying 100000 random
features in each table and measuring yields the following.

Without geom_area/geom_perimeter/geom_length
============================================

nvr:SCI_Rikstackande (forced to 2D)
Time: min=0.000s, max=0.668s, avg=0.003s

ren:riks_ren
Time: min=0.000s, max=0.090s, avg=0.012s

sametinget:flyttled
Time: min=0.000s, max=0.003s, avg=0.000s

sks:UtfordAvverk
Time: min=0.000s, max=0.180s, avg=0.001s

With geom_area/geom_perimeter/geom_length
=========================================

nvr:SCI_Rikstackande (forced to 2D)
Time: min=0.000s, max=1.242s, avg=0.004s

ren:riks_ren
Time: min=0.000s, max=0.092s, avg=0.012s

sametinget:flyttled
Time: min=0.000s, max=0.016s, avg=0.000s

sks:UtfordAvverk
Time: min=0.000s, max=0.190s, avg=0.001s

1 files changed, 40 insertions, 7 deletions

diff --git a/webmap-cgi b/webmap-cgi
index 019f0bf..c4c4142 100755
--- a/webmap-cgi
+++ b/webmap-cgi
@@ -60,7 +60,7 @@ def get_query_map(layernames : set[str]) -> dict[str,bytes]:
     # pylint: disable-next=no-member
     with PG_CONN.cursor(binary=True, scrollable=False, withhold=False) as cur:
         for layername in layernames:
-            cur.execute('SELECT f_geometry_column, coord_dimension '
+            cur.execute('SELECT f_geometry_column, coord_dimension, srid, type '
                           'FROM ' + common.escape_identifier(SCHEMA_NAME) + '.geometry_columns '
                           'WHERE f_table_schema = %s AND f_table_name = %s',
                         params=(SCHEMA_NAME, layername),
@@ -69,13 +69,35 @@ def get_query_map(layernames : set[str]) -> dict[str,bytes]:
             if resp is None:
                 continue
             geom_cols = [ resp[0] ]
-            force2d = resp[1] > 2
+            geom_dim = resp[1]
+            force2d = geom_dim > 2
+            geom_type = resp[3]
+            geom_srid = resp[2]
+            if geom_srid != 3006:
+                # If the SRS isn't projected and/or isn't in meter units then ST_Area() resp.
+                # ST_Length() aren't in m² resp. m.  We could reproject, but since the target
+                # SRS is SWEREF 99 we just warn on mismatch
+                logging.warning('Geometry column "%s" in table "%s" has SRID %d != 3006',
+                                geom_cols[0], layername, geom_srid)
+            if geom_type in ('POLYGON', 'MULTIPOLYGON',
+                             'POLYGONM', 'MULTIPOLYGONM'):
+                d = 2 # surface
+            elif geom_type in ('LINESTRING', 'MULTILINESTRING',
+                               'LINESTRINGM', 'MULTILINESTRINGM'):
+                d = 1 # curve
+            elif geom_type in ('POINT', 'MULTIPOINT',
+                               'POINTM', 'MULTIPOINTM'):
+                d = 0 # point
+            else:
+                logging.warning('Geometry column "%s" in table "%s" has unknown type %s',
+                                geom_cols[0], layername, geom_type)
+                d = -1
 
             resp = cur.fetchone()
             if resp is not None:
                 logging.warning('Table "%s" has multiple geometry colums, '
-                                'only considering "%s"',
-                                layername, geom_cols[0])
+                                'only considering "%s" (%s, SRID=%d, dim=%d)',
+                                layername, geom_cols[0], geom_type, geom_srid, geom_dim)
                 while resp is not None:
                     geom_cols.append( resp[0] )
                     resp = cur.fetchone()
@@ -109,9 +131,14 @@ def get_query_map(layernames : set[str]) -> dict[str,bytes]:
                           'WHERE table_schema = %s AND table_name = %s',
                         params=(SCHEMA_NAME, layername),
                         prepare=False)
+            # never empty since the we know the table exists and has a primary key
             resp = cur.fetchone()
             while resp is not None:
                 c = resp[0]
+                if (c in ('layer_group', 'layer') or
+                        (d == 2 and c in ('geom_area', 'geom_perimeter')) or
+                        (d == 1 and c == 'geom_length')):
+                    logging.warning('Duplicate column name "%s"', c)
                 if c != pkey_col and c not in geom_cols:
                     column_names.append(c)
                 resp = cur.fetchone()
@@ -122,10 +149,16 @@ def get_query_map(layernames : set[str]) -> dict[str,bytes]:
             for column_name in column_names:
                 query +=     'm.' + common.escape_identifier(column_name) + ','
             if force2d:
-                query +=     'ST_Force2D(m.' + common.escape_identifier(geom_cols[0]) + ') '
-                query +=         'AS ' + common.escape_identifier(geom_cols[0]) + ','
+                geom_col2d_esc = 'ST_Force2D(m.' + common.escape_identifier(geom_cols[0]) + ')'
+                query +=     geom_col2d_esc + ' AS ' + common.escape_identifier(geom_cols[0]) + ','
             else:
-                query +=     'm.' + common.escape_identifier(geom_cols[0]) + ','
+                geom_col2d_esc = 'm.' + common.escape_identifier(geom_cols[0])
+                query +=     geom_col2d_esc + ','
+            if d == 2:
+                query +=     'ST_Area(' + geom_col2d_esc +') AS geom_area,'
+                query +=     'ST_Perimeter(' + geom_col2d_esc +') AS geom_perimeter,'
+            elif d == 1:
+                query +=     'ST_Length(' + geom_col2d_esc +') AS geom_length,'
             query +=         '%s AS layer_group,%s AS layer '
             query +=      'FROM ' + common.escape_identifier(SCHEMA_NAME) + '.'
             query +=                common.escape_identifier(layername) + ' m '