[GRASS-SVN] r52017 - in grass/trunk/lib/vector: Vlib diglib

svn_grass at osgeo.org svn_grass at osgeo.org
Sat Jun 9 10:30:57 PDT 2012


Author: martinl
Date: 2012-06-09 10:30:56 -0700 (Sat, 09 Jun 2012)
New Revision: 52017

Modified:
   grass/trunk/lib/vector/Vlib/build_sfa.c
   grass/trunk/lib/vector/Vlib/open_pg.c
   grass/trunk/lib/vector/Vlib/pg_local_proto.h
   grass/trunk/lib/vector/Vlib/read_pg.c
   grass/trunk/lib/vector/diglib/spindex.c
Log:
vlib/pg: fix various issues when building areas/isles from PostGIS topology schema


Modified: grass/trunk/lib/vector/Vlib/build_sfa.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_sfa.c	2012-06-08 21:43:41 UTC (rev 52016)
+++ grass/trunk/lib/vector/Vlib/build_sfa.c	2012-06-09 17:30:56 UTC (rev 52017)
@@ -345,8 +345,8 @@
 	G_progress(iFeature + 1, 1e4);
 
 	/* cache feature (lines) */
-	if (SF_NONE == cache_feature(wkb_data,
-				     FALSE, &(pg_info->cache), &fparts)) {
+	if (SF_NONE == cache_feature(wkb_data, FALSE, FALSE,
+                                     &(pg_info->cache), &fparts)) {
 	    G_warning(_("Feature %d without geometry skipped"),
 		      iFeature + 1);
 	    continue;
@@ -681,9 +681,7 @@
 */
 int Vect__build_sfa(struct Map_info *Map, int build)
 {
-    int line;
     struct Plus_head *plus;
-    struct P_line *Line;
     
     plus = &(Map->plus);
     

Modified: grass/trunk/lib/vector/Vlib/open_pg.c
===================================================================
--- grass/trunk/lib/vector/Vlib/open_pg.c	2012-06-08 21:43:41 UTC (rev 52016)
+++ grass/trunk/lib/vector/Vlib/open_pg.c	2012-06-09 17:30:56 UTC (rev 52017)
@@ -29,6 +29,7 @@
     int end_node;
     int left_face;
     int right_face;
+    char *wkb_geom;
 };
 
 static char *get_key_column(struct Format_info_pg *);
@@ -40,8 +41,11 @@
 static int check_topo(struct Format_info_pg *, struct Plus_head *);
 static int parse_bbox(const char *, struct bound_box *);
 static int num_of_records(const struct Format_info_pg *, const char *);
-static struct P_node *read_p_node(struct Plus_head *, int, int, struct Format_info_pg *);
-static struct P_line *read_p_line(struct Plus_head *, int, const struct edge_data *);
+static struct P_node *read_p_node(struct Plus_head *, int, int,
+                                  const char *, struct Format_info_pg *);
+static struct P_line *read_p_line(struct Plus_head *, int,
+                                  const struct edge_data *,
+                                  struct Format_info_cache *);
 static int load_plus_head(struct Format_info_pg *, struct Plus_head *);
 static void notice_processor(void *, const char *);
 #endif
@@ -955,18 +959,21 @@
   \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 pg_info pointer to Format_info_pg sttucture
 
   \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,
+struct P_node *read_p_node(struct Plus_head *plus, int n,
+                           int id, const char *wkb_data,
                            struct Format_info_pg *pg_info)
 {
     int i, cnt;
     char stmt[DB_SQL_MAX];
     
     struct P_node *node;
+    struct line_pnts *points;
     
     PGresult *res;
     
@@ -1021,33 +1028,23 @@
     }
     PQclear(res);
     
-    /* ???
-    if (plus->with_z)
-    if (0 >= dig__fread_port_P(&n_edges, 1, fp))
-    */
+    /* get node coordinates */
+    if (SF_POINT != cache_feature(wkb_data, FALSE, FALSE,
+                                  &(pg_info->cache), NULL))
+        G_warning(_("Node %d: unexpected feature type %d"),
+                  n, pg_info->cache.sf_type);
     
-    /* get node coordinates */
-    sprintf(stmt,
-            "SELECT ST_X(geom),ST_Y(geom),ST_Z(geom) FROM \"%s\".node "
-            "WHERE node_id = %d",
-            pg_info->toposchema_name, id);
-    G_debug(2, "SQL: %s", stmt);
-    res = PQexec(pg_info->conn, stmt);
-    if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
-        PQntuples(res) != 1) {
-        G_warning(_("Unable to read node %d"), id);
-        if (res)
-            PQclear(res);
-        return NULL;
-    }
-    node->x = atof(PQgetvalue(res, 0, 0));
-    node->y = atof(PQgetvalue(res, 0, 1));
+    points = pg_info->cache.lines[0];
+    node->x = points->x[0];
+    node->y = points->y[0];
     if (plus->with_z)
-        node->z = atof(PQgetvalue(res, 0, 2));
+        node->z = points->z[0];
     else
-        node->z = 0;
-    PQclear(res);
-
+        node->z = 0.0;
+    
+    /* update spatial index */
+    dig_spidx_add_node(plus, n, node->x, node->y, node->z);
+    
     plus->Node[n] = node;
     
     return node;
@@ -1072,11 +1069,15 @@
   \return NULL on error
 */
 struct P_line *read_p_line(struct Plus_head *plus, int n,
-                           const struct edge_data *data)
+                           const struct edge_data *data,
+                           struct Format_info_cache *cache)
 {
-    int tp;
+    int tp, itype;
     struct P_line *line;
     
+    struct line_pnts *points;
+    struct bound_box box;
+    
     if (data->start_node == 0 && data->end_node == 0) {
         if (data->left_face == 0)
             tp = GV_POINT;
@@ -1131,11 +1132,22 @@
         else if (line->type == GV_CENTROID) {
             struct P_topo_c *topo = (struct P_topo_c *)line->topo;
             
-            topo->area  = data->left_face;
+            topo->area = data->left_face;
         }
         /* TODO: faces | kernels */
     }
 
+    /* update spatial index */
+    cache_feature(data->wkb_geom, FALSE, FALSE, cache, NULL);
+    itype = cache->lines_types[0];
+    if ((line->type & GV_POINTS && itype != GV_POINT) ||
+        (line->type & GV_LINES  && itype != GV_LINE))
+        G_warning(_("Line %d: unexpected feature type"), n);
+    
+    points = cache->lines[0];
+    dig_line_box(points, &box);
+    dig_spidx_add_line(plus, n, &box);
+    
     plus->Line[n] = line;
     
     return line;
@@ -1279,7 +1291,7 @@
        note: standalone nodes are ignored
     */
     sprintf(stmt,
-            "SELECT node_id FROM \"%s\".node WHERE node_id IN "
+            "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)",
@@ -1297,9 +1309,10 @@
 
     G_debug(3, "load_plus(): n_nodes = %d", plus->n_nodes);
     dig_alloc_nodes(plus, plus->n_nodes);
-    for (i = 1; i <= plus->n_nodes; i++) {
-        id = atoi(PQgetvalue(res, i - 1, 0));
-        read_p_node(plus, i, id, pg_info);
+    for (i = 0; i < plus->n_nodes; 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);
     }
     PQclear(res);
 
@@ -1314,7 +1327,7 @@
        -> points
     */
     sprintf(stmt,
-            "SELECT node_id FROM \"%s\".node WHERE node_id NOT IN "
+            "SELECT node_id,geom FROM \"%s\".node WHERE 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)",
@@ -1335,7 +1348,8 @@
     for (i = 0; i < ntuples; i++) {
         /* process standalone nodes (PostGIS Topo) */
         line_data.id = atoi(PQgetvalue(res, i, 0));
-        read_p_line(plus, i + 1, &line_data);
+        line_data.wkb_geom = (char *) PQgetvalue(res, i, 1);
+        read_p_line(plus, i + 1, &line_data, &(pg_info->cache));
     }
     PQclear(res);
     
@@ -1344,7 +1358,7 @@
        -> boundaries
     */
     sprintf(stmt,
-            "SELECT edge_id,start_node,end_node,left_face,right_face "
+            "SELECT edge_id,start_node,end_node,left_face,right_face,geom "
             "FROM \"%s\".edge",
             pg_info->toposchema_name);
     G_debug(2, "SQL: %s", stmt);
@@ -1365,9 +1379,10 @@
         line_data.end_node   = 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);
         
         id = plus->n_plines + i + 1; /* points already registered */
-        line = read_p_line(plus, id, &line_data);
+        line = read_p_line(plus, id, &line_data, &(pg_info->cache));
         if (line_data.left_face != 0 || line_data.right_face != 0) {
             /* boundary detected -> build area/isle on left and right*/
             
@@ -1391,6 +1406,7 @@
             if (line_data.right_face == 0)
                 line_data.right_face = -1;
             
+            /* check topo - left / right areas */
             topo = (struct P_topo_b *)line->topo;
             if (topo->left != line_data.left_face)
                 G_warning(_("Left area detected as %d (should be %d"),
@@ -1403,15 +1419,36 @@
     PQclear(res);
     
     /* attach centroids */
-    G_zero(&line_data, sizeof(struct edge_data));
-    for (i = 1; i <= plus->n_areas; i++) {
-        area = plus->Area[i];
+    if (plus->n_areas > 0) {
+        sprintf(stmt,
+                "SELECT ST_PointOnSurface(geom) AS geom FROM "
+                "ST_GetFaceGeometry('%s',"
+                "(SELECT face_id FROM \"%s\".face WHERE face_id > 0)) "
+                "AS geom",
+                pg_info->toposchema_name, pg_info->toposchema_name);
+        G_debug(2, "SQL: %s", stmt);
+        res = PQexec(pg_info->conn, stmt);
+        if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
+            PQntuples(res) > plus->n_areas) {
+            G_warning(_("Unable to attach centroids"));
+            if (res)
+                PQclear(res);
+            return -1;
+        }
         
-        line_data.id = plus->n_lines - plus->n_clines + i;
-        line_data.left_face = i;
-        read_p_line(plus, line_data.id, &line_data);
-        area->centroid = line_data.id;
+        G_zero(&line_data, sizeof(struct edge_data));
+        for (i = 1; i <= plus->n_areas; i++) {
+            area = plus->Area[i];
+            
+            id = plus->n_lines - plus->n_clines + i;
+            line_data.id = line_data.left_face = i;
+            line_data.wkb_geom = (char *)PQgetvalue(res, 0, 0);
+            
+            read_p_line(plus, id, &line_data, &(pg_info->cache));
+            area->centroid = line_data.id;
+        }
     }
+
     return 0;
 }
 

Modified: grass/trunk/lib/vector/Vlib/pg_local_proto.h
===================================================================
--- grass/trunk/lib/vector/Vlib/pg_local_proto.h	2012-06-08 21:43:41 UTC (rev 52016)
+++ grass/trunk/lib/vector/Vlib/pg_local_proto.h	2012-06-09 17:30:56 UTC (rev 52017)
@@ -49,7 +49,7 @@
 
 /* functions used in *_pg.c files */
 int execute(PGconn *, const char *);
-SF_FeatureType cache_feature(const char *, int,
+SF_FeatureType cache_feature(const char *, int, int,
 			     struct Format_info_cache *,
 			     struct feat_parts *);
 int set_initial_query(struct Format_info_pg *, int);

Modified: grass/trunk/lib/vector/Vlib/read_pg.c
===================================================================
--- grass/trunk/lib/vector/Vlib/read_pg.c	2012-06-08 21:43:41 UTC (rev 52016)
+++ grass/trunk/lib/vector/Vlib/read_pg.c	2012-06-09 17:30:56 UTC (rev 52017)
@@ -467,7 +467,8 @@
 {
     char *data;
     char stmt[DB_SQL_MAX];
-
+    int left_face, right_face;
+    
     if (!topotable_name && !pg_info->geom_column) {
         G_warning(_("No geometry column defined"));
         return -1;
@@ -500,7 +501,8 @@
         else {
             /* topological access */
             sprintf(stmt,
-                    "DECLARE %s_%s%p CURSOR FOR SELECT geom FROM \"%s\".\"%s\" "
+                    "DECLARE %s_%s%p CURSOR FOR SELECT geom,left_face,right_face "
+                    " FROM \"%s\".\"%s\" "
                     "WHERE %s_id = %d", pg_info->schema_name, pg_info->table_name,
                     pg_info->conn, pg_info->toposchema_name,
                     topotable_name, topotable_name, fid);
@@ -557,9 +559,21 @@
         return -2;
     }
     data = (char *)PQgetvalue(pg_info->res, pg_info->next_line, 0);
-
+    left_face = right_face = 0;
+    if (pg_info->toposchema_name) { /* -> PostGIS topology access */
+        if (fid < 0) { /* sequential access */
+            left_face  = atoi(PQgetvalue(pg_info->res, pg_info->next_line, 2));
+            right_face = atoi(PQgetvalue(pg_info->res, pg_info->next_line, 3));
+        }
+        else {         /* random access */
+            left_face  = atoi(PQgetvalue(pg_info->res, pg_info->next_line, 1));
+            right_face = atoi(PQgetvalue(pg_info->res, pg_info->next_line, 2));
+        }
+    }
     pg_info->cache.sf_type =
-        cache_feature(data, FALSE, &(pg_info->cache), NULL);
+        cache_feature(data, FALSE,
+                      left_face != 0 || right_face != 0 ? TRUE : FALSE,
+                      &(pg_info->cache), NULL);
     if (fid < 0) {
         pg_info->cache.fid =
             atoi(PQgetvalue(pg_info->res, pg_info->next_line, 1));
@@ -633,6 +647,7 @@
 
    \param data HEX data
    \param skip_polygon skip polygons (level 1)
+   \param force_boundary force GV_BOUNDARY feature type (for PostGIS topology)
    \param[out] cache lines cache
    \param[out] fparts used for building pseudo-topology (or NULL)
 
@@ -640,7 +655,8 @@
    \return SF_UNKNOWN on error
  */
 SF_FeatureType cache_feature(const char *data, int skip_polygon,
-                             struct Format_info_cache * cache,
+                             int force_boundary,
+                             struct Format_info_cache *cache,
                              struct feat_parts * fparts)
 {
     int ret, byte_order, nbytes, is3D;
@@ -738,7 +754,10 @@
     }
     else if (ftype == SF_LINESTRING) {
         cache->lines_num = 1;
-        cache->lines_types[0] = GV_LINE;
+        if (force_boundary)
+            cache->lines_types[0] = GV_BOUNDARY;
+        else
+            cache->lines_types[0] = GV_LINE;
         ret = linestring_from_wkb(wkb_data, nbytes, byte_order,
                                   is3D, cache->lines[0], FALSE);
         add_fpart(fparts, ftype, 0, 1);
@@ -1119,12 +1138,14 @@
         sprintf(stmt,
                 "DECLARE %s_%s%p CURSOR FOR "
                 "SELECT geom,row_number() OVER "
-                "(ORDER BY ST_GeometryType(geom) DESC) AS fid FROM ("
-                "SELECT geom FROM \"%s\".node WHERE node_id NOT IN "
+                "(ORDER BY ST_GeometryType(geom) DESC) AS fid,"
+                "left_face,right_face FROM ("
+                "SELECT geom,0 AS left_face, 0 AS right_face FROM "
+                "\"%s\".node WHERE 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 SELECT geom FROM \"%s\".edge) AS foo",
+                "UNION ALL SELECT geom,left_face,right_face FROM \"%s\".edge) AS foo",
                 pg_info->schema_name, pg_info->table_name, pg_info->conn,
                 pg_info->toposchema_name, pg_info->toposchema_name,
                 pg_info->toposchema_name, pg_info->toposchema_name); 
@@ -1257,7 +1278,7 @@
     
     sprintf(stmt,
             "SELECT ST_PointOnSurface(geom) AS geom FROM "
-            "ST_GetFaceGeometry('%s', %d) as geom",
+            "ST_GetFaceGeometry('%s', %d) AS geom",
             pg_info->toposchema_name, centroid);
     res = PQexec(pg_info->conn, stmt);
     if (!res || PQresultStatus(res) != PGRES_TUPLES_OK ||
@@ -1272,7 +1293,7 @@
     data = (char *)PQgetvalue(res, 0, 0);
     PQclear(res);
     
-    if (GV_POINT != cache_feature(data, FALSE, &(pg_info->cache), NULL))
+    if (GV_POINT != cache_feature(data, FALSE, FALSE, &(pg_info->cache), NULL))
         return -1;
     
     Vect_append_points(line_p, pg_info->cache.lines[0], GV_FORWARD);

Modified: grass/trunk/lib/vector/diglib/spindex.c
===================================================================
--- grass/trunk/lib/vector/diglib/spindex.c	2012-06-08 21:43:41 UTC (rev 52016)
+++ grass/trunk/lib/vector/diglib/spindex.c	2012-06-09 17:30:56 UTC (rev 52017)
@@ -287,7 +287,7 @@
 
    \return 0
  */
-int dig_spidx_add_line(struct Plus_head *Plus, int line, struct bound_box * box)
+int dig_spidx_add_line(struct Plus_head *Plus, int line, const struct bound_box * box)
 {
     static struct RTree_Rect rect;
     static int rect_init = 0;
@@ -320,7 +320,7 @@
 
    \return 0
  */
-int dig_spidx_add_area(struct Plus_head *Plus, int area, struct bound_box * box)
+int dig_spidx_add_area(struct Plus_head *Plus, int area, const struct bound_box * box)
 {
     static struct RTree_Rect rect;
     static int rect_init = 0;
@@ -354,7 +354,7 @@
    \return 0
  */
 
-int dig_spidx_add_isle(struct Plus_head *Plus, int isle, struct bound_box * box)
+int dig_spidx_add_isle(struct Plus_head *Plus, int isle, const struct bound_box * box)
 {
     static struct RTree_Rect rect;
     static int rect_init = 0;
@@ -433,7 +433,6 @@
 int dig_spidx_del_line(struct Plus_head *Plus, int line, double x, double y, double z)
 {
     int ret;
-    struct P_line *Line;
     static struct RTree_Rect rect;
     static int rect_init = 0;
 
@@ -445,8 +444,6 @@
 
     G_debug(3, "dig_spidx_del_line(): line = %d", line);
 
-    Line = Plus->Line[line];
-
     rect.boundary[0] = x;
     rect.boundary[1] = y;
     rect.boundary[2] = z;



More information about the grass-commit mailing list