build_sfa.c

Go to the documentation of this file.

/*!
   \file lib/vector/Vlib/build_sfa.c
 
   \brief Vector library - Building pseudo-topology for simple feature access
 
   Higher level functions for reading/writing/manipulating vectors.
 
   Line offset is
   - centroids   : FID
   - other types : index of the first record (which is FID) in offset array.
 
   (C) 2001-2012 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 Radim Blazek
   \author Piero Cavalieri
   \author Various updates for GRASS 7 by Martin Landa <landa.martin gmail.com>
 */
 
#include <stdlib.h>
 
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/glocale.h>
 
#include <ogr_api.h>
 
/*!
   \brief This structure keeps info about geometry parts above current
   geometry, path to current geometry in the feature. First 'part' number
   however is feature id */
struct geom_parts {
    int *part;
    int a_parts;
    int n_parts;
};
 
static void init_parts(struct geom_parts *);
static void reset_parts(struct geom_parts *);
static void free_parts(struct geom_parts *);
static void add_part(struct geom_parts *, int);
static void del_part(struct geom_parts *);
static void add_parts_to_offset(struct Format_info_offset *,
                                struct geom_parts *);
static int add_line(struct Plus_head *, struct Format_info_offset *, int,
                    struct line_pnts *, int, struct geom_parts *);
 
#ifdef HAVE_POSTGRES
#include "pg_local_proto.h"
 
static int add_geometry_pg(struct Plus_head *, struct Format_info_pg *,
                           struct feat_parts *, int, int, int,
                           struct geom_parts *);
static void build_pg(struct Map_info *, int);
#endif
 
static int add_geometry_ogr(struct Plus_head *, struct Format_info_ogr *,
                            OGRGeometryH, int, int, struct geom_parts *);
 
static void build_ogr(struct Map_info *, int);
 
/*!
   \brief Init parts
 */
void init_parts(struct geom_parts *parts)
{
    G_zero(parts, sizeof(struct geom_parts));
}
 
/*!
   \brief Reset parts
 */
void reset_parts(struct geom_parts *parts)
{
    parts->n_parts = 0;
}
 
/*!
   \brief Free parts
 */
void free_parts(struct geom_parts *parts)
{
    G_free(parts->part);
    G_zero(parts, sizeof(struct geom_parts));
}
 
/*!
   \brief Add new part number to parts
 */
void add_part(struct geom_parts *parts, int part)
{
    if (parts->a_parts == parts->n_parts) {
        parts->a_parts += 10;
        parts->part =
            (int *)G_realloc((void *)parts->part, parts->a_parts * sizeof(int));
    }
    parts->part[parts->n_parts] = part;
    parts->n_parts++;
}
 
/*!
   \brief Remove last part
 */
void del_part(struct geom_parts *parts)
{
    parts->n_parts--;
}
 
/*!
   \brief Add parts to offset
 */
void add_parts_to_offset(struct Format_info_offset *offset,
                         struct geom_parts *parts)
{
    int i, j;
 
    if (offset->array_num + parts->n_parts >= offset->array_alloc) {
        offset->array_alloc += parts->n_parts + 1000;
        offset->array =
            (int *)G_realloc(offset->array, offset->array_alloc * sizeof(int));
    }
    j = offset->array_num;
    for (i = 0; i < parts->n_parts; i++) {
        G_debug(4, "add offset %d", parts->part[i]);
        offset->array[j] = parts->part[i];
        j++;
    }
    offset->array_num += parts->n_parts;
}
 
/*!
   \brief Add line to support structures
 */
int add_line(struct Plus_head *plus, struct Format_info_offset *offset,
             int type, struct line_pnts *Points, int FID,
             struct geom_parts *parts)
{
    int line;
    long offset_value;
    struct bound_box box;
 
    if (type != GV_CENTROID) {
        /* beginning in the offset array */
        offset_value = offset->array_num;
    }
    else {
        /* TODO : could be used to statore category ? */
        /* because centroids are read from topology, not from layer */
        offset_value = FID;
    }
 
    G_debug(4, "Register line: FID = %d offset = %ld", FID, offset_value);
    dig_line_box(Points, &box);
    line = dig_add_line(plus, type, Points, &box, offset_value);
    G_debug(4, "Line registered with line = %d", line);
 
    /* Set box */
    if (line == 1)
        Vect_box_copy(&(plus->box), &box);
    else
        Vect_box_extend(&(plus->box), &box);
 
    if (type != GV_BOUNDARY) {
        dig_cidx_add_cat(plus, 1, (int)FID, line, type);
    }
    else {
        dig_cidx_add_cat(plus, 0, 0, line, type);
    }
 
    /* because centroids are read from topology, not from layer */
    if (type != GV_CENTROID)
        add_parts_to_offset(offset, parts);
 
    return line;
}
 
#ifdef HAVE_POSTGRES
/*!
   \brief Recursively add geometry (PostGIS) to topology
 */
int add_geometry_pg(struct Plus_head *plus, struct Format_info_pg *pg_info,
                    struct feat_parts *fparts, int ipart, int FID, int build,
                    struct geom_parts *parts)
{
    int line, i, idx, area, isle, outer_area, ret;
    int lines[1];
    double area_size, x, y;
    SF_FeatureType ftype;
    struct bound_box box;
    struct Format_info_offset *offset;
    struct line_pnts *line_i;
 
    ftype = fparts->ftype[ipart];
 
    G_debug(4, "add_geometry_pg() FID = %d ftype = %d", FID, ftype);
 
    offset = &(pg_info->offset);
 
    outer_area = 0;
 
    switch (ftype) {
    case SF_POINT:
        G_debug(4, "Point");
        line_i = pg_info->cache.lines[fparts->idx[ipart]];
        add_line(plus, offset, GV_POINT, line_i, FID, parts);
        break;
    case SF_LINESTRING:
        G_debug(4, "LineString");
        line_i = pg_info->cache.lines[fparts->idx[ipart]];
        add_line(plus, offset, GV_LINE, line_i, FID, parts);
        break;
    case SF_POLYGON:
        G_debug(4, "Polygon");
 
        /* register boundaries */
        idx = fparts->idx[ipart];
        for (i = 0; i < fparts->nlines[ipart]; i++) {
            line_i = pg_info->cache.lines[idx++];
            G_debug(4, "part %d", i);
            add_part(parts, i);
            line = add_line(plus, offset, GV_BOUNDARY, line_i, FID, parts);
            del_part(parts);
 
            if (build < GV_BUILD_AREAS)
                continue;
 
            /* add area (each inner ring is also area) */
            dig_line_box(line_i, &box);
            dig_find_area_poly(line_i, &area_size);
 
            if (area_size > 0) /* area clockwise */
                lines[0] = line;
            else
                lines[0] = -line;
 
            area = dig_add_area(plus, 1, lines, &box);
 
            /* each area is also isle */
            lines[0] = -lines[0]; /* island is counter clockwise */
 
            isle = dig_add_isle(plus, 1, lines, &box);
 
            if (build < GV_BUILD_ATTACH_ISLES)
                continue;
 
            if (i == 0) { /* outer ring */
                outer_area = area;
            }
            else { /* inner ring */
                struct P_isle *Isle;
 
                Isle = plus->Isle[isle];
                Isle->area = outer_area;
 
                dig_area_add_isle(plus, outer_area, isle);
            }
        }
 
        if (build >= GV_BUILD_CENTROIDS) {
            /* create virtual centroid */
            ret = Vect_get_point_in_poly_isl(
                (const struct line_pnts *)
                    pg_info->cache.lines[fparts->idx[ipart]],
                (const struct line_pnts **)&pg_info->cache
                        .lines[fparts->idx[ipart]] +
                    1,
                fparts->nlines[ipart] - 1, &x, &y);
            if (ret < -1) {
                G_warning(_("Unable to calculate centroid for area %d"),
                          outer_area);
            }
            else {
                struct P_area *Area;
                struct P_topo_c *topo;
                struct P_line *Line;
                struct line_pnts *line_c;
 
                G_debug(4, "  Centroid: %f, %f", x, y);
                line_c = Vect_new_line_struct();
                Vect_append_point(line_c, x, y, 0.0);
                line = add_line(plus, offset, GV_CENTROID, line_c, FID, parts);
 
                Line = plus->Line[line];
                topo = (struct P_topo_c *)Line->topo;
                topo->area = outer_area;
 
                /* register centroid to area */
                Area = plus->Area[outer_area];
                Area->centroid = line;
                Vect_destroy_line_struct(line_c);
            }
        }
        break;
    default:
        G_warning(_("Feature type %d not supported"), ftype);
        break;
    }
 
    return 0;
}
 
/*!
   \brief Build pseudo-topology for PostGIS layers
 */
void build_pg(struct Map_info *Map, int build)
{
    int iFeature, ipart, fid, nrecords, npoints;
    char *wkb_data;
 
    struct Format_info_pg *pg_info;
 
    struct feat_parts fparts;
    struct geom_parts parts;
 
    pg_info = &(Map->fInfo.pg);
 
    /* initialize data structures */
    init_parts(&parts);
    G_zero(&fparts, sizeof(struct feat_parts));
 
    /* get all features */
    if (Vect__open_cursor_next_line_pg(pg_info, TRUE, Map->plus.built) != 0)
        return;
 
    /* scan records */
    npoints = 0;
    nrecords = PQntuples(pg_info->res);
    G_debug(4, "build_pg(): nrecords = %d", nrecords);
    G_message(_("Registering primitives..."));
    for (iFeature = 0; iFeature < nrecords; iFeature++) {
        /* get feature id */
        fid = atoi(PQgetvalue(pg_info->res, iFeature, 1));
        if (fid < 1)
            continue; /* PostGIS Topology: skip features with negative
                       * fid (isles, universal face, ...) */
 
        wkb_data = PQgetvalue(pg_info->res, iFeature, 0);
 
        G_progress(iFeature + 1, 1e4);
 
        /* cache feature (lines) */
        if (SF_NONE == Vect__cache_feature_pg(wkb_data, FALSE, FALSE,
                                              &(pg_info->cache), &fparts)) {
            G_warning(_("Feature %d without geometry skipped"), iFeature + 1);
            continue;
        }
 
        /* register all parts */
        reset_parts(&parts);
        add_part(&parts, fid);
        for (ipart = 0; ipart < fparts.n_parts; ipart++) {
            if (fparts.nlines[ipart] < 1) {
                G_warning(_("Feature %d without geometry skipped"), fid);
                continue;
            }
 
            npoints += pg_info->cache.lines[ipart]->n_points;
 
            G_debug(4, "Feature: fid = %d part = %d", fid, ipart);
 
            if (fparts.n_parts > 1)
                add_part(&parts, ipart);
            add_geometry_pg(&(Map->plus), pg_info, &fparts, ipart, fid, build,
                            &parts);
            if (fparts.n_parts > 1)
                del_part(&parts);
        }
 
        /* read next feature from cache */
        pg_info->cache.lines_next = 0;
    }
    G_progress(1, 1);
 
    G_message(n_("One primitive registered", "%d primitives registered",
                 Map->plus.n_lines),
              Map->plus.n_lines);
    G_message(n_("One vertex registered", "%d vertices registered", npoints),
              npoints);
 
    Map->plus.built = GV_BUILD_BASE;
 
    PQclear(pg_info->res);
    pg_info->res = NULL;
 
    /* free allocated space */
    free_parts(&parts);
}
#endif /* HAVE_POSTGRES */
 
/*!
   \brief Recursively add geometry (OGR) to topology
 */
int add_geometry_ogr(struct Plus_head *plus, struct Format_info_ogr *ogr_info,
                     OGRGeometryH hGeom, int FID, int build,
                     struct geom_parts *parts)
{
    int i, ret, npoints, line;
    int area, isle, outer_area;
    int lines[1];
    double area_size, x, y;
    int eType, nRings, iPart, nParts, nPoints;
 
    struct bound_box box;
    struct P_line *Line;
    struct Format_info_offset *offset;
 
    OGRGeometryH hGeom2, hRing;
 
    G_debug(4, "add_geometry_ogr() FID = %d", FID);
 
    offset = &(ogr_info->offset);
 
    /* allocate space in cache */
    if (!ogr_info->cache.lines) {
        ogr_info->cache.lines_alloc = 1;
        ogr_info->cache.lines =
            (struct line_pnts **)G_malloc(sizeof(struct line_pnts *));
 
        ogr_info->cache.lines_types = (int *)G_malloc(sizeof(int));
        ogr_info->cache.lines[0] = Vect_new_line_struct();
        ogr_info->cache.lines_types[0] = -1;
    }
 
    npoints = outer_area = 0;
    eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
    G_debug(4, "OGR type = %d", eType);
 
    switch (eType) {
    case wkbPoint:
        G_debug(4, "Point");
 
        ogr_info->cache.lines_types[0] = GV_POINT;
        Vect_reset_line(ogr_info->cache.lines[0]);
        Vect_append_point(ogr_info->cache.lines[0], OGR_G_GetX(hGeom, 0),
                          OGR_G_GetY(hGeom, 0), OGR_G_GetZ(hGeom, 0));
        add_line(plus, offset, GV_POINT, ogr_info->cache.lines[0], FID, parts);
        npoints += ogr_info->cache.lines[0]->n_points;
        break;
 
    case wkbLineString:
        G_debug(4, "LineString");
 
        ogr_info->cache.lines_types[0] = GV_LINE;
        nPoints = OGR_G_GetPointCount(hGeom);
        Vect_reset_line(ogr_info->cache.lines[0]);
        for (i = 0; i < nPoints; i++) {
            Vect_append_point(ogr_info->cache.lines[0], OGR_G_GetX(hGeom, i),
                              OGR_G_GetY(hGeom, i), OGR_G_GetZ(hGeom, i));
        }
        add_line(plus, offset, GV_LINE, ogr_info->cache.lines[0], FID, parts);
        npoints += ogr_info->cache.lines[0]->n_points;
        break;
 
    case wkbPolygon:
        G_debug(4, "Polygon");
 
        nRings = OGR_G_GetGeometryCount(hGeom);
        G_debug(4, "Number of rings: %d", nRings);
 
        /* alloc space for islands if needed */
        if (nRings > ogr_info->cache.lines_alloc) {
            ogr_info->cache.lines_alloc += nRings;
            ogr_info->cache.lines = (struct line_pnts **)G_realloc(
                ogr_info->cache.lines,
                ogr_info->cache.lines_alloc * sizeof(struct line_pnts *));
            ogr_info->cache.lines_types =
                (int *)G_realloc(ogr_info->cache.lines_types,
                                 ogr_info->cache.lines_alloc * sizeof(int));
 
            for (i = ogr_info->cache.lines_alloc - nRings;
                 i < ogr_info->cache.lines_alloc; i++) {
                ogr_info->cache.lines[i] = Vect_new_line_struct();
                ogr_info->cache.lines_types[i] = -1;
            }
        }
 
        /* go through rings */
        for (iPart = 0; iPart < nRings; iPart++) {
            ogr_info->cache.lines_types[iPart] = GV_BOUNDARY;
            hRing = OGR_G_GetGeometryRef(hGeom, iPart);
            nPoints = OGR_G_GetPointCount(hRing);
            G_debug(4, "  ring %d : nPoints = %d", iPart, nPoints);
 
            Vect_reset_line(ogr_info->cache.lines[iPart]);
            for (i = 0; i < nPoints; i++) {
                Vect_append_point(ogr_info->cache.lines[iPart],
                                  OGR_G_GetX(hRing, i), OGR_G_GetY(hRing, i),
                                  OGR_G_GetZ(hRing, i));
            }
            npoints += ogr_info->cache.lines[iPart]->n_points;
 
            /* register boundary */
            add_part(parts, iPart);
            line = add_line(plus, offset, GV_BOUNDARY,
                            ogr_info->cache.lines[iPart], FID, parts);
            del_part(parts);
 
            if (build < GV_BUILD_AREAS)
                continue;
 
            /* add area (each inner ring is also area) */
            dig_line_box(ogr_info->cache.lines[iPart], &box);
            dig_find_area_poly(ogr_info->cache.lines[iPart], &area_size);
 
            if (area_size > 0) /* area clockwise */
                lines[0] = line;
            else
                lines[0] = -line;
 
            area = dig_add_area(plus, 1, lines, &box);
 
            /* each area is also isle */
            lines[0] = -lines[0]; /* island is counter clockwise */
 
            isle = dig_add_isle(plus, 1, lines, &box);
 
            if (build < GV_BUILD_ATTACH_ISLES)
                continue;
 
            if (iPart == 0) { /* outer ring */
                outer_area = area;
            }
            else { /* inner ring */
                struct P_isle *Isle;
 
                Isle = plus->Isle[isle];
                Isle->area = outer_area;
 
                dig_area_add_isle(plus, outer_area, isle);
            }
        }
 
        if (build >= GV_BUILD_CENTROIDS) {
            /* create virtual centroid */
            ret = Vect_get_point_in_poly_isl(
                (const struct line_pnts *)ogr_info->cache.lines[0],
                (const struct line_pnts **)ogr_info->cache.lines + 1,
                nRings - 1, &x, &y);
            if (ret < -1) {
                G_warning(_("Unable to calculate centroid for area %d"),
                          outer_area);
            }
            else {
                struct P_area *Area;
                struct P_topo_c *topo;
 
                G_debug(4, "  Centroid: %f, %f", x, y);
                Vect_reset_line(ogr_info->cache.lines[0]);
                Vect_append_point(ogr_info->cache.lines[0], x, y, 0.0);
                line = add_line(plus, offset, GV_CENTROID,
                                ogr_info->cache.lines[0], FID, parts);
 
                Line = plus->Line[line];
                topo = (struct P_topo_c *)Line->topo;
                topo->area = outer_area;
 
                /* register centroid to area */
                Area = plus->Area[outer_area];
                Area->centroid = line;
            }
        }
        break;
 
    case wkbMultiPoint:
    case wkbMultiLineString:
    case wkbMultiPolygon:
    case wkbGeometryCollection:
        nParts = OGR_G_GetGeometryCount(hGeom);
        G_debug(4, "%d geoms -> next level", nParts);
 
        /* alloc space for parts if needed */
        if (nParts > ogr_info->cache.lines_alloc) {
            ogr_info->cache.lines_alloc += nParts;
            ogr_info->cache.lines = (struct line_pnts **)G_realloc(
                ogr_info->cache.lines,
                ogr_info->cache.lines_alloc * sizeof(struct line_pnts *));
            ogr_info->cache.lines_types =
                (int *)G_realloc(ogr_info->cache.lines_types,
                                 ogr_info->cache.lines_alloc * sizeof(int));
 
            for (i = ogr_info->cache.lines_alloc - nParts;
                 i < ogr_info->cache.lines_alloc; i++) {
                ogr_info->cache.lines[i] = Vect_new_line_struct();
                ogr_info->cache.lines_types[i] = -1;
            }
        }
 
        /* go through all parts */
        for (i = 0; i < nParts; i++) {
            add_part(parts, i);
            hGeom2 = OGR_G_GetGeometryRef(hGeom, i);
            npoints +=
                add_geometry_ogr(plus, ogr_info, hGeom2, FID, build, parts);
            del_part(parts);
        }
        break;
 
    default:
        G_warning(_("OGR feature type %d not supported"), eType);
        break;
    }
 
    return npoints;
}
 
void build_ogr(struct Map_info *Map, int build)
{
    int iFeature, FID, npoints, nskipped;
 
    struct Format_info_ogr *ogr_info;
 
    OGRFeatureH hFeature;
    OGRGeometryH hGeom;
 
    struct geom_parts parts;
 
    ogr_info = &(Map->fInfo.ogr);
 
    /* initialize data structures */
    init_parts(&parts);
 
    /* Note: Do not use OGR_L_GetFeatureCount (it may scan all features) */
    OGR_L_ResetReading(ogr_info->layer);
    if (ogr_info->where)
        /* set attribute filter if where sql statement defined */
        OGR_L_SetAttributeFilter(ogr_info->layer, ogr_info->where);
    npoints = iFeature = nskipped = 0;
    G_message(_("Registering primitives..."));
    while ((hFeature = OGR_L_GetNextFeature(ogr_info->layer)) != NULL) {
        G_debug(3, "   Feature %d", iFeature);
 
        G_progress(++iFeature, 1e4);
 
        hGeom = OGR_F_GetGeometryRef(hFeature);
        if (hGeom == NULL) {
            G_debug(3, "Feature %d without geometry skipped", iFeature);
            OGR_F_Destroy(hFeature);
            nskipped++;
            continue;
        }
 
        FID = (int)OGR_F_GetFID(hFeature);
        if (FID == OGRNullFID) {
            G_debug(3, "OGR feature %d without ID skipped", iFeature);
            OGR_F_Destroy(hFeature);
            nskipped++;
            continue;
        }
        G_debug(4, "    FID = %d", FID);
 
        reset_parts(&parts);
        add_part(&parts, FID);
        npoints +=
            add_geometry_ogr(&(Map->plus), ogr_info, hGeom, FID, build, &parts);
 
        OGR_F_Destroy(hFeature);
    } /* while */
    G_progress(1, 1);
 
    G_message(n_("One primitive registered", "%d primitives registered",
                 Map->plus.n_lines),
              Map->plus.n_lines);
    G_message(n_("One vertex registered", "%d vertices registered", npoints),
              npoints);
 
    if (nskipped > 0)
        G_warning(n_("One feature without geometry skipped",
                     "%d features without geometry skipped", nskipped),
                  nskipped);
 
    Map->plus.built = GV_BUILD_BASE;
 
    free_parts(&parts);
}
 
/*!
   \brief Build pseudo-topology (for simple features) - internal use only
 
   See Vect_build_ogr() and Vect_build_pg() for implementation issues.
 
   Build levels:
   - GV_BUILD_NONE
   - GV_BUILD_BASE
   - GV_BUILD_ATTACH_ISLES
   - GV_BUILD_CENTROIDS
   - GV_BUILD_ALL
 
   \param Map pointer to Map_info structure
   \param build build level
 
   \return 1 on success
   \return 0 on error
 */
int Vect__build_sfa(struct Map_info *Map, int build)
{
    struct Plus_head *plus;
 
    plus = &(Map->plus);
 
    /* check if upgrade or downgrade */
    if (build < plus->built) {
        /* -> downgrade */
        Vect__build_downgrade(Map, build);
        return 1;
    }
 
    /* -> upgrade */
    if (plus->built < GV_BUILD_BASE) {
        if (Map->format == GV_FORMAT_OGR ||
            Map->format == GV_FORMAT_OGR_DIRECT) {
            build_ogr(Map, build);
        }
        else if (Map->format == GV_FORMAT_POSTGIS) {
#ifdef HAVE_POSTGRES
            build_pg(Map, build);
#else
            G_fatal_error(_("GRASS is not compiled with PostgreSQL support"));
#endif
        }
        else {
            G_fatal_error(_("%s: Native format unsupported"),
                          "Vect__build_sfa()");
        }
    }
 
    plus->built = build;
 
    return 1;
}
 
/*!
   \brief Dump feature index to file
 
   \param Map pointer to Map_info struct
   \param out file for output (stdout/stderr for example)
 
   \return 1 on success
   \return 0 on error
 */
int Vect_fidx_dump(struct Map_info *Map, FILE *out)
{
    int i;
    const struct Format_info_offset *offset;
 
    if (Map->format != GV_FORMAT_OGR && Map->format != GV_FORMAT_POSTGIS) {
        G_warning(_("Feature index is built only for non-native formats. "
                    "Nothing to dump."));
        return 0;
    }
 
    if (Map->format == GV_FORMAT_OGR)
        offset = &(Map->fInfo.ogr.offset);
    else
        offset = &(Map->fInfo.pg.offset);
 
    fprintf(out, "---------- FEATURE INDEX DUMP ----------\n");
 
    fprintf(out, "format: %s\n", Vect_maptype_info(Map));
    if (Vect_maptype(Map) == GV_FORMAT_POSTGIS && Map->fInfo.pg.toposchema_name)
        fprintf(out, "topology: PostGIS\n");
    else
        fprintf(out, "topology: pseudo\n");
    fprintf(out, "feature type: %s\n", Vect_get_finfo_geometry_type(Map));
    fprintf(out,
            "number of features: %d\n\noffset : value (fid or part idx):\n",
            Vect_get_num_lines(Map));
    for (i = 0; i < offset->array_num; i++) {
        fprintf(out, "%6d : %d\n", i, offset->array[i]);
    }
 
    return 1;
}