open_pg.c

Go to the documentation of this file.

/*!
   \file lib/vector/Vlib/open_pg.c
 
   \brief Vector library - Open PostGIS layer as vector map layer
 
   Higher level functions for reading/writing/manipulating vectors.
 
   (C) 2011-2013 by the GRASS Development Team
 
   This program is free software under the GNU General Public License
   (>=v2). Read the file COPYING that comes with GRASS for details.
 
   \author Martin Landa <landa.martin gmail.com>
 */
 
#include <string.h>
#include <stdlib.h>
 
#include <grass/vector.h>
#include <grass/dbmi.h>
#include <grass/glocale.h>
 
#include "local_proto.h"
 
#ifdef HAVE_POSTGRES
#include "pg_local_proto.h"
 
struct line_data {
    int id;
    int fid;
    int start_node;
    int end_node;
    int left_face;
    int right_face;
    char *wkb_geom;
};
 
static char *get_key_column(struct Format_info_pg *);
static SF_FeatureType ftype_from_string(const char *);
static int drop_table(struct Format_info_pg *);
static void connect_db(struct Format_info_pg *);
static int check_topo(struct Format_info_pg *, struct Plus_head *);
static int parse_bbox(const char *, struct bound_box *);
static struct P_node *read_p_node(struct Plus_head *, int, int, const char *,
                                  const char *, const char *,
                                  struct Format_info_pg *, int);
static struct P_line *read_p_line(struct Plus_head *, int,
                                  const struct line_data *, int,
                                  struct Format_info_cache *);
static struct P_area *read_p_area(struct Plus_head *, int, const char *, int,
                                  const char *);
static struct P_isle *read_p_isle(struct Plus_head *, int, const char *, int);
static void notice_processor(void *, const char *);
static char **scan_array(const char *);
static int remap_node(const struct Format_info_offset *, int);
static int remap_line(const struct Plus_head *, off_t, int);
 
#define NOPG_UNUSED
#else
#define NOPG_UNUSED UNUSED
#endif
 
/*!
   \brief Open vector map - PostGIS feature table on non-topological
   level
 
   \param[in,out] Map pointer to Map_info structure
   \param update TRUE for write mode, otherwise read-only
 
   \return 0 success
   \return -1 error
 */
int V1_open_old_pg(struct Map_info *Map NOPG_UNUSED, int update NOPG_UNUSED)
{
#ifdef HAVE_POSTGRES
    int found;
 
    char stmt[DB_SQL_MAX];
 
    PGresult *res;
 
    struct Format_info_pg *pg_info;
 
    G_debug(2, "V1_open_old_pg(): update = %d", update);
 
    pg_info = &(Map->fInfo.pg);
    if (!pg_info->conninfo) {
        G_warning(_("Connection string not defined"));
        return -1;
    }
 
    if (!pg_info->table_name) {
        G_warning(_("PostGIS feature table not defined"));
        return -1;
    }
 
    G_debug(1, "V1_open_old_pg(): conninfo='%s' table='%s'", pg_info->conninfo,
            pg_info->table_name);
 
    /* connect database */
    if (!pg_info->conn)
        connect_db(pg_info);
 
    /* get fid and geometry column */
    sprintf(stmt,
            "SELECT f_geometry_column, coord_dimension, srid, type "
            "FROM geometry_columns WHERE f_table_schema = '%s' AND "
            "f_table_name = '%s'",
            pg_info->schema_name, pg_info->table_name);
    G_debug(2, "SQL: %s", stmt);
 
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK)
        G_fatal_error("%s\n%s", _("No feature tables found in database."),
                      PQresultErrorMessage(res));
 
    found = PQntuples(res) > 0 ? TRUE : FALSE;
    if (found) {
        /* geometry column */
        pg_info->geom_column = G_store(PQgetvalue(res, 0, 0));
        G_debug(3, "\t-> table = %s column = %s", pg_info->table_name,
                pg_info->geom_column);
        /* fid column */
        pg_info->fid_column = get_key_column(pg_info);
        /* coordinates dimension */
        pg_info->coor_dim = atoi(PQgetvalue(res, 0, 1));
        /* SRS ID */
        pg_info->srid = atoi(PQgetvalue(res, 0, 2));
        /* feature type */
        pg_info->feature_type = ftype_from_string(PQgetvalue(res, 0, 3));
    }
    PQclear(res);
 
    /* fid -1 for empty / -2 for full cache */
    pg_info->cache.fid = pg_info->cache.ctype != CACHE_MAP ? -1 : -2;
 
    if (!found) {
        G_warning(_("Feature table <%s> not found in 'geometry_columns'"),
                  pg_info->table_name);
        return 0; /* avoid calling G_fatal_error() */
    }
 
    /* check for topo schema */
    check_topo(pg_info, &(Map->plus));
 
    return 0;
#else
    G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
    return -1;
#endif
}
 
/*!
   \brief Open vector map - PostGIS feature table on topological level
 
   Simple feature access:
   - open feature index file
   PostGIS Topology:
   - check if topological schema exists
 
   \param[in,out] Map pointer to Map_info structure
 
   \return 0 success
   \return -1 error
 */
int V2_open_old_pg(struct Map_info *Map NOPG_UNUSED)
{
#ifdef HAVE_POSTGRES
    struct Format_info_pg *pg_info;
 
    PGresult *res;
 
    G_debug(3, "V2_open_old_pg(): name = %s mapset = %s", Map->name,
            Map->mapset);
 
    pg_info = &(Map->fInfo.pg);
 
    if (pg_info->toposchema_name) {
        char stmt[DB_SQL_MAX];
 
        /* get topo schema id */
        sprintf(stmt, "SELECT id FROM topology.topology WHERE name = '%s'",
                pg_info->toposchema_name);
        res = PQexec(pg_info->conn, stmt);
        if (!res || PQresultStatus(res) != PGRES_TUPLES_OK) {
            G_warning("%s\n%s", _("Topology schema not found."),
                      PQresultErrorMessage(res));
            if (res)
                PQclear(res);
            return -1;
        }
        pg_info->toposchema_id = atoi(PQgetvalue(res, 0, 0));
        PQclear(res);
    }
    else {
        /* fidx file needed only for simple features access */
        if (Vect_open_fidx(Map, &(pg_info->offset)) != 0) {
            G_warning(
                _("Unable to open feature index file for vector map <%s>"),
                Vect_get_full_name(Map));
            G_zero(&(pg_info->offset), sizeof(struct Format_info_offset));
        }
    }
    return 0;
#else
    G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
    return -1;
#endif
}
 
/*!
   \brief Prepare PostGIS database for creating new feature table
   (level 1)
 
   New PostGIS table is created when writing features by
   Vect_wrile_line().
 
   \param[out] Map pointer to Map_info structure
   \param name name of PostGIS feature table to create
   \param with_z WITH_Z for 3D vector data otherwise WITHOUT_Z
 
   \return 0 success
   \return -1 error
 */
int V1_open_new_pg(struct Map_info *Map NOPG_UNUSED,
                   const char *name NOPG_UNUSED, int with_z NOPG_UNUSED)
{
#ifdef HAVE_POSTGRES
    char stmt[DB_SQL_MAX];
 
    struct Format_info_pg *pg_info;
 
    PGresult *res;
 
    G_debug(2, "V1_open_new_pg(): name = %s with_z = %d", name, with_z);
 
    pg_info = &(Map->fInfo.pg);
    if (!pg_info->conninfo) {
        G_warning(_("Connection string not defined"));
        return -1;
    }
 
    if (!pg_info->table_name) {
        G_warning(_("PostGIS feature table not defined"));
        return -1;
    }
 
    G_debug(1, "V1_open_new_pg(): conninfo='%s' table='%s'", pg_info->conninfo,
            pg_info->table_name);
 
    /* connect database & get DB name */
    connect_db(pg_info);
 
    /* if schema not defined, use 'public' */
    if (!pg_info->schema_name)
        pg_info->schema_name = G_store("public");
 
    /* if fid_column not defined, use 'fid' */
    if (!pg_info->fid_column)
        pg_info->fid_column = G_store(GV_PG_FID_COLUMN);
 
    /* if geom_column not defined, use 'geom' */
    if (!pg_info->geom_column)
        pg_info->geom_column = G_store(GV_PG_GEOMETRY_COLUMN);
 
    /* check if feature table already exists */
    sprintf(stmt,
            "SELECT * FROM pg_tables "
            "WHERE schemaname = '%s' AND tablename = '%s'",
            pg_info->schema_name, pg_info->table_name);
    G_debug(2, "SQL: %s", stmt);
 
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK)
        G_fatal_error("%s\n%s", _("No feature tables found in database."),
                      PQresultErrorMessage(res));
 
    if (PQntuples(res) > 0) {
        /* table found */
        if (G_get_overwrite()) {
            G_warning(_("PostGIS layer <%s.%s> already exists and will be "
                        "overwritten"),
                      pg_info->schema_name, pg_info->table_name);
            if (drop_table(pg_info) == -1) {
                G_warning(_("Unable to delete PostGIS layer <%s>"),
                          pg_info->table_name);
                return -1;
            }
        }
        else {
            G_warning(
                _("PostGIS layer <%s.%s> already exists in database '%s'"),
                pg_info->schema_name, pg_info->table_name, pg_info->db_name);
            return -1;
        }
    }
 
    /* no feature in cache */
    pg_info->cache.fid = -1;
 
    /* unknown feature type */
    pg_info->feature_type = SF_GEOMETRY;
 
    PQclear(res);
 
    return 0;
#else
    G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
    return -1;
#endif
}
 
/*!
   \brief Read full-topology for PostGIS links
 
   Note: Only 2D topological primitives are currently supported
 
   \param[in,out] Map pointer to Map_info structure
   \param head_only TRUE to read only header
   \param update TRUE to clean GRASS topology in update mode
 
   \return 0 on success
   \return 1 topology layer does not exist
   \return -1 on error
 */
int Vect__open_topo_pg(struct Map_info *Map NOPG_UNUSED,
                       int head_only NOPG_UNUSED, int update NOPG_UNUSED)
{
#ifdef HAVE_POSTGRES
    int ret;
    struct Plus_head *plus;
    struct Format_info_pg *pg_info;
 
    Map->open = VECT_OPEN_CODE; /* needed by load_plus */
 
    plus = &(Map->plus);
    pg_info = &(Map->fInfo.pg);
 
    /* check for topo schema */
    if (check_topo(pg_info, plus) != 0)
        return 1;
 
    /* free and init plus structure, update spatial and category
     * indices  */
    dig_init_plus(plus);
    plus->Spidx_new = TRUE; /* force building spatial and category indices */
    plus->update_cidx = TRUE;
    Map->support_updated = FALSE; /* don't write support files */
 
    ret = Vect__load_plus_pg(Map, head_only);
 
    if (update)
        Vect__clean_grass_db_topo(pg_info);
 
    plus->cidx_up_to_date = TRUE; /* category index built from topo */
 
    return ret;
#else
    G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
    return -1;
#endif
}
 
#ifdef HAVE_POSTGRES
/*!
   \brief Get key column for feature table
 
   \param pg_info pointer to Format_info_pg
 
   \return string buffer with key column name
   \return NULL on missing key column
 */
char *get_key_column(struct Format_info_pg *pg_info)
{
    char *key_column;
    char stmt[DB_SQL_MAX];
 
    PGresult *res;
 
    sprintf(stmt,
            "SELECT kcu.column_name "
            "FROM INFORMATION_SCHEMA.TABLES t "
            "LEFT JOIN INFORMATION_SCHEMA.TABLE_CONSTRAINTS tc "
            "ON tc.table_catalog = t.table_catalog "
            "AND tc.table_schema = t.table_schema "
            "AND tc.table_name = t.table_name "
            "AND tc.constraint_type = 'PRIMARY KEY' "
            "LEFT JOIN INFORMATION_SCHEMA.KEY_COLUMN_USAGE kcu "
            "ON kcu.table_catalog = tc.table_catalog "
            "AND kcu.table_schema = tc.table_schema "
            "AND kcu.table_name = tc.table_name "
            "AND kcu.constraint_name = tc.constraint_name "
            "WHERE t.table_schema = '%s' AND t.table_name = '%s'",
            pg_info->schema_name, pg_info->table_name);
    G_debug(2, "SQL: %s", stmt);
 
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK || PQntuples(res) != 1 ||
        strlen(PQgetvalue(res, 0, 0)) < 1) {
        G_warning(_("No key column detected."));
        if (res)
            PQclear(res);
        return NULL;
    }
    key_column = G_store(PQgetvalue(res, 0, 0));
 
    PQclear(res);
 
    return key_column;
}
 
/*!
   \brief Get simple feature type from string
 
   \param type string
 
   \return SF type
 */
SF_FeatureType ftype_from_string(const char *type)
{
    SF_FeatureType sf_type;
 
    if (G_strcasecmp(type, "POINT") == 0)
        return SF_POINT;
    else if (G_strcasecmp(type, "LINESTRING") == 0)
        return SF_LINESTRING;
    else if (G_strcasecmp(type, "POLYGON") == 0)
        return SF_POLYGON;
    else if (G_strcasecmp(type, "MULTIPOINT") == 0)
        return SF_MULTIPOINT;
    else if (G_strcasecmp(type, "MULTILINESTRING") == 0)
        return SF_MULTILINESTRING;
    else if (G_strcasecmp(type, "MULTIPOLYGON") == 0)
        return SF_MULTIPOLYGON;
    else if (G_strcasecmp(type, "GEOMETRYCOLLECTION") == 0)
        return SF_GEOMETRYCOLLECTION;
    else
        return SF_GEOMETRY;
 
    G_debug(3, "ftype_from_string(): type='%s' -> %d", type, sf_type);
 
    return sf_type;
}
 
/*!
   \brief Drop feature table and topology schema if exists
 
   \param pg_info pointer to Format_info_pg
 
   \return -1 on error
   \return 0 on success
 */
int drop_table(struct Format_info_pg *pg_info)
{
    int i;
    char stmt[DB_SQL_MAX];
    char *topo_schema;
 
    PGresult *result, *result_drop;
 
    /* check if topology schema exists */
    sprintf(stmt,
            "SELECT COUNT(*) FROM pg_tables WHERE schemaname = 'topology'");
    if (Vect__execute_get_value_pg(pg_info->conn, stmt) != 0) {
        /* drop topology schema(s) related to the feature table */
        sprintf(stmt,
                "SELECT t.name FROM topology.layer AS l JOIN "
                "topology.topology AS t ON l.topology_id = t.id "
                "WHERE l.table_name = '%s'",
                pg_info->table_name);
        G_debug(2, "SQL: %s", stmt);
 
        result = PQexec(pg_info->conn, stmt);
        if (!result || PQresultStatus(result) != PGRES_TUPLES_OK) {
            G_warning(_("Execution failed: %s"), PQerrorMessage(pg_info->conn));
            PQclear(result);
            return -1;
        }
        for (i = 0; i < PQntuples(result); i++) {
            topo_schema = PQgetvalue(result, i, 0);
            sprintf(stmt, "SELECT topology.DropTopology('%s')", topo_schema);
            G_debug(2, "SQL: %s", stmt);
 
            result_drop = PQexec(pg_info->conn, stmt);
            if (!result_drop || PQresultStatus(result_drop) != PGRES_TUPLES_OK)
                G_warning(_("Execution failed: %s"),
                          PQerrorMessage(pg_info->conn));
 
            G_verbose_message(_("PostGIS topology schema <%s> dropped"),
                              topo_schema);
            PQclear(result_drop);
        }
        PQclear(result);
    }
 
    /* drop feature table */
    sprintf(stmt, "DROP TABLE \"%s\".\"%s\"", pg_info->schema_name,
            pg_info->table_name);
    G_debug(2, "SQL: %s", stmt);
 
    if (Vect__execute_pg(pg_info->conn, stmt) == -1) {
        return -1;
    }
 
    return 0;
}
 
/*!
   \brief Establish PG connection (pg_info->conninfo)
 
   Check if DB is spatial as defined by PostGIS.
 
   \param pg_info pointer to Format_info_pg
 */
void connect_db(struct Format_info_pg *pg_info)
{
    char stmt[DB_SQL_MAX];
 
    /* check if connection string already contains user/passwd */
    if (!strstr(pg_info->conninfo, "user")) {
        char dbname[GNAME_MAX];
        char *p;
        const char *user, *passwd, *host, *port;
 
        dbname[0] = '\0';
        p = strstr(pg_info->conninfo, "dbname");
        if (p) {
            int i;
 
            p += strlen("dbname") + 1; /* dbname= */
 
            for (i = 0; *p && *p != ' '; i++, p++)
                dbname[i] = *p;
        }
 
        /* try connection settings for given database first, then try
         * any settings defined for pg driver */
        db_get_login("pg", dbname, &user, &passwd, &host, &port);
        /* any settings defined for pg driver disabled - can cause
           problems when running multiple local/remote db clusters
           if (strlen(dbname) > 0 && !user && !passwd)
           db_get_login("pg", NULL, &user, &passwd, &host, &port);
         */
        if (user || passwd || host || port) {
            char conninfo[DB_SQL_MAX];
 
            sprintf(conninfo, "%s", pg_info->conninfo);
            if (user) {
                strcat(conninfo, " user=");
                strcat(conninfo, user);
            }
            if (passwd) {
                strcat(conninfo, " password=");
                strcat(conninfo, passwd);
            }
            if (host) {
                strcat(conninfo, " host=");
                strcat(conninfo, host);
            }
            if (port) {
                strcat(conninfo, " port=");
                strcat(conninfo, port);
            }
            G_free(pg_info->conninfo);
            pg_info->conninfo = G_store(conninfo);
        }
    }
 
    pg_info->conn = PQconnectdb(pg_info->conninfo);
    G_debug(1, "   PQconnectdb(): %s", pg_info->conninfo);
    if (PQstatus(pg_info->conn) == CONNECTION_BAD)
        G_fatal_error("%s\n%s",
                      _("Connection to PostgreSQL database failed. "
                        "Try to set up username/password by db.login."),
                      PQerrorMessage(pg_info->conn));
 
    /* get DB name */
    pg_info->db_name = G_store(PQdb(pg_info->conn));
    if (!pg_info->db_name)
        G_warning(_("Unable to get database name"));
 
    sprintf(
        stmt,
        "SELECT COUNT(*) FROM pg_tables WHERE tablename = 'spatial_ref_sys'");
    if (Vect__execute_get_value_pg(pg_info->conn, stmt) != 1) {
        PQfinish(pg_info->conn);
        G_fatal_error(_("<%s> is not PostGIS database. DB table "
                        "'spatial_ref_sys' not found."),
                      pg_info->db_name ? pg_info->db_name : pg_info->conninfo);
    }
 
    if (pg_info->toposchema_name) {
        /* check if topology schema exists */
        sprintf(stmt,
                "SELECT COUNT(*) FROM pg_tables WHERE schemaname = 'topology'");
        if (Vect__execute_get_value_pg(pg_info->conn, stmt) == 0) {
            PQfinish(pg_info->conn);
            G_fatal_error(
                _("PostGIS Topology extension not found in the database <%s>"),
                pg_info->db_name);
        }
    }
 
    /* print notice messages only on verbose level */
    PQsetNoticeProcessor(pg_info->conn, notice_processor, NULL);
}
 
/*!
   \brief Check for topology schema (pg_info->toposchema_name)
 
   \param pg_info pointer to Format_info_pg
 
   \return 0 schema exists
   \return 1 schema doesn't exists
 */
int check_topo(struct Format_info_pg *pg_info, struct Plus_head *plus)
{
    char stmt[DB_SQL_MAX];
 
    PGresult *res;
 
    /* connect database */
    if (!pg_info->conn)
        connect_db(pg_info);
 
    if (pg_info->toposchema_name)
        return 0;
 
    /* check if topology layer/schema exists */
    sprintf(stmt,
            "SELECT t.id,t.name,t.hasz,l.feature_column FROM topology.layer "
            "AS l JOIN topology.topology AS t ON l.topology_id = t.id "
            "WHERE schema_name = '%s' AND table_name = '%s'",
            pg_info->schema_name, pg_info->table_name);
    G_debug(2, "SQL: %s", stmt);
 
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK || PQntuples(res) != 1) {
        G_debug(1, "Topology layers for '%s.%s' not found (%s)",
                pg_info->schema_name, pg_info->table_name,
                PQerrorMessage(pg_info->conn));
        if (res)
            PQclear(res);
        return 1;
    }
 
    pg_info->toposchema_id = atoi(PQgetvalue(res, 0, 0));
    pg_info->toposchema_name = G_store(PQgetvalue(res, 0, 1));
    pg_info->topogeom_column = G_store(PQgetvalue(res, 0, 3));
 
    /* check extra GRASS tables */
    sprintf(stmt,
            "SELECT COUNT(*) FROM pg_tables WHERE schemaname = '%s' "
            "AND tablename LIKE '%%_grass'",
            pg_info->toposchema_name);
    if (Vect__execute_get_value_pg(pg_info->conn, stmt) != TOPO_TABLE_NUM)
        pg_info->topo_geo_only = TRUE;
 
    G_debug(
        1,
        "PostGIS topology detected: schema = %s column = %s topo_geo_only = %d",
        pg_info->toposchema_name, pg_info->topogeom_column,
        pg_info->topo_geo_only);
 
    /* check for 3D */
    if (strcmp(PQgetvalue(res, 0, 2), "t") == 0)
        plus->with_z = WITH_Z;
    PQclear(res);
 
    return 0;
}
 
/*!
   \brief Parse BBOX string
 
   \param value string buffer
   \param[out] bbox pointer to output bound_box struct
 
   \return 0 on success
   \return -1 on error
 */
int parse_bbox(const char *value, struct bound_box *bbox)
{
    unsigned int i;
    size_t length, prefix_length;
    char **tokens, **tokens_coord, *coord;
 
    if (strlen(value) < 1) {
        G_warning(_("Empty bounding box"));
        return -1;
    }
 
    prefix_length = strlen("box3d(");
    if (G_strncasecmp(value, "box3d(", prefix_length) != 0)
        return -1;
 
    /* strip off "bbox3d(...)" */
    length = strlen(value);
    coord = G_malloc(length - prefix_length);
    for (i = prefix_length; i < length; i++)
        coord[i - prefix_length] = value[i];
    coord[length - prefix_length - 1] = '\0';
 
    tokens = G_tokenize(coord, ",");
    G_free(coord);
 
    if (G_number_of_tokens(tokens) != 2) {
        G_free_tokens(tokens);
        return -1;
    }
 
    /* parse bbox LL corner */
    tokens_coord = G_tokenize(tokens[0], " ");
    if (G_number_of_tokens(tokens_coord) != 3) {
        G_free_tokens(tokens);
        G_free_tokens(tokens_coord);
    }
    bbox->W = atof(tokens_coord[0]);
    bbox->S = atof(tokens_coord[1]);
    bbox->B = atof(tokens_coord[2]);
 
    G_free_tokens(tokens_coord);
 
    /* parse bbox UR corner */
    tokens_coord = G_tokenize(tokens[1], " ");
    if (G_number_of_tokens(tokens_coord) != 3) {
        G_free_tokens(tokens);
        G_free_tokens(tokens_coord);
    }
    bbox->E = atof(tokens_coord[0]);
    bbox->N = atof(tokens_coord[1]);
    bbox->T = atof(tokens_coord[2]);
 
    G_free_tokens(tokens_coord);
    G_free_tokens(tokens);
 
    return 0;
}
 
/*!
   \brief Read P_node structure
 
   See dig_Rd_P_node() for reference.
 
   \param plus pointer to Plus_head structure
   \param n index (starts at 1)
   \param id node id (table "node")
   \param wkb_data geometry data (wkb)
   \param lines_data lines array or NULL
   \param angles_data angles array or NULL
   \param pg_info pointer to Format_info_pg sttucture
   \param geom_only TRUE to read node's geometry
 
   \return pointer to new P_node struct
   \return NULL on error
 */
struct P_node *read_p_node(struct Plus_head *plus, int n, int id,
                           const char *wkb_data, const char *lines_data,
                           const char *angles_data,
                           struct Format_info_pg *pg_info, int geom_only)
{
    int i, cnt;
    char **lines, **angles;
 
    struct P_node *node;
    struct line_pnts *points;
 
    PGresult *res;
 
    /* get lines connected to the node */
    cnt = 0;
    lines = angles = NULL;
    if (!geom_only) {
        if (!lines_data && !angles_data) { /* pg_info->topo_geo_only == TRUE */
            char stmt[DB_SQL_MAX];
 
            sprintf(
                stmt,
                "SELECT edge_id,'s' as node,"
                "ST_Azimuth(ST_StartPoint(geom), ST_PointN(geom, 2)) AS angle"
                " FROM \"%s\".edge WHERE start_node = %d UNION ALL "
                "SELECT edge_id,'e' as node,"
                "ST_Azimuth(ST_EndPoint(geom), ST_PointN(geom, "
                "ST_NumPoints(geom) - 1)) AS angle"
                " FROM \"%s\".edge WHERE end_node = %d"
                " ORDER BY angle DESC",
                pg_info->toposchema_name, id, pg_info->toposchema_name, id);
            G_debug(2, "SQL: %s", stmt);
            res = PQexec(pg_info->conn, stmt);
            if (!res || PQresultStatus(res) != PGRES_TUPLES_OK) {
                G_warning(
                    _("Inconsistency in topology: unable to read node %d"), id);
                if (res)
                    PQclear(res);
                return NULL;
            }
            cnt = PQntuples(res);
        }
        else { /* pg_info->topo_geo_only != TRUE */
            lines = scan_array(lines_data);
            angles = scan_array(angles_data);
 
            cnt = G_number_of_tokens(lines);
            if (cnt != G_number_of_tokens(angles))
                return NULL; /* 'lines' and 'angles' array must have the same
                                size */
        }
 
        if (cnt == 0) { /* dead */
            plus->Node[n] = NULL;
            return NULL;
        }
    }
 
    node = dig_alloc_node();
    node->n_lines = cnt;
    G_debug(4, "read_p_node(): id = %d, n_lines = %d", id, cnt);
 
    if (!geom_only) {
        if (dig_node_alloc_line(node, node->n_lines) == -1)
            return NULL;
 
        /* lines / angles */
        if (lines) {
            for (i = 0; i < node->n_lines; i++) {
                node->lines[i] = atoi(lines[i]);
                node->angles[i] = atof(angles[i]);
 
                G_debug(5, "\tline = %d angle = %f", node->lines[i],
                        node->angles[i]);
            }
 
            G_free_tokens(lines);
            G_free_tokens(angles);
        }
        else {
            for (i = 0; i < node->n_lines; i++) {
                node->lines[i] = atoi(PQgetvalue(res, i, 0));
                if (strcmp(PQgetvalue(res, i, 1), "s") != 0) {
                    /* end node */
                    node->lines[i] *= -1;
                }
                node->angles[i] = M_PI / 2 - atof(PQgetvalue(res, i, 2));
                /* angles range <-PI; PI> */
                if (node->angles[i] > M_PI)
                    node->angles[i] = node->angles[i] - 2 * M_PI;
                if (node->angles[i] < -1.0 * M_PI)
                    node->angles[i] = node->angles[i] + 2 * M_PI;
                G_debug(5, "\tline = %d angle = %f", node->lines[i],
                        node->angles[i]);
            }
            PQclear(res);
        }
    }
 
    /* get node coordinates */
    if (SF_POINT !=
        Vect__cache_feature_pg(wkb_data, FALSE, FALSE, &(pg_info->cache), NULL))
        G_warning(_("Inconsistency in topology: node %d - unexpected feature "
                    "type %d"),
                  n, pg_info->cache.sf_type);
 
    points = pg_info->cache.lines[0];
    node->x = points->x[0];
    node->y = points->y[0];
    if (plus->with_z)
        node->z = points->z[0];
    else
        node->z = 0.0;
 
    /* update spatial index */
    if (plus->Spidx_new)
        dig_spidx_add_node(plus, n, node->x, node->y, node->z);
 
    if (plus->uplist.do_uplist)
        /* collect updated nodes if requested */
        dig_node_add_updated(plus, n);
 
    plus->Node[n] = node;
 
    return node;
}
 
/*!
   \brief Read P_line structure
 
   See dig_Rd_P_line() for reference.
 
   Supported feature types:
   - GV_POINT
   - GV_LINE
   - GV_BOUNDARY
 
   \param plus pointer to Plus_head structure
   \param n index (starts at 1)
   \param data edge data (id, start/end node, left/right face, ...)
   \param topo_geo_only TRUE for topo-geo-only mode
   \param cache pointer to Format_info_cache structure
 
   \return pointer to P_line struct
   \return NULL on error
 */
struct P_line *read_p_line(struct Plus_head *plus, int n,
                           const struct line_data *data, int topo_geo_only,
                           struct Format_info_cache *cache)
{
    int tp, cat;
    struct P_line *line;
 
    if (data->start_node == 0 && data->end_node == 0) {
        if (data->left_face == 0)
            tp = GV_POINT;
        else
            tp = GV_CENTROID;
    }
    else if (data->left_face == 0 && data->right_face == 0) {
        tp = GV_LINE;
    }
    else {
        tp = GV_BOUNDARY;
    }
 
    if (tp == 0) { /* dead ??? */
        plus->Line[n] = NULL;
        return NULL;
    }
 
    line = dig_alloc_line();
 
    /* type & offset ( = id) */
    line->type = tp;
    line->offset = data->id;
    G_debug(4, "read_p_line(): id/offset = %d type = %d", data->id, line->type);
 
    /* topo */
    if (line->type == GV_POINT) {
        line->topo = NULL;
    }
    else {
        line->topo = dig_alloc_topo(line->type);
 
        if ((line->type & GV_LINES) &
            (data->start_node < 0 || data->end_node < 0))
            return NULL;
 
        /* lines */
        if (line->type == GV_LINE) {
            struct P_topo_l *topo = (struct P_topo_l *)line->topo;
 
            topo->N1 = data->start_node;
            topo->N2 = data->end_node;
        }
        /* boundaries */
        else if (line->type == GV_BOUNDARY) {
            struct P_topo_b *topo = (struct P_topo_b *)line->topo;
 
            topo->N1 = data->start_node;
            topo->N2 = data->end_node;
 
            if (topo_geo_only) {
                topo->left = topo->right = 0;
            }
            else {
                topo->left = data->left_face;
                topo->right = data->right_face;
            }
        }
        /* centroids */
        else if (line->type == GV_CENTROID) {
            struct P_topo_c *topo = (struct P_topo_c *)line->topo;
 
            topo->area = data->left_face;
        }
    }
 
    Vect__cache_feature_pg(data->wkb_geom, FALSE, tp, cache, NULL);
    cat = cache->lines_cats[cache->lines_num - 1] =
        data->fid > 0 ? data->fid : -1;
 
    /* update spatial index */
    if (plus->Spidx_new) {
        struct line_pnts *points;
        struct bound_box box;
 
        points = cache->lines[cache->lines_num - 1];
        dig_line_box(points, &box);
        dig_spidx_add_line(plus, n, &box);
    }
 
    /* update category index */
    if (plus->update_cidx)
        dig_cidx_add_cat(plus, cat > 0 ? 1 : 0, cat > 0 ? cat : 0, n, tp);
 
    if (plus->uplist.do_uplist) {
        /* collect updated lines if requested */
        dig_line_add_updated(plus, n, line->offset);
    }
 
    plus->Line[n] = line;
 
    return line;
}
 
/*!
   \brief Read P_area structure
 
   \param plus pointer to Plus_head structure
   \param n index (starts at 1)
   \param lines_data lines array (see P_area struct)
   \param centroid centroid id (see P_area struct)
   \param isles_data lines array (see P_area struct)
 
   \return pointer to P_area struct
   \return NULL on error
 */
struct P_area *read_p_area(struct Plus_head *plus, int n,
                           const char *lines_data, int centroid,
                           const char *isles_data)
{
    int i;
    int nlines, nisles;
    char **lines, **isles;
 
    struct P_area *area;
 
    lines = scan_array(lines_data);
    nlines = G_number_of_tokens(lines);
    isles = scan_array(isles_data);
    nisles = G_number_of_tokens(isles);
 
    if (nlines < 1) {
        G_warning(_("Area %d without boundary detected"), n);
        return NULL;
    }
 
    G_debug(3, "read_p_area(): n = %d nlines = %d nisles = %d", n, nlines,
            nisles);
 
    /* allocate area */
    area = dig_alloc_area();
    dig_area_alloc_line(area, nlines);
    dig_area_alloc_isle(area, nisles);
 
    /* set lines */
    area->n_lines = nlines;
    for (i = 0; i < nlines; i++) {
        area->lines[i] = atoi(lines[i]);
    }
 
    /* set isles */
    area->n_isles = nisles;
    for (i = 0; i < nisles; i++) {
        area->isles[i] = atoi(isles[i]);
    }
 
    /* set centroid */
    area->centroid = remap_line(plus, centroid, GV_CENTROID);
 
    G_free_tokens(lines);
    G_free_tokens(isles);
 
    plus->Area[n] = area;
 
    return area;
}
 
/*!
   \brief Read P_isle structure
 
   \param plus pointer to Plus_head structure
   \param n index (starts at 1)
   \param lines_data lines array (see P_isle struct)
   \param area area id (see P_isle struct)
 
   \return pointer to P_isle struct
   \return NULL on error
 */
struct P_isle *read_p_isle(struct Plus_head *plus, int n,
                           const char *lines_data, int area)
{
    int i;
    int nlines;
    char **lines;
 
    struct P_isle *isle;
 
    lines = scan_array(lines_data);
    nlines = G_number_of_tokens(lines);
 
    if (nlines < 1) {
        G_warning(_("Isle %d without boundary detected"), n);
        return NULL;
    }
 
    G_debug(3, "read_p_isle(): n = %d nlines = %d", n, nlines);
 
    /* allocate isle */
    isle = dig_alloc_isle();
    dig_isle_alloc_line(isle, nlines);
 
    /* set lines */
    isle->n_lines = nlines;
    for (i = 0; i < nlines; i++) {
        isle->lines[i] = atoi(lines[i]);
    }
 
    /* set area */
    isle->area = area;
 
    G_free_tokens(lines);
 
    plus->Isle[n] = isle;
 
    return isle;
}
 
/*!
   \brief Read topo from PostGIS topology schema -- header info only
 
   \param[in,out] plus pointer to Plus_head struct
 
   \return 0 on success
   \return -1 on error
 */
int Vect__load_plus_head(struct Map_info *Map)
{
    char stmt[DB_SQL_MAX];
 
    struct Format_info_pg *pg_info;
    struct Plus_head *plus;
 
    PGresult *res;
 
    plus = &(Map->plus);
    pg_info = &(Map->fInfo.pg);
 
    plus->off_t_size = -1;
 
    /* get map bounding box
       first try to get info from 'topology.grass' table */
    sprintf(stmt, "SELECT %s FROM \"%s\".\"%s\" WHERE %s = %d", TOPO_BBOX,
            TOPO_SCHEMA, TOPO_TABLE, TOPO_ID, pg_info->toposchema_id);
    G_debug(2, "SQL: %s", stmt);
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK || PQntuples(res) != 1) {
        PQclear(res);
 
        /* otherwise try to calculate bbox from TopoGeometry elements */
        sprintf(stmt, "SELECT ST_3DExtent(%s) FROM \"%s\".\"%s\"",
                pg_info->topogeom_column, pg_info->schema_name,
                pg_info->table_name);
        G_debug(2, "SQL: %s", stmt);
        res = PQexec(pg_info->conn, stmt);
        if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
            PQntuples(res) != 1 || strlen(PQgetvalue(res, 0, 0)) < 1) {
            G_warning(_("Unable to get map bounding box from topology"));
            PQclear(res);
            return -1;
        }
    }
 
    if (parse_bbox(PQgetvalue(res, 0, 0), &(plus->box)) != 0) {
        G_warning(_("Unable to parse map bounding box:\n%s"),
                  PQgetvalue(res, 0, 0));
        return -1;
    }
    PQclear(res);
 
    /* get number of topological elements */
 
    /* nodes
       note: isolated nodes are registered in GRASS Topology model */
    sprintf(stmt,
            "SELECT COUNT(DISTINCT node) FROM (SELECT start_node AS node "
            "FROM \"%s\".edge GROUP BY start_node UNION ALL SELECT end_node "
            "AS node FROM \"%s\".edge GROUP BY end_node) AS foo",
            pg_info->toposchema_name, pg_info->toposchema_name);
    plus->n_nodes = Vect__execute_get_value_pg(pg_info->conn, stmt);
    if (!pg_info->topo_geo_only) {
        int n_nodes;
 
        /* check nodes consistency */
        sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".%s",
                pg_info->toposchema_name, TOPO_TABLE_NODE);
        n_nodes = Vect__execute_get_value_pg(pg_info->conn, stmt);
        if (n_nodes != plus->n_nodes) {
            G_warning(_("Different number of nodes detected (%d, %d)"),
                      plus->n_nodes, n_nodes);
            return -1;
        }
    }
    G_debug(3, "Vect_open_topo_pg(): n_nodes=%d", plus->n_nodes);
 
    /* lines (edges in PostGIS Topology model) */
    sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".edge", pg_info->toposchema_name);
    /* + isolated nodes as points
       + centroids */
    plus->n_lines = Vect__execute_get_value_pg(pg_info->conn, stmt);
 
    /* areas (faces with face_id > 0 in PostGIS Topology model) */
    sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".face WHERE face_id > 0",
            pg_info->toposchema_name);
    plus->n_areas = Vect__execute_get_value_pg(pg_info->conn, stmt);
    if (!pg_info->topo_geo_only) {
        int n_areas;
 
        /* check areas consistency */
        sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".%s",
                pg_info->toposchema_name, TOPO_TABLE_AREA);
        n_areas = Vect__execute_get_value_pg(pg_info->conn, stmt);
        if (n_areas != plus->n_areas) {
            G_warning(_("Different number of areas detected (%d, %d)"),
                      plus->n_areas, n_areas);
            return -1;
        }
    }
    G_debug(3, "Vect_open_topo_pg(): n_areas=%d", plus->n_areas);
 
    /* isles (faces with face_id <=0 in PostGIS Topology model)
       note: universal face is represented in GRASS Topology model as isle
       (area=0)
     */
    sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".face WHERE face_id < 0",
            pg_info->toposchema_name);
    plus->n_isles = Vect__execute_get_value_pg(pg_info->conn, stmt);
    if (!pg_info->topo_geo_only) {
        int n_isles;
 
        /* check areas consistency */
        sprintf(stmt, "SELECT COUNT(*) FROM \"%s\".%s",
                pg_info->toposchema_name, TOPO_TABLE_ISLE);
        n_isles = Vect__execute_get_value_pg(pg_info->conn, stmt);
        if (n_isles != plus->n_isles) {
            G_warning(_("Different number of areas detected (%d, %d)"),
                      plus->n_isles, n_isles);
            return -1;
        }
    }
    G_debug(3, "Vect_open_topo_pg(): n_isles=%d", plus->n_isles);
 
    /* number of features according the type */
 
    /* points */
    sprintf(stmt,
            "SELECT COUNT(*) FROM \"%s\".node WHERE containing_face "
            "IS NULL AND node_id NOT IN "
            "(SELECT node FROM (SELECT start_node AS node FROM \"%s\".edge "
            "GROUP BY start_node UNION ALL SELECT end_node AS node FROM "
            "\"%s\".edge GROUP BY end_node) AS foo)",
            pg_info->toposchema_name, pg_info->toposchema_name,
            pg_info->toposchema_name);
    plus->n_plines = Vect__execute_get_value_pg(pg_info->conn, stmt);
    G_debug(3, "Vect_open_topo_pg(): n_plines=%d", plus->n_plines);
 
    /* lines */
    sprintf(stmt,
            "SELECT COUNT(*) FROM \"%s\".edge WHERE "
            "left_face = 0 AND right_face = 0",
            pg_info->toposchema_name);
    plus->n_llines = Vect__execute_get_value_pg(pg_info->conn, stmt);
    G_debug(3, "Vect_open_topo_pg(): n_llines=%d", plus->n_llines);
 
    /* boundaries */
    sprintf(stmt,
            "SELECT COUNT(*) FROM \"%s\".edge WHERE "
            "left_face != 0 OR right_face != 0",
            pg_info->toposchema_name);
    plus->n_blines = Vect__execute_get_value_pg(pg_info->conn, stmt);
    G_debug(3, "Vect_open_topo_pg(): n_blines=%d", plus->n_blines);
 
    /* centroids */
    sprintf(stmt,
            "SELECT COUNT(*) FROM \"%s\".node WHERE containing_face "
            "IS NOT NULL AND node_id NOT IN "
            "(SELECT node FROM (SELECT start_node AS node FROM \"%s\".edge "
            "GROUP BY start_node UNION ALL SELECT end_node AS node FROM "
            "\"%s\".edge GROUP BY end_node) AS foo)",
            pg_info->toposchema_name, pg_info->toposchema_name,
            pg_info->toposchema_name);
    plus->n_clines = Vect__execute_get_value_pg(pg_info->conn, stmt);
    G_debug(3, "Vect_open_topo_pg(): n_clines=%d", plus->n_clines);
 
    /* update number of lines - add points and centroids */
    plus->n_lines += plus->n_plines + plus->n_clines;
    G_debug(3, "Vect_open_topo_pg(): n_lines=%d", plus->n_lines);
 
    return 0;
}
 
/*!
   \brief Read topo info from PostGIS topology schema
 
   \param pg_info pointer to Format_info_pg
   \param[in,out] plus pointer to Plus_head struct
   \param head_only TRUE to read only header info
 
   \return 0 on success
   \return -1 on error
 */
int Vect__load_plus_pg(struct Map_info *Map, int head_only)
{
    int i, side, line;
    char stmt[DB_SQL_MAX];
 
    struct Format_info_pg *pg_info;
    struct Plus_head *plus;
    struct P_line *Line;
    struct line_pnts *Points;
    struct ilist *List;
    struct bound_box box;
 
    PGresult *res;
 
    pg_info = &(Map->fInfo.pg);
    plus = &(Map->plus);
 
    if (Vect__load_plus_head(Map) != 0)
        return -1;
 
    if (head_only)
        return 0;
 
    Points = Vect_new_line_struct();
    List = Vect_new_list();
 
    /* read nodes (GRASS Topo)
       note: standalone nodes (ie. points/centroids) are ignored
     */
    Vect__load_map_nodes_pg(Map, FALSE);
 
    /* read lines (GRASS Topo)
       - standalone nodes -> points|centroids
       - edges -> lines/boundaries
     */
    Vect__free_cache(&(pg_info->cache));
    pg_info->cache.ctype = CACHE_MAP;
 
    Vect__load_map_lines_pg(Map);
 
    /* build areas */
    if (pg_info->topo_geo_only) {
        /* build areas for boundaries
           reset values -> build from scratch */
        plus->n_areas = plus->n_isles = 0;
        for (line = 1; line <= plus->n_lines; line++) {
            Line = plus->Line[line]; /* centroids: Line is NULL */
            if (!Line || Line->type != GV_BOUNDARY)
                continue;
 
            for (i = 0; i < 2; i++) { /* for both sides build an area/isle */
                side = i == 0 ? GV_LEFT : GV_RIGHT;
 
                G_debug(3, "Build area for line = %d, side = %d", i, side);
                Vect_build_line_area(Map, line, side);
            }
        }
    }
    else {
        int cat;
 
        /* read areas from 'area_grass' table */
        sprintf(stmt,
                "SELECT area_id,lines,centroid,isles FROM \"%s\".%s ORDER BY "
                "area_id",
                pg_info->toposchema_name, TOPO_TABLE_AREA);
        G_debug(2, "SQL: %s", stmt);
 
        res = PQexec(pg_info->conn, stmt);
        if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
            plus->n_areas != PQntuples(res)) {
            if (res)
                PQclear(res);
            return -1;
        }
 
        dig_alloc_areas(plus, plus->n_areas);
        G_zero(plus->Area, sizeof(struct P_area *) *
                               (plus->n_areas + 1)); /* index starts at 1 */
        G_debug(3, "Vect_open_topo_pg(): n_areas=%d", plus->n_areas);
 
        for (i = 0; i < plus->n_areas; i++) {
            read_p_area(plus, i + 1, (char *)PQgetvalue(res, i, 1),
                        atoi(PQgetvalue(res, i, 2)),
                        (char *)PQgetvalue(res, i, 3));
 
            if (plus->Spidx_new) {
                /* update spatial index */
                Vect_get_area_points(Map, i + 1, Points);
                dig_line_box(Points, &box);
                dig_spidx_add_area(&(Map->plus), i + 1, &box);
            }
 
            if (plus->update_cidx) {
                /* update category index */
                cat =
                    pg_info->cache.lines_cats[plus->Area[i + 1]->centroid - 1];
                dig_cidx_add_cat(plus, cat > 0 ? 1 : 0, cat > 0 ? cat : 0,
                                 i + 1, GV_AREA);
            }
        }
        PQclear(res);
    }
    plus->built = GV_BUILD_AREAS;
 
    /* attach isles */
    if (pg_info->topo_geo_only) {
        plus->n_isles = 0; /* reset isles */
        G_warning(
            _("To be implemented: isles not attached in Topo-Geo-only mode"));
    }
    else {
        /* read isles from 'isle_grass' table */
        sprintf(stmt,
                "SELECT isle_id,lines,area FROM \"%s\".%s ORDER BY isle_id",
                pg_info->toposchema_name, TOPO_TABLE_ISLE);
        G_debug(2, "SQL: %s", stmt);
 
        res = PQexec(pg_info->conn, stmt);
        if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
            plus->n_isles != PQntuples(res)) {
            if (res)
                PQclear(res);
            return -1;
        }
 
        dig_alloc_isles(plus, plus->n_isles);
        G_zero(plus->Isle, sizeof(struct P_isle *) *
                               (plus->n_isles + 1)); /* index starts at 1 */
        G_debug(3, "Vect_open_topo_pg(): n_isles=%d", plus->n_isles);
 
        for (i = 0; i < plus->n_isles; i++) {
            read_p_isle(plus, i + 1, (char *)PQgetvalue(res, i, 1),
                        atoi(PQgetvalue(res, i, 2)));
 
            if (plus->Spidx_new) {
                /* update spatial index */
                Vect_get_isle_points(Map, i + 1, Points);
                dig_line_box(Points, &box);
                dig_spidx_add_isle(&(Map->plus), i + 1, &box);
            }
        }
        PQclear(res);
    }
    plus->built = GV_BUILD_ATTACH_ISLES;
 
    /* attach centroids */
    if (pg_info->topo_geo_only && plus->n_areas > 0) {
        int area;
        struct P_area *Area;
        struct P_topo_c *topo;
 
        for (line = 1; line <= plus->n_lines; line++) {
            Line = plus->Line[line];
            if (Line->type != GV_CENTROID)
                continue;
 
            Vect_read_line(Map, Points, NULL, line);
            area = Vect_find_area(Map, Points->x[0], Points->y[0]);
            topo = (struct P_topo_c *)Line->topo;
            topo->area = area;
            Area = plus->Area[topo->area];
            Area->centroid = Line->offset;
        }
    }
    plus->built = GV_BUILD_CENTROIDS;
 
    /* done */
    plus->built = GV_BUILD_ALL;
 
    Vect_destroy_line_struct(Points);
    Vect_destroy_list(List);
 
    return 0;
}
 
/*!
   \brief Read nodes from DB
 
   \param Map pointer to Map_info struct
   \param geom_only read only node's geometry
 
   \return number of nodes
 */
int Vect__load_map_nodes_pg(struct Map_info *Map, int geom_only)
{
    int i, id, n_nodes;
    char stmt[DB_SQL_MAX];
    struct Plus_head *plus;
    struct Format_info_pg *pg_info;
    struct Format_info_offset *offset;
 
    PGresult *res;
 
    plus = &(Map->plus);
    pg_info = &(Map->fInfo.pg);
    offset = &(pg_info->offset);
 
    if (pg_info->topo_geo_only || geom_only)
        sprintf(stmt,
                "SELECT node_id,geom FROM \"%s\".node WHERE node_id IN "
                "(SELECT node FROM (SELECT start_node AS node FROM \"%s\".edge "
                "GROUP BY start_node UNION ALL SELECT end_node AS node FROM "
                "\"%s\".edge GROUP BY end_node) AS foo) ORDER BY node_id",
                pg_info->toposchema_name, pg_info->toposchema_name,
                pg_info->toposchema_name);
    else
        sprintf(
            stmt,
            "SELECT node.node_id,geom,lines,angles FROM \"%s\".node AS node "
            "JOIN \"%s\".%s AS node_grass ON node.node_id = node_grass.node_id "
            "ORDER BY node_id",
            pg_info->toposchema_name, pg_info->toposchema_name,
            TOPO_TABLE_NODE);
    G_debug(2, "SQL: %s", stmt);
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
        (!geom_only && PQntuples(res) != plus->n_nodes)) {
        G_warning(_("Inconsistency in topology: number of "
                    "nodes %d (should be %d)"),
                  PQntuples(res), plus->n_nodes);
        if (res)
            PQclear(res);
        return -1;
    }
 
    n_nodes = PQntuples(res);
    G_debug(3, "load_plus(): n_nodes = %d", n_nodes);
    dig_alloc_nodes(plus, n_nodes);
    offset->array = (int *)G_malloc(sizeof(int) * n_nodes);
    offset->array_num = offset->array_alloc = n_nodes;
    for (i = 0; i < n_nodes; i++) {
        G_debug(5, "node: %d", i);
        id = atoi(PQgetvalue(res, i, 0));
        read_p_node(
            plus, i + 1, /* node index starts at 1 */
            id, (const char *)PQgetvalue(res, i, 1),
            !pg_info->topo_geo_only ? (const char *)PQgetvalue(res, i, 2)
                                    : NULL,
            !pg_info->topo_geo_only ? (const char *)PQgetvalue(res, i, 3)
                                    : NULL,
            pg_info, geom_only);
        /* update offset */
        offset->array[i] = id;
    }
    PQclear(res);
 
    return n_nodes;
}
 
/*!
   \brief Read features from DB
 
   \param Map pointer to Map_info struct
 
   \return number of features
 */
int Vect__load_map_lines_pg(struct Map_info *Map)
{
    int i, id, ntuples;
 
    char stmt[DB_SQL_MAX];
 
    struct Plus_head *plus;
    struct Format_info_pg *pg_info;
    struct line_data line_data;
    struct Format_info_offset *offset;
 
    PGresult *res;
 
    plus = &(Map->plus);
    pg_info = &(Map->fInfo.pg);
    offset = &(pg_info->offset);
 
    dig_alloc_lines(plus, plus->n_lines);
    G_zero(plus->Line, sizeof(struct P_line *) *
                           (plus->n_lines + 1)); /* index starts at 1 */
 
    /* read PostGIS Topo standalone nodes (containing_face is null)
       -> points
     */
    if (pg_info->topo_geo_only)
        sprintf(stmt,
                "SELECT tt.node_id,tt.geom,ft.%s FROM \"%s\".node AS tt "
                "LEFT JOIN \"%s\".\"%s\" AS ft ON "
                "(%s).type = 1 AND (%s).id = node_id WHERE containing_face "
                "IS NULL AND node_id NOT IN "
                "(SELECT node FROM (SELECT start_node AS node FROM \"%s\".edge "
                "GROUP BY start_node UNION ALL SELECT end_node AS node FROM "
                "\"%s\".edge GROUP BY end_node) AS foo) ORDER BY node_id",
                pg_info->fid_column, pg_info->toposchema_name,
                pg_info->schema_name, pg_info->table_name,
                pg_info->topogeom_column, pg_info->topogeom_column,
                pg_info->toposchema_name, pg_info->toposchema_name);
    else
        sprintf(stmt,
                "SELECT tt.node_id,tt.geom,ft.%s "
                "FROM \"%s\".node AS tt LEFT JOIN \"%s\".\"%s\" AS ft ON "
                "(%s).type = 1 AND (%s).id = node_id WHERE node_id NOT IN "
                "(SELECT node_id FROM \"%s\".%s) AND containing_face IS NULL "
                "ORDER BY node_id",
                pg_info->fid_column, pg_info->toposchema_name,
                pg_info->schema_name, pg_info->table_name,
                pg_info->topogeom_column, pg_info->topogeom_column,
                pg_info->toposchema_name, TOPO_TABLE_NODE);
    G_debug(2, "SQL: %s", stmt);
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
        PQntuples(res) > plus->n_plines) {
        G_warning(_("Inconsistency in topology: number of "
                    "points %d (should be %d)"),
                  PQntuples(res), plus->n_plines);
        if (res)
            PQclear(res);
        return -1;
    }
 
    ntuples = PQntuples(res); /* plus->n_plines */
    G_zero(&line_data, sizeof(struct line_data));
    for (i = 0; i < ntuples; i++) {
        /* process standalone nodes (PostGIS Topo) */
        line_data.id = atoi(PQgetvalue(res, i, 0));
        line_data.wkb_geom = (char *)PQgetvalue(res, i, 1);
        line_data.fid = atoi(PQgetvalue(res, i, 2)); /* feature id */
 
        read_p_line(plus, i + 1, &line_data, pg_info->topo_geo_only,
                    &(pg_info->cache));
    }
    PQclear(res);
 
    /* read PostGIS Topo edges
       -> lines
       -> boundaries
     */
    if (pg_info->topo_geo_only)
        sprintf(stmt,
                "SELECT edge_id,start_node,end_node,left_face,right_face AS "
                "right_area,tt.geom,ft.%s "
                "FROM \"%s\".edge AS tt LEFT JOIN \"%s\".\"%s\" AS ft ON "
                "(%s).type = 2 AND "
                "(%s).id = edge_id ORDER BY edge_id",
                pg_info->fid_column, pg_info->toposchema_name,
                pg_info->schema_name, pg_info->table_name,
                pg_info->topogeom_column, pg_info->topogeom_column);
    else
        sprintf(
            stmt,
            "SELECT "
            "edge_id,start_node,end_node,left_area,right_area,tt.geom,ft.%s "
            "FROM \"%s\".edge AS tt LEFT JOIN \"%s\".\"%s\" ON "
            "edge_id = line_id LEFT JOIN \"%s\".\"%s\" AS ft ON (%s).type = 2 "
            "AND "
            "(%s).id = edge_id ORDER BY edge_id",
            pg_info->fid_column, pg_info->toposchema_name,
            pg_info->toposchema_name, TOPO_TABLE_LINE, pg_info->schema_name,
            pg_info->table_name, pg_info->topogeom_column,
            pg_info->topogeom_column);
 
    G_debug(2, "SQL: %s", stmt);
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
        PQntuples(res) > plus->n_lines) {
        G_warning(_("Inconsistency in topology: number of "
                    "lines %d (should be %d)"),
                  PQntuples(res), plus->n_lines);
        if (res)
            PQclear(res);
        return -1;
    }
 
    /* process edges (PostGIS Topo) */
    ntuples = PQntuples(res);
    for (i = 0; i < ntuples; i++) {
        line_data.id = atoi(PQgetvalue(res, i, 0));
        line_data.start_node = remap_node(offset, atoi(PQgetvalue(res, i, 1)));
        line_data.end_node = remap_node(offset, atoi(PQgetvalue(res, i, 2)));
        line_data.left_face = atoi(PQgetvalue(res, i, 3));
        line_data.right_face = atoi(PQgetvalue(res, i, 4));
        line_data.wkb_geom = (char *)PQgetvalue(res, i, 5);
        line_data.fid = atoi(PQgetvalue(res, i, 6)); /* feature id */
 
        id = plus->n_plines + i + 1; /* points already registered */
        read_p_line(plus, id, &line_data, pg_info->topo_geo_only,
                    &(pg_info->cache));
    }
    PQclear(res);
 
    /* read PostGIS Topo standalone nodes (containing_face is not null)
       -> centroids
     */
    if (pg_info->topo_geo_only)
        sprintf(
            stmt,
            "SELECT node_id,tt.geom,containing_face,ft.%s FROM "
            "\"%s\".node AS tt LEFT JOIN \"%s\".\"%s\" AS ft ON "
            "(%s).type = 3 AND (%s).id = containing_face WHERE containing_face "
            "IS NOT NULL AND node_id NOT IN "
            "(SELECT node FROM (SELECT start_node AS node FROM \"%s\".edge "
            "GROUP BY start_node UNION ALL SELECT end_node AS node FROM "
            "\"%s\".edge GROUP BY end_node) AS foo) ORDER BY node_id",
            pg_info->fid_column, pg_info->toposchema_name, pg_info->schema_name,
            pg_info->table_name, pg_info->topogeom_column,
            pg_info->topogeom_column, pg_info->toposchema_name,
            pg_info->toposchema_name);
    else
        sprintf(stmt,
                "SELECT tt.node_id,tt.geom,containing_face,ft.%s FROM "
                "\"%s\".node AS tt LEFT JOIN \"%s\".\"%s\" AS ft ON "
                "(%s).type = 3 AND (%s).id = containing_face WHERE "
                "node_id NOT IN (SELECT node_id FROM \"%s\".%s) AND "
                "containing_face "
                "IS NOT NULL ORDER BY node_id",
                pg_info->fid_column, pg_info->toposchema_name,
                pg_info->schema_name, pg_info->table_name,
                pg_info->topogeom_column, pg_info->topogeom_column,
                pg_info->toposchema_name, TOPO_TABLE_NODE);
    G_debug(2, "SQL: %s", stmt);
    res = PQexec(pg_info->conn, stmt);
    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
        PQntuples(res) != plus->n_clines) {
        G_warning(_("Inconsistency in topology: number of "
                    "centroids %d (should be %d)"),
                  PQntuples(res), plus->n_clines);
        if (res)
            PQclear(res);
        return -1;
    }
 
    G_zero(&line_data, sizeof(struct line_data));
    id = plus->n_plines + plus->n_llines + plus->n_blines + 1;
    for (i = 0; i < plus->n_clines; i++) {
        line_data.id = atoi(PQgetvalue(res, i, 0));
        line_data.wkb_geom = (char *)PQgetvalue(res, i, 1);
        line_data.left_face = atoi(PQgetvalue(res, i, 2)); /* face id */
        line_data.fid = atoi(PQgetvalue(res, i, 3));       /* feature id */
        /* area id and face id can be different */
 
        read_p_line(plus, id + i, &line_data, pg_info->topo_geo_only,
                    &(pg_info->cache));
    }
    PQclear(res);
 
    plus->built = GV_BUILD_BASE;
 
    return plus->n_lines;
}
 
/*
   \brief PostgreSQL notice processor
 
   Print out NOTICE message only on verbose level
 */
void notice_processor(void *arg UNUSED, const char *message)
{
    if (G_verbose() > G_verbose_std()) {
        fprintf(stderr, "%s", message);
    }
}
 
/*!
   \brief Scan string array
 
   Creates tokens based on string array, eg. '{1, 2, 3}' become
   [1,2,3].
 
   Allocated tokes should be freed by G_free_tokens().
 
   \param sArray string array
 
   \return tokens
 */
char **scan_array(const char *sarray)
{
    char *buf, **tokens;
    int i, len;
 
    /* remove '{}' */
    len = strlen(sarray) - 1; /* skip '}' */
    buf = (char *)G_malloc(len);
 
    for (i = 1; i < len; i++)
        buf[i - 1] = sarray[i];
    buf[len - 1] = '\0';
 
    tokens = G_tokenize(buf, ",");
    G_free(buf);
 
    return tokens;
}
 
/*!
   \brief Get node id from offset
 
   \param offset pointer to Format_info_offset struct
   \param node node to find
 
   \return node id
   \return -1 not found
 */
int remap_node(const struct Format_info_offset *offset, int node)
{
    /* probably not needed
       int i;
       for (i = node-1; i < offset->array_num; i++) {
       if (offset->array[i] == node)
       return i + 1;
       }
 
       return -1;
     */
 
    return offset->array[node - 1];
}
 
/*!
   \brief Get line id from offset
 
   \param plus pointer to Plus_head struct
   \param offset line offset
   \param type line type
 
   \return line id
   \return -1 not found
 */
int remap_line(const struct Plus_head *plus, off_t offset, int type)
{
    int i;
 
    struct P_line *Line;
 
    for (i = (int)offset; i <= plus->n_lines; i++) {
        Line = plus->Line[i];
 
        if (!Line || Line->type != type)
            continue;
 
        if ((int)Line->offset == offset)
            return i;
    }
 
    return -1;
}
#endif