vlib/pg: implement creating simple features from topogeometry data

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@58363 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/Vlib/build_pg.c

@@ -35,6 +35,7 @@ static int write_lines(const struct Plus_head *, const struct Format_info_pg *);
 static int write_areas(const struct Plus_head *, const struct Format_info_pg *);
 static int write_isles(const struct Plus_head *, const struct Format_info_pg *);
 static void build_stmt_id(const void *, int, int, const struct Plus_head *, char **, size_t *);
+static int create_simple_feature_from_topo(struct Map_info *);
 #endif
 
 /*!
@@ -136,6 +137,7 @@ int build_topo(struct Map_info *Map, int build)
     int area, nareas, isle, nisles;
     int face[2];
     char stmt[DB_SQL_MAX];
+    char *def_file;
     
     struct Plus_head *plus;
     struct Format_info_pg *pg_info;
@@ -395,6 +397,31 @@ int build_topo(struct Map_info *Map, int build)
     if (Vect__execute_pg(pg_info->conn, "COMMIT") == -1)
         return 0;
 
+    /* check if we want to create simple features from topogeometry
+       data */
+    def_file = getenv("GRASS_VECTOR_PGFILE");
+    
+    if (G_find_file2("", def_file ? def_file : "PG", G_mapset())) {
+        FILE *fp;
+        const char *p;
+        
+        struct Key_Value *key_val;
+        fp = G_fopen_old("", def_file ? def_file : "PG", G_mapset());
+        if (!fp) {
+            G_fatal_error(_("Unable to open PG file"));
+        }
+        key_val = G_fread_key_value(fp);
+        fclose(fp);
+        
+        /* build simple features from topogeometry data */
+        p = G_find_key_value("simple_feature", key_val);
+        if (p && G_strcasecmp(p, "yes") == 0)
+            if (create_simple_feature_from_topo(Map) != 0)
+                return 0;
+
+        G_free_key_value(key_val);
+    }
+    
     return 1;
 }
 
@@ -948,4 +975,48 @@ int Vect__clean_grass_db_topo(struct Format_info_pg *pg_info)
     return 0;
 }
 
+/*!
+  \brief Create simple features geometry from topogeometry data
+
+  \param Map pointer to Map_info struct 
+
+  \return 0 on success
+  \return -1 on error
+*/
+int create_simple_feature_from_topo(struct Map_info *Map)
+{
+    char stmt[DB_SQL_MAX];
+
+    struct Format_info_pg *pg_info;
+
+    pg_info = &(Map->fInfo.pg);
+
+    G_debug(1, "build_simple_feature_from_topo(): %d", pg_info->feature_type);
+    
+    G_message(_("Create simple features topology from topogeometry data..."));
+    Vect__execute_pg(pg_info->conn, "BEGIN");
+    if (pg_info->feature_type == SF_POINT ||
+        pg_info->feature_type == SF_LINESTRING) {
+        sprintf(stmt, "UPDATE \"%s\".\"%s\" SET %s = (SELECT geom FROM \"%s\".node "
+                "WHERE node_id = (%s).id)", pg_info->schema_name, pg_info->table_name,
+                pg_info->geom_column, pg_info->toposchema_name, pg_info->topogeom_column);
+
+        if(Vect__execute_pg(pg_info->conn, stmt) == -1) {
+            Vect__execute_pg(pg_info->conn, "ROLLBACK");
+            return -1;
+
+        }
+    }
+    else if (pg_info->feature_type == SF_POLYGON) {
+        Vect__copy_areas(Map, 1, Map);
+    }
+    else {
+        G_warning(_("Unable to build simple features from topogeometry data. "
+                    "Unsupported type %d."), pg_info->feature_type);
+    }
+ 
+    Vect__execute_pg(pg_info->conn, "COMMIT");
+
+    return 0;
+}
 #endif

lib/vector/Vlib/build_sfa.c

@@ -329,7 +329,7 @@ void build_pg(struct Map_info *Map, int build)
     G_zero(&fparts, sizeof(struct feat_parts));
     
     /* get all features */
-    if (Vect__open_cursor_next_line_pg(pg_info, TRUE) != 0)
+    if (Vect__open_cursor_next_line_pg(pg_info, TRUE, Map->plus.built) != 0)
         return;
     
     /* scan records */

lib/vector/Vlib/copy.c

@@ -43,7 +43,6 @@ static int copy_nodes(const struct Map_info *, struct Map_info *);
 static int copy_line_nodes(const struct Map_info *, int, int, struct line_pnts *,
                            struct Map_info *);
 static int is_isle(const struct Map_info *, int);
-static int copy_areas(const struct Map_info *, int, struct Map_info *);
 
 /*!
    \brief Copy all alive vector features from input vector map to
@@ -125,7 +124,7 @@ int Vect_copy_map_lines_field(struct Map_info *In, int field,
 
         if (topo == TOPO_NONE) {
             /* copy areas - external formats and simple features access only */
-            ret += copy_areas(In, field, Out);
+            ret += Vect__copy_areas(In, field, Out);
         }
     }
     else {
@@ -460,7 +459,7 @@ int is_isle(const struct Map_info *Map, int area)
   \return 0 on success
   \return 1 on error
 */
-int copy_areas(const struct Map_info *In, int field, struct Map_info *Out)
+int Vect__copy_areas(const struct Map_info *In, int field, struct Map_info *Out)
 {
     int i, area, nareas, cat, isle, nisles, nparts_alloc, nskipped;
     struct line_pnts **Points;
@@ -526,11 +525,22 @@ int copy_areas(const struct Map_info *In, int field, struct Map_info *Out)
             Vect_get_isle_points(In, isle, Points[i + 1]);
         }
         
-        if (0 > V2__write_area_sfa(Out, (const struct line_pnts **) Points,
-                                   nisles + 1, Cats)) {
-            G_warning(_("Writing area %d failed"), area);
-            return -1;
+        if (In != Out) {
+            if (0 > V2__write_area_sfa(Out, (const struct line_pnts **) Points,
+                                       nisles + 1, Cats)) {
+                G_warning(_("Writing area %d failed"), area);
+                return -1;
+            }
+        }
+#ifdef HAVE_POSTGRES
+        else { /* building simple features geometry from topogeometry data */
+            if (0 > V2__update_area_pg(Out, (const struct line_pnts **) Points,
+                                        nisles + 1, cat)) {
+                G_warning(_("Writing area %d failed"), area);
+                return -1;
+            }
         }
+#endif
     }
 
     if (nskipped > 0)

lib/vector/Vlib/local_proto.h

@@ -18,6 +18,9 @@ int Vect__get_area_points_nat(const struct Map_info *, const plus_t *, int, stru
 void Vect__free_cache(struct Format_info_cache *);
 void Vect__free_offset(struct Format_info_offset *);
 
+/* copy.c */
+int Vect__copy_areas(const struct Map_info *, int, struct Map_info *);
+
 /* map.c */
 int Vect__delete(const char *, int);
 

lib/vector/Vlib/pg_local_proto.h

@@ -68,7 +68,7 @@ int Vect__clean_grass_db_topo(struct Format_info_pg *);
 SF_FeatureType Vect__cache_feature_pg(const char *, int, int,
                                       struct Format_info_cache *,
                                       struct feat_parts *);
-int Vect__open_cursor_next_line_pg(struct Format_info_pg *, int);
+int Vect__open_cursor_next_line_pg(struct Format_info_pg *, int, int);
 int Vect__open_cursor_line_pg(struct Format_info_pg *, int, int);
 int Vect__close_cursor_pg(struct Format_info_pg *);
 int Vect__select_line_pg(struct Format_info_pg *, int, int);
@@ -80,6 +80,9 @@ void Vect__reallocate_cache(struct Format_info_cache *, int, int);
 off_t V2__write_node_pg(struct Map_info *, const struct line_pnts *);
 off_t V2__write_area_pg(struct Map_info *, const struct line_pnts **, int,
                         const struct line_cats *);
+int V2__update_area_pg(struct Map_info *, const struct line_pnts **, int,
+                       int);
+
 int Vect__insert_face_pg(struct Map_info *, int);
 
 /* open_pg.c */

lib/vector/Vlib/read_pg.c

@@ -41,7 +41,7 @@ static unsigned int wkb_data_length;
 
 static int read_next_line_pg(struct Map_info *,
                              struct line_pnts *, struct line_cats *, int);
-SF_FeatureType get_feature(struct Format_info_pg *, int, int);
+SF_FeatureType get_feature(struct Map_info *, int, int);
 static unsigned char *hex_to_wkb(const char *, int *);
 static int point_from_wkb(const unsigned char *, int, int, int,
                           struct line_pnts *);
@@ -266,7 +266,7 @@ int V1_read_line_pg(struct Map_info *Map,
 
         G_debug(3, "read (%s) feature (fid = %ld) to cache",
                 pg_info->table_name, fid);
-        get_feature(pg_info, fid, -1);
+        get_feature(Map, fid, -1);
 
         if (pg_info->cache.sf_type == SF_NONE) {
             G_warning(_("Feature %ld without geometry skipped"), fid);
@@ -369,7 +369,7 @@ int V2_read_line_pg(struct Map_info *Map, struct line_pnts *line_p,
                       line, pg_info->cache.lines_types[cache_idx], Line->type);
     }
     else {
-        get_feature(pg_info, fid, Line->type);
+        get_feature(Map, fid, Line->type);
         cache_idx = 0;
     }
 
@@ -465,7 +465,7 @@ int read_next_line_pg(struct Map_info *Map,
             }
             
             /* cache feature -> line_p & line_c */
-            sf_type = get_feature(pg_info, -1, -1);
+            sf_type = get_feature(Map, -1, -1);
             
             if (sf_type == SF_NONE) {
                 G_warning(_("Feature %ld without geometry skipped"), pg_info->cache.fid);
@@ -554,26 +554,30 @@ int read_next_line_pg(struct Map_info *Map,
 
    Geometry is stored in lines cache.
 
-   \param[in,out] pg_info pointer to Format_info_pg struct
+   \param[in,out] Map pointer to Map_info struct
    \param fid feature id to be read (-1 for next)
    \param type feature type (GV_POINT, GV_LINE, ...) - use only for topological access
    
    \return simple feature type (SF_POINT, SF_LINESTRING, ...)
    \return -1 on error
  */
-SF_FeatureType get_feature(struct Format_info_pg *pg_info, int fid, int type)
+SF_FeatureType get_feature(struct Map_info *Map, int fid, int type)
 {
     int seq_type;
     int force_type; /* force type (GV_BOUNDARY or GV_CENTROID) for topo access only */
     char *data;
     
+    struct Format_info_pg *pg_info;
+
+    pg_info = &(Map->fInfo.pg);
+    
     if (!pg_info->geom_column && !pg_info->topogeom_column) {
         G_warning(_("No geometry or topo geometry column defined"));
         return -1;
     }
     if (fid < 1) { /* sequantial access */
         if (pg_info->cursor_name == NULL &&
-            Vect__open_cursor_next_line_pg(pg_info, FALSE) != 0)
+            Vect__open_cursor_next_line_pg(pg_info, FALSE, Map->plus.built) != 0)
         return -1;
     }
     else {         /* random access */
@@ -670,15 +674,17 @@ SF_FeatureType get_feature(struct Format_info_pg *pg_info, int fid, int type)
     
     /* cache also categories (only for PostGIS Topology) */
     if (pg_info->toposchema_name) {
-        int col_idx;
+        int cat, col_idx;
 
-        col_idx = type & GV_POINTS ? 2 : 3;
+        col_idx = fid < 0 ? 3 : 2; /* TODO: dermine col_idx for random access */
             
         if (!PQgetisnull(pg_info->res, pg_info->next_line, col_idx))
-            pg_info->cache.lines_cats[pg_info->cache.lines_next] =
-                atoi(PQgetvalue(pg_info->res, pg_info->next_line, col_idx)); 
+            cat = atoi(PQgetvalue(pg_info->res, pg_info->next_line, col_idx)); 
         else
-            pg_info->cache.lines_cats[pg_info->cache.lines_next] = -1; /* no cat */
+            cat = -1; /* no cat */
+        pg_info->cache.lines_cats[pg_info->cache.lines_next] = cat;
+        G_debug(3, "line=%d, type=%d -> cat=%d", pg_info->cache.lines_next+1,
+                pg_info->cache.lines_types[pg_info->cache.lines_next], cat);
     }
 
     /* set feature id */
@@ -1223,7 +1229,7 @@ int error_corrupted_data(const char *msg)
   \return 0 on success
   \return -1 on failure
 */
-int Vect__open_cursor_next_line_pg(struct Format_info_pg *pg_info, int fetch_all)
+int Vect__open_cursor_next_line_pg(struct Format_info_pg *pg_info, int fetch_all, int built_level)
 {
     char stmt[DB_SQL_MAX];
     
@@ -1243,7 +1249,7 @@ int Vect__open_cursor_next_line_pg(struct Format_info_pg *pg_info, int fetch_all
                 pg_info->table_name, pg_info->fid_column);
     }
     else {
-        /* topology access (geom,fid,type) */
+        /* topology access (geom,id,fid,type) */
         /* TODO: optimize SQL statement (for points/centroids) */
         sprintf(stmt,
                 "DECLARE %s CURSOR FOR "
@@ -1254,7 +1260,7 @@ int Vect__open_cursor_next_line_pg(struct Format_info_pg *pg_info, int fetch_all
                 "(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) UNION ALL "
                 "SELECT tt.node_id AS id,tt.geom, %d AS type, ft.%s AS fid FROM \"%s\".node AS tt "
-                "LEFT JOIN \"%s\".\"%s\" AS ft ON (%s).type = 3 AND (%s).id = containing_face "
+                "LEFT JOIN \"%s\".\"%s\" AS ft ON (%s).type = 3 AND (%s).id = %s "
                 "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) UNION ALL "
@@ -1268,7 +1274,9 @@ int Vect__open_cursor_next_line_pg(struct Format_info_pg *pg_info, int fetch_all
                 GV_POINT, 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,
                 GV_CENTROID, 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,
+                pg_info->topogeom_column, pg_info->topogeom_column, 
+                built_level >= GV_BUILD_CENTROIDS ? "containing_face" : "node_id",
+                pg_info->toposchema_name, pg_info->toposchema_name,
                 GV_LINE, pg_info->fid_column, pg_info->toposchema_name, pg_info->schema_name, pg_info->table_name,
                 pg_info->topogeom_column, pg_info->topogeom_column,
                 GV_BOUNDARY, pg_info->fid_column, pg_info->toposchema_name, pg_info->schema_name, pg_info->table_name,

lib/vector/Vlib/write_pg.c

@@ -516,6 +516,61 @@ off_t V2__write_area_pg(struct Map_info *Map,
 }
 
 /*!
+   \brief Updates simple features geometry from GRASS-like topo
+
+   \param Map pointer to Map_info structure
+   \param points feature geometry (exterior + interior rings)
+   \param nparts number of parts including exterior ring
+   \param cat area category
+   
+   \return 0 on success
+   \return -1 on error
+*/
+int V2__update_area_pg(struct Map_info *Map, 
+                         const struct line_pnts **points, int nparts,
+                         int cat)
+{
+    int part, npoints;
+    char *stmt, *geom_data;
+
+    struct Format_info_pg *pg_info;
+
+    pg_info = &(Map->fInfo.pg);
+
+    for (part = 0; part < nparts; part++) { 
+        npoints = points[part]->n_points - 1;
+        if (points[part]->x[0] != points[part]->x[npoints] ||
+            points[part]->y[0] != points[part]->y[npoints] ||
+            points[part]->z[0] != points[part]->z[npoints]) {
+            G_warning(_("Boundary is not closed. Skipping."));
+            return -1;
+        }
+    }
+    
+    geom_data = line_to_wkb(pg_info, points, nparts, GV_AREA, Vect_is_3d(Map));
+    if (!geom_data)
+        return -1;
+    
+    stmt = NULL;
+    G_asprintf(&stmt, "UPDATE \"%s\".\"%s\" SET %s = '%s'::GEOMETRY WHERE %s = %d",
+            pg_info->schema_name, pg_info->table_name, pg_info->geom_column,
+            geom_data, pg_info->fid_column,
+            cat);
+    if (Vect__execute_pg(pg_info->conn, stmt) == -1) {
+        /* rollback transaction */
+        Vect__execute_pg(pg_info->conn, "ROLLBACK"); 
+        G_free(geom_data);
+        G_free(stmt);
+        return -1;
+    }
+
+    G_free(geom_data);
+    G_free(stmt);
+
+    return 0;
+}
+
+/*!
   \brief Create new feature table
 
   \param pg_info pointer to Format_info_pg
@@ -559,6 +614,8 @@ int create_table(struct Format_info_pg *pg_info, const struct field_info *Fi)
             p = G_find_key_value("primary_key", key_val);
             if (p && G_strcasecmp(p, "no") == 0)
                 primary_key = FALSE;
+
+            G_free_key_value(key_val);
         }
     }
     
@@ -824,6 +881,14 @@ int create_topo_schema(struct Format_info_pg *pg_info, int with_z)
         if (p && G_strcasecmp(p, "yes") == 0)
             pg_info->topo_geo_only = TRUE;
         G_debug(1, "PG: topo_geo_only :%d", pg_info->topo_geo_only);
+
+        /* build simple features from topogeometry data */
+        p = G_find_key_value("simple_feature", key_val);
+        if (p && G_strcasecmp(p, "yes") == 0)
+            pg_info->topo_geo_only = TRUE;
+        G_debug(1, "PG: topo_geo_only :%d", pg_info->topo_geo_only);
+
+        G_free_key_value(key_val);
     }
 
     /* begin transaction (create topo schema) */
@@ -1689,8 +1754,8 @@ char *line_to_wkb(struct Format_info_pg *pg_info,
         wkb_data = point_to_wkb(byte_order, points[0], with_z, &nbytes);
     else if (type == GV_LINE)
         wkb_data = linestring_to_wkb(byte_order, points[0], with_z, &nbytes);
-    else if (type & (GV_BOUNDARY | GV_FACE)) {
-        if (!pg_info->toposchema_name) {
+    else if (type & (GV_BOUNDARY | GV_FACE | GV_AREA)) {
+        if (!pg_info->toposchema_name || type == GV_AREA) {
             /* PostGIS simple feature access */
             wkb_data = polygon_to_wkb(byte_order, points, nparts,
                                       with_z, &nbytes);
@@ -1839,6 +1904,7 @@ int write_feature(struct Map_info *Map, int line, int type,
         return -1;
     }
     G_free(geom_data);
+    G_free(stmt);
     
     return pg_info->toposchema_name ? topo_id : 0;
 }